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

Go to the documentation of this file.

subroutine continuation

use global

use krylov

use mf



implicit none


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

integer valsize,valsizeMax,ptrSize,indSize,numberOfColors

real(8),allocatable :: jacVal(:)
integer,allocatable :: jaccolInd(:),jacrowPtr(:),uptr(:)

real(8) inexactLinearTolerance,initial_norm

real(8) backup,myEps


integer :: LinearIterations

allocate(delta(globalcurrentsetlength))
allocate(F0(globalcurrentsetlength))

allocate(resVec(globalcurrentsetlength))



if (bootStrapping.eqv.on) return


print*,'"""""""""""""""""""""""""""""""""""""'
print*,'continuation'
print*,'"""""""""""""""""""""""""""""""""""""'


F0=0.d0
delta=0.d0
resVEc=0.d0



myEps=dsqrt(epsilon(0.d0))
!myEps=epsilon(0.0)


inexactLinearTolerance=0.00001d0


myModel%physics(1)%physicalParameters(1)=myModel%physics(1)%physicalParametersInitial(1)


F0=discrete(0)



myModel%physics(1)%physicalParameters(1)=myModel%physics(1)%physicalParameters(1)+myEps

!same y can be used since y is not modified, just the parameter is perturbed

F0=(discrete(0)-F0)/myEps ! delF/delPi

myModel%physics(1)%physicalParameters(1)=myModel%physics(1)%physicalParametersInitial(1)


call symbolic(valsizeMax)
ptrSize=globalcurrentsetlength+1


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



initial_norm=1.d0

delta=-F0-matvec(delta)


initial_norm=dsqrt(dot_product(delta,delta))
!initial_norm=dsqrt(dot_product(F0,F0))

!print*,initial_norm
!stop


!inexactLinearTolerance=inexactLinearTolerance*initial_norm


!print*,inexactLinearTolerance

!print*,f0(133:134)
!stop


!print*,f0(67),dexp(y(67)),myModel%physics(1)%physicalParametersFinal(1)-myModel%physics(1)%physicalParametersInitial(1)
!stop


delta=0.d0


call gmres(-F0,delta,200,20,inexactLinearTolerance,resVec,JacVal(1:valsize),jaccolInd(1:valsize),jacrowPtr,uptr,valsize,LinearIterations,2)

!print*,delta(67)
!stop

!stop
tempy=y


y(setVec)=y(setvec)+bootdamping*delta*(myModel%physics(1)%physicalParametersFinal(1)-myModel%physics(1)%physicalParametersInitial(1))


!print*,y(67)
!stop


myModel%physics(1)%physicalParameters(1)=myModel%physics(1)%physicalParametersFinal(1)


end subroutine continuation

---