D:/demona/2.0/manual/doc/2-31set/processor/newton/newton1set.f90

Go to the documentation of this file.

subroutine newtonSet

use global

use model

use krylov

use demonaParameters

implicit none

real(8),allocatable :: deltax(:),F0(:),resVec(:),setF0(:)

real(8) normF,normFrms

integer newtonIter


real(8) newtonTolerance

real(8) inexactLinearTolerance

integer linearIterations

integer i

!<--------------------------------------------
!<crs dependent declarations
integer indSize
integer ptrSize
integer valSize,valSizeMax

real(8), allocatable :: JacVal(:),luval(:)
integer, allocatable :: JacColInd(:),JacRowPtr(:),iw(:),uptr(:)

integer numberofcolors






!######################################################

newtonTolerance=1e-5
inexactLinearTolerance=1e-5


newtonSuccess=0

Open (2,File='newtonResidual.dat',Status='unknown')
!print*,setvec

!dof=>globalcurrentsetlength



!allocate(deltax(globalcurrentsetlength))
allocate(F0(globalcurrentsetlength))
!allocate(resVec(globalcurrentsetlength))
allocate(setF0(globalcurrentsetlength))
!deltax=0.d0




F0=discrete(0)
setF0=discrete(0)

!print*,F0(69),y(69)


!stop

normF=dsqrt(dot_product(F0,F0))

print*,'initial norm:',normF

write(2,'(F20.12,1X)')normF


if (normF.le.newtonTolerance) then

        print*,'zero initial norm'
        print*,'newton return'
        newtonSuccess=1

        return


end if


call symbolic(valsizeMax)
ptrSize=globalcurrentsetlength+1


totalLinearIterations=0

!print*,discrete(0)

tempy=y

!print*,y
!pause


do newtonIter=1,maxNewtonIter


        if (allocated(deltax)) deallocate(deltax)
        if (allocated(resVec)) deallocate(resVec)
                
        allocate(deltax(globalcurrentsetlength))
        allocate(resVec(globalcurrentsetlength))
        !allocate(setF0(globalcurrentsetlength))

        deltax=0.d0
        resVec=0.d0
        !setF0=discrete(1)

        print*,'newton iter and dof=',newtonIter,globalcurrentsetlength


        !inexactLinearTolerance=1e-5
        inexactLinearTolerance=newtonTolerance
        inexactLinearTolerance=1e-1*normF

        deltax=0.d0


                        ! compute the jacobian
                        if (allocated(JacVal)) deallocate(JacVal)
                        if (allocated(JacColInd)) deallocate(JacColInd)
                        if (allocated(uptr)) deallocate(uptr)
                        if (allocated(JacRowPtr)) deallocate(JacRowPtr)

                        !discOrder=lowOrder

                        allocate(JacVal(valSizeMax),JacColInd(valSizeMax),JacRowPtr(ptrSize),uptr(globalcurrentsetlength))
                        
                        !symbolic=0
                        JacVal=0.d0
                        JacColInd=0
                        uptr=0
                        JacRowPtr=0

                        tempy=y

                        call colorJacobianCRS(JacVal,JacColInd,JacRowPtr,globalcurrentsetlength,1,valSize,indSize,uptr,0,numberOfColors)

                        
                        
                        !if (allocated(luval)) deallocate(luval)
                        !allocate(luVal(valSize))
                        !luval=0.d0

                        !luval=JacVal(1:valsize)

                        ! compute ilu(0)
                        !call basicilu0CRS(luval,JacColInd(1:valsize),JacRowPtr,globalcurrentsetlength,valSize,uptr)
                        call basicilu0CRS(jacval(1:valsize),JacColInd(1:valsize),JacRowPtr,globalcurrentsetlength,valSize,uptr)

                        !luval=JacVal(1:valsize)

                        !discOrder=myModel%physics(currentPhysics)%discOrder

                        JacVal(1:valsize)=JacVal(1:valsize)+epsilon(0.d0)
                        JacVal(1:valsize)=JacVal(1:valsize)-epsilon(0.d0)






        !print*,size(deltax),size(f0)
        !stop


        !call bicgstabRP(-F0(setVec),deltax,inexactLinearTolerance,1000,resVec,luval,JacColInd,JacRowPtr,uptr,valsize,LinearIterations)
        !call gmres(-setF0,deltax,500,2,inexactLinearTolerance,resVec,luval,jaccolInd,jacrowPtr,uptr,valsize,LinearIterations)

        !call gmres(-F0(setVec),deltax,500,3,inexactLinearTolerance,resVec,luval,jaccolInd,jacrowPtr,uptr,valsize,LinearIterations)
        call gmres(-F0,deltax,200,20,inexactLinearTolerance,resVec,JacVal(1:valsize),jaccolInd(1:valsize),jacrowPtr,uptr,valsize,LinearIterations,2)
        
        tempy(setVec)=tempy(setVec)+deltax
        !y(setVec)=y(setVec)+deltax

        !print*,deltax

        !myModel%physics(currentPhysics)%domain(currentDomain)%dofMatrix(:,:,:,flag)=1

        !call currentSet


        !F0(setVec)=discrete(1)  ! not setvec but some  relative index should be used




        
        F0=discrete(1)  ! not setvec but some  relative index should be used

        !if (newtoniter.eq.2) then
        !       do i=1,17*17
        !               print*,F0(i),tempy(i)
        !       end do
        !       stop
        !end if
        !print*,'------------'
        !pause

        normF=dsqrt(dot_product(F0,F0))

        normFrms=normF/dsqrt(size(F0)*1.d0)
        print*,'norm=',normF
        print*,'rms norm=',normFrms


        totalLinearIterations=totalLinearIterations+LinearIterations
        

        y=tempy
        write(2,'(F20.12,1X)')normF

        
        !tempy(1:33*33)=F0
        !call tecFD
        !pause
        !tempy=y



        if (normF.le.newtonTolerance) then
        !if (normFrms.le.newtonTolerance) then

                
                print*,'newton converged at iter',newtonIter
                print*,'residual:',normF
                print*,'total linear iterations:',totalLinearIterations

                        !tempy(1:33*33)=F0
                        !call tecFD
        !               pause
                        !tempy=y

                !y=tempy

                exit


        end if

        



        
        if (newtonIter.eq.7) then


        !tempy(1:9*9)=F0

        !call tecFD

        !stop

        end if
        
        
        ! if step is accepted
        




        
        !print*,y(67)


        if ((resSet.eqv.yes).and.(newtonIter.gt.0)) then

                
                !call fullSet !< else the set keep reducing with no improvement since zeroed points cannot be turned into 1

                !myModel%physics(currentPhysics)%domain(currentDomain)%dofMatrix(:,:,:,flag)=1

                !call currentSet

                !call setRule(F0(setVec),setVec,normF/size(F0))
                !call setRule(F0,setVec,normF/dsqrt(size(F0)*1.d0))
                call setRule(F0,size(F0),normF,normFrms)

                call currentSet


                !dof=>globalcurrentsetlength
                
                !deallocate(deltax,resvec,setF0)
                deallocate(deltax,resvec)
                

                !allocate(deltax(globalcurrentsetlength))
                !allocate(resVec(globalcurrentsetlength))
                !allocate(setF0(globalcurrentsetlength))

                !deltax=0.d0
                !resVec=0.d0
                !setF0=0.d0

                !setF0=discrete(1)

                !print*,globalcurrentsetlength !,setVEc(1:70)
                !print*,setVEc

                !tempy(1:5*5)=F0
                !call tecFD
                !pause
                !tempy=y

                        !do i=1,9
                        !       print*,setVec(i),setF0(i) !,tempy(i)
                        !end do
                        !stop

                !print*,size(setVec)

                !stop




        end if








end do

close(2)



        !Open (4400,File='x.dat',Status='unknown')  !<--- physics and domain and subdomain based numbering
        
        
        
        910 FORMAT((E18.10))


        !WRITE(4400,910)y
        !close(4400)





end subroutine newtonSet

---