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

Go to the documentation of this file.

subroutine newton

use global

use model

use krylov

use demonaParameters

use DFPORT

implicit none

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

real(8) normF

integer newtonIter


real(8) newtonTolerance

real(8) inexactLinearTolerance

integer linearIterations

real(8) newtonConvergenceTolerance



!inexact newton declarations

real(8) eta


!jacobian related declarations

integer myJacInterval


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

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


integer numberOfColors

!<------------- LSBT

real(8) lambda



!<-------------- IN


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


myJacInterval=jacInterval !< this way the original jacInterval is protected

if (freezedJacobian.eqv.on) then

        
        myJacInterval=maxNewtonIter !< so the jacobian is calculated only at the first iteration 

        if (gsSolver.eqv.no) print*,'freezed jacobian'

end if



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


eta=1e-2

!print*,setvec

!dof=>globalcurrentsetlength



allocate(deltax(globalcurrentsetlength))
allocate(F0(globalcurrentsetlength))
allocate(resVec(globalcurrentsetlength))

deltax=0.d0




F0=discrete(0)

totalLinearIterations=0

!print*,f0
!stop
!pause



normF=dsqrt(dot_product(F0,F0))



newtonConvergenceTolerance=newtonTolerance















if (gsSolver.eqv.no) print*,'initial norm:',normF


if (normF.eq.0.d0) then

        if (gsSolver.eqv.no) print*,'zero initial norm'
        if (gsSolver.eqv.no) print*,'newton return'
        return


end if

!print*,f0
!stop

!print*,globalcurrentsetlength
!stop


!< symbolic jacobian evaluation do detect the size of the crs-jacobian
call symbolic(valsizeMax)
ptrSize=globalcurrentsetlength+1



do newtonIter=1,maxNewtonIter


                if (mod(newtonIter+myJacInterval-1,myJacInterval).eq.0) then

                        if (normF.le.JacSkipTol*newtonConvergenceTolerance) then

                                if(jacSkip.eqv.on) then
                                        print*,'almost the same jacobian probably'
                                        print*,'its evaluation is skipped.'
                                        goto 1
                                end if

                        end if

                        ! 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)

                        tempy=y
                        
                        !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,JacColInd(1:valsize),JacRowPtr,globalcurrentsetlength,valSize,uptr)

                        !luval=JacVal(1:valsize)

                        !discOrder=myModel%physics(currentPhysics)%discOrder

                end if

                1 continue
                

        if (gsSolver.eqv.no) print*,'newton iter=',newtonIter

        inexactLinearTolerance=1e-5

        if (inexactNewton.eqv.On) inexactLinearTolerance=eta*normF

        deltax=0.d0

        !call bicgstabRP(-F0,deltax,inexactLinearTolerance,1000,resVec,luval,JacColInd,JacRowPtr,uptr,valsize,LinearIterations)
        call gmres(-F0,deltax,500,20,inexactLinearTolerance,resVec,JacVal(1:valsize),jaccolInd(1:valsize),jacrowPtr,uptr,valsize,LinearIterations,2)
        
        
        lambda=1.d0
        !lambda=0.5d0
        
        tempy(setVec)=tempy(setVec)+lambda*deltax
        !y(setVec)=y(setVec)+deltax

        F0=discrete(1)

        normF=dsqrt(dot_product(F0,F0))
        
        totalLinearIterations=totalLinearIterations+LinearIterations

        if (gsSolver.eqv.no) print*,normF
        
        if (normF.le.newtonConvergenceTolerance) then

                
                if (gsSolver.eqv.no) print*,'newton converged at iter',newtonIter
                if (gsSolver.eqv.no) print*,'residual:',normF
                if (gsSolver.eqv.no) print*,'total linear iterations:',totalLinearIterations

                y=tempy

                exit


        end if

        
        ! if step is accepted
        
        y=tempy


        if (resSet.eqv.yes) then

                call setRule

                call set


        end if


end do



end subroutine newton

---