demona: D:/demona/2.0/manual/doc/2-31set/processor/newton/newton2.f90 Source File

00001 subroutine newton
00002 
00003 use global
00004 
00005 use model
00006 
00007 use krylov
00008 
00009 use demonaParameters
00010 
00011 use DFPORT
00012 
00013 implicit none
00014 
00015 real(8),allocatable :: deltax(:),F0(:),resVec(:),F2(:),wM(:)
00016 
00017 real(8) stepWNorm
00018 
00019 integer newtonIter,stepLengthConverged
00020 
00021 integer j
00022 
00023 integer ii,jj
00024 
00025 
00026 real(8) newtonTolerance
00027 
00028 real(8) inexactLinearTolerance
00029 
00030 integer linearIterations
00031 
00032 real(8) newtonConvergenceTolerance
00033 
00034 
00035 
00036 !inexact newton declarations
00037 
00038 real(8) eta
00039 
00040 
00041 !jacobian related declarations
00042 
00043 integer myJacInterval
00044 
00045 
00046 !<--------------------------------------------
00047 !<crs dependent declarations
00048 integer indSize
00049 integer ptrSize
00050 integer valSize,valSizeMax
00051 
00052 real(8), allocatable :: JacVal(:),luval(:)
00053 integer, allocatable :: JacColInd(:),JacRowPtr(:),iw(:),uptr(:)
00054 
00055 
00056 integer numberOfColors
00057 
00058 !<------------- LSBT
00059 
00060 real(8) lambda,lambda2
00061 
00062 real(8) fZero, fTwo
00063 
00064 real(8) slope
00065 
00066 real(8) thetaINMin,thetaINMax
00067 
00068 !<-------------- IN
00069 
00070 real(8) etaZero,etaPrev
00071 
00072 real(8) etaVec(maxNewtonIter)
00073 
00074 real(8) etaMax
00075 
00076 real(8) normFcurrent,normNewFcurrent,normFprev,normRes
00077 
00078 real(8) initialNorm
00079 
00080 real(8) testLHS,testRHS
00081 
00082 real(8) tIn,alphaIn,gammaIN,compareIN
00083 
00084 integer choiceIn
00085 
00086 !######################################################
00087 
00088 Open (2,File='newtonResidual.dat',Status='unknown')
00089 
00090 
00091 
00092 myJacInterval=jacInterval !< this way the original jacInterval is protected
00093 
00094 if (freezedJacobian.eqv.on) then
00095 
00096         
00097         myJacInterval=maxNewtonIter !< so the jacobian is calculated only at the first iteration 
00098 
00099         if (gsSolver.eqv.no) print*,'freezed jacobian'
00100 
00101 end if
00102 
00103 
00104 
00105 newtonTolerance=1e-5
00106 newtonTolerance=0.01d0
00107 inexactLinearTolerance=1e-5
00108 
00109 
00110 eta=0.01d0
00111 etaZero=0.01d0
00112 
00113 tIN=0.0001d0
00114 
00115 thetaINMin=0.1d0
00116 thetaINMax=0.5d0
00117 
00118 alphaIN=1.5d0
00119 alphaIN=2.d0
00120 
00121 gammaIN=0.9d0
00122 gammaIN=1.d0
00123 
00124 etaMax=0.9d0
00125 
00126 choiceIN=1
00127 
00128 newtonSuccess=0
00129 
00130 !print*,setvec
00131 
00132 !dof=>globalcurrentsetlength
00133 
00134 
00135 
00136 allocate(deltax(globalcurrentsetlength))
00137 allocate(F0(globalcurrentsetlength))
00138 allocate(F2(globalcurrentsetlength))
00139 allocate(resVec(globalcurrentsetlength))
00140 allocate(wM(globalcurrentsetlength))
00141 
00142 deltax=0.d0
00143 
00144 
00145 
00146 
00147 F0=discrete(0)
00148 
00149 totalLinearIterations=0
00150 
00151 !print*,f0
00152 !stop
00153 !pause
00154 
00155 
00156 
00157 normNewFcurrent=dsqrt(dot_product(F0,F0))
00158 initialNorm=normNewFcurrent
00159 write(2,'(F20.12,1X)')normNewFcurrent
00160 
00161 newtonConvergenceTolerance=newtonTolerance
00162 newtonConvergenceTolerance=newtonTolerance*initialNorm
00163 
00164 inexactLinearTolerance=newtonConvergenceTolerance
00165 
00166 
00167 
00168 
00169 if (gsSolver.eqv.no) print*,'initial norm:',normNewFcurrent
00170 
00171 !stop
00172 if (normNewFcurrent.le.newtonConvergenceTolerance) then
00173 
00174         if (gsSolver.eqv.no) print*,'zero initial norm'
00175         if (gsSolver.eqv.no) print*,'newton return'
00176         return
00177 
00178 
00179 end if
00180 
00181 !print*,f0
00182 !stop
00183 
00184 !print*,globalcurrentsetlength
00185 !stop
00186 
00187 
00188 !< symbolic jacobian evaluation do detect the size of the crs-jacobian
00189 call symbolic(valsizeMax)
00190 ptrSize=globalcurrentsetlength+1
00191 
00192 eta=etaZero
00193 
00194 
00195 !print*,'etaprev!!!!!!!!!!!!!!!!!!!'
00196 !stop
00197 
00198 
00199 etaPrev=eta
00200 
00201 
00202 
00203 
00204 do newtonIter=1,maxNewtonIter
00205 
00206         !if (newtonIter.eq.2) test=1
00207 
00208         !print*,newtoniter
00209 
00210 
00211         etaVec(newtonIter)=eta
00212 
00213         if (inexactNewton.eqv.On) inexactLinearTolerance= normNewFcurrent*eta
00214 
00215 
00216                 if (mod(newtonIter+myJacInterval-1,myJacInterval).eq.0) then
00217 
00218                         if (normFcurrent.le.JacSkipTol*newtonConvergenceTolerance) then
00219 
00220                                 if(jacSkip.eqv.on) then
00221                                         print*,'almost the same jacobian probably'
00222                                         print*,'its evaluation is skipped.'
00223                                         goto 1
00224                                 end if
00225 
00226                         end if
00227 
00228                         ! compute the jacobian
00229                         if (allocated(JacVal)) deallocate(JacVal)
00230                         if (allocated(JacColInd)) deallocate(JacColInd)
00231                         if (allocated(uptr)) deallocate(uptr)
00232                         if (allocated(JacRowPtr)) deallocate(JacRowPtr)
00233 
00234                         !discOrder=lowOrder
00235 
00236                         allocate(JacVal(valSizeMax),JacColInd(valSizeMax),JacRowPtr(ptrSize),uptr(globalcurrentsetlength))
00237                         
00238                         !symbolic=0
00239                         JacVal=0.d0
00240                         JacColInd=0
00241                         uptr=0
00242                         JacRowPtr=0
00243 
00244                         tempy=y
00245 
00246                         myModel%physics(currentPhysics)%evalOrder=>jacDiscOrder
00247 
00248                         call colorJacobianCRS(JacVal,JacColInd,JacRowPtr,globalcurrentsetlength,1,valSize,indSize,uptr,0,numberOfColors)
00249 
00250                         myModel%physics(currentPhysics)%evalOrder=>myModel%physics(currentPhysics)%discOrder
00251 
00252                         !open(9,file='row.dat',status='unknown')
00253                         !write(9,'(F20.12,1X)')JacVal(JacRowPtr(37):JacRowPtr(38)-1)
00254                         !close(9)
00255                         !if (newtoniter.eq.2) stop      
00256                         
00257 
00258                         !if (allocated(luval)) deallocate(luval)
00259                         !allocate(luVal(valSize))
00260                         !luval=0.d0
00261 
00262                         !luval=JacVal(1:valsize)
00263 
00264 
00265                         wM=0.d0
00266 
00267                         do jj=1,globalcurrentsetlength
00268 
00269                                 do ii=jacrowPtr(jj),jacrowPtr(jj+1)-1
00270                                         
00271                                         wM(jj)=dabs(jacVal(ii))+resvec(jj) ! absolute row sum of the jacobian 
00272 
00273 
00274                                 end do
00275 
00276 
00277                         end do
00278 
00279 
00280 
00281                         ! compute ilu(0)
00282                         !call basicilu0CRS(luval,JacColInd(1:valsize),JacRowPtr,globalcurrentsetlength,valSize,uptr)
00283                         call basicilu0CRS(jacval(1:valsize),JacColInd(1:valsize),JacRowPtr,globalcurrentsetlength,valSize,uptr)
00284 
00285                         !luval=JacVal(1:valsize)
00286 
00287                         !discOrder=myModel%physics(currentPhysics)%discOrder
00288 
00289                         JacVal(1:valsize)=JacVal(1:valsize)+epsilon(0.d0)
00290                         JacVal(1:valsize)=JacVal(1:valsize)-epsilon(0.d0)
00291 
00292                 end if
00293 
00294                 1 continue
00295 
00296                 !print*,'==============',valsize
00297                 !pause
00298                 
00299         tempy=y
00300 
00301 
00302 
00303 
00304         if (gsSolver.eqv.no) print*,'newton iter=',newtonIter
00305 
00306         !inexactLinearTolerance=1e-5
00307 
00308         !if (gsSolver.eqv.no) print*,'linear',inexactLinearTolerance
00309 
00310         
00311         !eta=1e-2
00312         
00313         !inexactLinearTolerance= normNewFcurrent*eta
00314 
00315 
00316         !print*,inexactLinearTolerance,normNewFcurrent,eta
00317         !pause
00318 
00319         deltax=0.d0
00320 
00321         !inexactLinearTolerance=eta
00322 
00323 
00324         !wM=jacval(uptr) !should be the diagonals
00325 
00326         wM=1.d0
00327 
00328         resVec=wM
00329 
00330         !print*,maxval(setvec)
00331         !stop
00332 
00333         !call bicgstabRP(-F0,deltax,inexactLinearTolerance,1000,resVec,luval,JacColInd,JacRowPtr,uptr,valsize,LinearIterations)
00334         call gmres(-F0,deltax,500,40,inexactLinearTolerance,resVec,JacVal(1:valsize),jaccolInd(1:valsize),jacrowPtr,uptr,valsize,LinearIterations,preconditioner)
00335         
00336         
00337         !stop
00338         !wM=1.d0
00339         
00340         !print*,tempy(133),y(133),deltax(133)
00341 
00342         !totalLinearIterations=totalLinearIterations+LinearIterations
00343 
00344     lambda=1.d0
00345     lambda2=lambda
00346         
00347         F0=discrete(0)
00348 
00349         resVec=-F0-matVec(deltax)
00350 
00351         normRes=dsqrt(dot_product(resVec,resVec))
00352         !normRes=wNorm(resVec,wM)
00353 
00354 
00355         tempy(setVec)=tempy(setVec)+lambda*deltax
00356         !y(setVec)=y(setVec)+deltax
00357 
00358         !print*,tempy(setVec(133)) !+lambda*deltax(133)
00359         !stop
00360 
00361 
00362         F2=discrete(1)
00363 
00364 
00365         !print*,tempy(133),y(133),deltax(133)
00366         !stop
00367 
00368 
00369         fZero=0.5d0*dot_product(F0,F0)
00370         fzero=0.5d0*dot_product(invscaler(F0,wM),invScaler(F0,wM))
00371     
00372         fTwo=0.5d0*dot_product(F2,F2)
00373         fTwo=0.5d0*dot_product(invscaler(F2,wM),invScaler(F2,wM))
00374 
00375         slope=-dot_product(F0,F0) -dot_product(F0,resVec)
00376         slope=-dot_product(invscaler(F0,wM),invScaler(F0,wM)) -dot_product(invscaler(F0,wM),invScaler(resVec,wM))
00377 
00378         normFcurrent=dsqrt(dot_product(F0,F0))
00379         !normFcurrent=wNorm(F0,wM)
00380 
00381         normNewFcurrent=dsqrt(dot_product(F2,F2))
00382         !normNewFcurrent=wNorm(F2,wM)
00383         
00384         totalLinearIterations=totalLinearIterations+LinearIterations
00385 
00386 
00387         if (doBacktrack.eqv.on) then
00388                 if (slope.ge.0.d0) then
00389 
00390                         !print*,'positive slope'
00391                         !print*,'line search fails'
00392                         
00393                         !print*,slope,dot_product(F0,F0),dot_product(F0,resVec)
00394                         !stop
00395                         !exit
00396 
00397                 end if
00398         end if
00399         
00400         !print*,normNewFcurrent,normFcurrent
00401         !stop
00402 
00403         
00404         if (doBacktrack.eqv.off) then
00405 
00406                 y(setVec)=y(setVec)+deltax
00407                 !print*,normNewFcurrent,normFcurrent
00408                 !stop
00409 
00410         else
00411 
00412 backtrack:      do j=1,maxBacktrackCount
00413                                 
00414 
00415                                 
00416                         testLHS=normNewFcurrent
00417                         testRHS=(1.d0-tIn*(1.d0-eta))*normFcurrent
00418 
00419                         !print*,testLHS,testRHS
00420                         !stop
00421 
00422                         if (testLHS.le.testRHS) then
00423 
00424                                 normFprev=normFcurrent
00425                                 F0=F0+matVec(deltaX) ! use F0 as temporary vector so dot_product is simple
00426                                 
00427                                 normRes=dsqrt(dot_product(F0,F0))
00428                                 !normRes=wNorm(F0,wM)
00429 
00430                                 y(setVec)=y(setVec)+deltax
00431 
00432                                 F2=discrete(0)
00433 
00434                                 normNewFcurrent=dsqrt(dot_product(F2,F2))
00435                                 !normNewFcurrent=wNorm(F2,wM)
00436 
00437 
00438                                 !print*,normNewFcurrent,normFcurrent
00439                                 !stop
00440 
00441                                 exit
00442 
00443                         else
00444 
00445 
00446                         if (j.eq.1)  then
00447                                 
00448                                 lambda=quadratic(slope,fZero,fTwo)
00449 
00450                         else
00451 
00452                                 !lambda=cubic(slope,fZero,fOne,fTwo,lambda1,lambda2)
00453                                 lambda=quadratic(slope,fZero,fTwo)
00454                          
00455                         end if
00456 
00457 
00458 
00459                         !print*,lambda
00460                         !stop
00461 
00462 
00463                         if (lambda.lt.thetaINMin) then
00464                                 lambda=thetaINMin
00465                         elseif (lambda.gt.thetaINMax) then 
00466                                 lambda=thetaINMax
00467                         end if
00468 
00469 
00470                                 deltaX=lambda*deltaX
00471                                         
00472                                 eta=1.d0-lambda*(1.d0-eta)
00473 
00474                                 tempy=y
00475 
00476                                 tempy(setVec)=tempy(setVEc)+deltax
00477 
00478                                 F2=discrete(1)
00479 
00480                                 normNewFcurrent=dsqrt(dot_product(F2,F2))
00481                                 normNewFcurrent=wNorm(F2,wM)
00482 
00483                                 fTwo=0.5d0*normNewFcurrent**2.d0
00484                                 fTwo=0.5d0*dot_product(invscaler(F2,wM),invScaler(F2,wM))
00485 
00486 
00487                                 if (j.eq.maxBacktrackCount) then
00488 
00489                                         !if (gsSolver.eqv.no) print*,'not enough backtracks'
00490                                         !print*,'not enough backtracks'
00491                                         exit
00492 
00493                                 else
00494 
00495                                         cycle backtrack
00496 
00497                                 end if
00498 
00499                         end if
00500                         
00501 
00502 
00503 
00504 
00505                 end do backtrack
00506 
00507         end if
00508 
00509 
00510         !print*,normNewFcurrent,normRes,normFcurrent
00511         !stop
00512 
00513 
00514         wM=0.d0
00515 
00516 
00517         do j=1,globalcurrentsetlength
00518 
00519                 wM(j)=1.d0/(1e-3*dabs(y(setvec(j)))+1e-8)
00520 
00521         end do
00522 
00523         wM=scaler(deltax,wM)
00524 
00525 
00526         stepWNorm=dsqrt(dot_product(wM,wM))/globalcurrentsetlength
00527 
00528         if (stepWNorm.lt.1) then
00529 
00530                 stepLengthConverged=1
00531 
00532         else
00533 
00534                 stepLengthConverged=0
00535 
00536         end if
00537 
00538 
00539         if (gsSolver.eqv.no) print*,'stepWNorm=',stepWNorm
00540 
00541 
00542 
00543         select case(choiceIN)
00544 
00545                 case(1)
00546 
00547                         eta=dabs( normNewFcurrent  -  normRes)/normFcurrent
00548 
00549 
00550 
00551                         compareIN=etaPrev**((1.d0+dsqrt(5.d0))/2.d0)
00552 
00553                         if (compareIN.gt.0.1d0) then
00554                 
00555                                 eta=max(eta,compareIN)
00556 
00557                         end if
00558 
00559                         
00560                 case(2)
00561 
00562                         eta=gammaIN*(normNewFcurrent/normFcurrent)**alphaIN
00563 
00564                         compareIN=gammaIN*etaPrev**alphaIN
00565                 
00566                         if (compareIN.gt.0.1d0) then
00567                                         
00568                                 eta=max(eta,compareIN)
00569 
00570                         end if
00571 
00572 
00573                 end select
00574 
00575 
00576                 !if (gsSolver.eqv.no) print*,eta
00577                 !stop
00578 
00579         !print*,eta
00580         !pause
00581 
00582                 
00583         !eta=eta*1e1
00584 
00585 
00586         eta=min(eta,etaMax)
00587 
00588         
00589         etaPrev=eta
00590 
00591         if (gsSolver.eqv.no) print*,'current eta is=',eta
00592 
00593 
00594 
00595         if (gsSolver.eqv.no) print*,'normNewFcurrent=',normNewFcurrent,' <-------------'
00596         if (gsSolver.eqv.no) print*,' against conv tol=',newtonConvergenceTolerance
00597 
00598 
00599         write(2,'(F20.12,1X)')normNewFcurrent
00600 
00601 
00602         if ((normNewFcurrent.le.newtonConvergenceTolerance).and.(stepLengthConverged.eq.1)) then
00603 
00604                 if (gsSolver.eqv.no) print*,'newton converged'
00605                 if (gsSolver.eqv.no) print*,'LinearIterations=',totalLinearIterations
00606 
00607                 totalNonLinearIterations=newtonIter
00608 
00609                 newtonSuccess=1
00610                 exit
00611 
00612         elseif (normNewFcurrent.le.newtonConvergenceTolerance) then
00613 
00614                 if (gsSolver.eqv.no) print*,'step length NOT converged'
00615 
00616         elseif (stepLengthConverged.eq.1) then
00617 
00618                 if (gsSolver.eqv.no) print*,'F NOT converged'
00619 
00620 
00621         end if
00622 
00623 
00624         F0=F2
00625 
00626         !print*,F2
00627         !stop
00628 
00629         
00630         if (newtonIter.eq.maxNewtonIter) then
00631 
00632                 print*,'newton failed in',maxNewtonIter,' iterations'
00633                 print*,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
00634 
00635         end if
00636 
00637 
00638 
00639 
00640         
00641         ! if step is accepted
00642         
00643         !y=tempy
00644 
00645 
00646         if (resSet.eqv.yes) then
00647 
00648                 !call setRule
00649 
00650                 !call set
00651 
00652 
00653         end if
00654 
00655 
00656 end do
00657 
00658 close(2)
00659 
00660 end subroutine newton