demona: D:/demona/2.0/manual/doc/2-31set/processor/newton/newton2Set.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 real(8) setMin
00089 
00090 real(8) normNewFcurrentRMS
00091 
00092 integer mylength
00093 
00094 
00095 Open (2,File='newtonResidual.dat',Status='unknown')
00096 
00097 
00098 
00099 myJacInterval=jacInterval !< this way the original jacInterval is protected
00100 
00101 if (freezedJacobian.eqv.on) then
00102 
00103         
00104         myJacInterval=maxNewtonIter !< so the jacobian is calculated only at the first iteration 
00105 
00106         if (gsSolver.eqv.no) print*,'freezed jacobian'
00107 
00108 end if
00109 
00110 
00111 
00112 newtonTolerance=1e-5
00113 newtonTolerance=0.01d0
00114 !newtonTolerance=0.00001d0
00115 inexactLinearTolerance=1e-5
00116 
00117 
00118 eta=0.01d0
00119 etaZero=0.01d0
00120 
00121 tIN=0.0001d0
00122 
00123 thetaINMin=0.1d0
00124 thetaINMax=0.5d0
00125 
00126 alphaIN=1.5d0
00127 alphaIN=2.d0
00128 
00129 gammaIN=0.9d0
00130 gammaIN=1.d0
00131 
00132 etaMax=0.9d0
00133 
00134 choiceIN=1
00135 
00136 newtonSuccess=0
00137 
00138 !print*,setvec
00139 
00140 !dof=>globalcurrentsetlength
00141 
00142 
00143 
00144 allocate(deltax(globalcurrentsetlength))
00145 allocate(F0(globalcurrentsetlength))
00146 allocate(F2(globalcurrentsetlength))
00147 allocate(resVec(globalcurrentsetlength))
00148 allocate(wM(globalcurrentsetlength))
00149 
00150 deltax=0.d0
00151 
00152 
00153 
00154 
00155 F0=discrete(0)
00156 
00157 totalLinearIterations=0
00158 
00159 !print*,f0
00160 !stop
00161 !pause
00162 
00163 
00164 
00165 normNewFcurrent=dsqrt(dot_product(F0,F0))
00166 initialNorm=normNewFcurrent
00167 write(2,'(F20.12,1X)')normNewFcurrent
00168 
00169 newtonConvergenceTolerance=newtonTolerance
00170 newtonConvergenceTolerance=newtonTolerance*initialNorm
00171 
00172 inexactLinearTolerance=newtonConvergenceTolerance
00173 
00174 
00175 
00176 
00177 if (gsSolver.eqv.no) print*,'initial norm:',normNewFcurrent
00178 
00179 !stop
00180 if (normNewFcurrent.le.newtonConvergenceTolerance) then
00181 
00182         if (gsSolver.eqv.no) print*,'zero initial norm'
00183         if (gsSolver.eqv.no) print*,'newton return'
00184         return
00185 
00186 
00187 end if
00188 
00189 !print*,f0
00190 !stop
00191 
00192 !print*,globalcurrentsetlength
00193 !stop
00194 
00195 
00196 !< symbolic jacobian evaluation do detect the size of the crs-jacobian
00197 call symbolic(valsizeMax)
00198 ptrSize=globalcurrentsetlength+1
00199 
00200 eta=etaZero
00201 
00202 
00203 !print*,'etaprev!!!!!!!!!!!!!!!!!!!'
00204 !stop
00205 
00206 
00207 etaPrev=eta
00208 
00209 !print*,F0
00210 !print*,'----'
00211 normNewFcurrentRMS=normNewFcurrent/dsqrt(globalcurrentsetlength*1.d0)
00212 
00213 setMin=normNewFcurrentRMS
00214 
00215 mylength=globalcurrentsetlength
00216 
00217 do newtonIter=1,maxNewtonIter
00218 
00219         !if (newtonIter.eq.2) test=1
00220 
00221         !print*,newtoniter
00222 
00223 
00224         etaVec(newtonIter)=eta
00225 
00226         if (inexactNewton.eqv.On) inexactLinearTolerance= normNewFcurrent*eta
00227 
00228 
00229                 if (mod(newtonIter+myJacInterval-1,myJacInterval).eq.0) then
00230 
00231                         if (normFcurrent.le.JacSkipTol*newtonConvergenceTolerance) then
00232 
00233                                 if(jacSkip.eqv.on) then
00234                                         print*,'almost the same jacobian probably'
00235                                         print*,'its evaluation is skipped.'
00236                                         goto 1
00237                                 end if
00238 
00239                         end if
00240 
00241                         ! compute the jacobian
00242                         if (allocated(JacVal)) deallocate(JacVal)
00243                         if (allocated(JacColInd)) deallocate(JacColInd)
00244                         if (allocated(uptr)) deallocate(uptr)
00245                         if (allocated(JacRowPtr)) deallocate(JacRowPtr)
00246 
00247                         !discOrder=lowOrder
00248 
00249                         allocate(JacVal(valSizeMax),JacColInd(valSizeMax),JacRowPtr(ptrSize),uptr(globalcurrentsetlength))
00250                         
00251                         !symbolic=0
00252                         JacVal=0.d0
00253                         JacColInd=0
00254                         uptr=0
00255                         JacRowPtr=0
00256 
00257                         tempy=y
00258 
00259                         myModel%physics(currentPhysics)%evalOrder=>jacDiscOrder
00260 
00261                         call colorJacobianCRS(JacVal,JacColInd,JacRowPtr,globalcurrentsetlength,1,valSize,indSize,uptr,0,numberOfColors)
00262 
00263                         myModel%physics(currentPhysics)%evalOrder=>myModel%physics(currentPhysics)%discOrder
00264 
00265                         !open(9,file='row.dat',status='unknown')
00266                         !write(9,'(F20.12,1X)')JacVal(JacRowPtr(37):JacRowPtr(38)-1)
00267                         !close(9)
00268                         !if (newtoniter.eq.2) stop      
00269                         
00270 
00271                         !if (allocated(luval)) deallocate(luval)
00272                         !allocate(luVal(valSize))
00273                         !luval=0.d0
00274 
00275                         !luval=JacVal(1:valsize)
00276 
00277 
00278                         wM=0.d0
00279 
00280                         do jj=1,globalcurrentsetlength
00281 
00282                                 do ii=jacrowPtr(jj),jacrowPtr(jj+1)-1
00283                                         
00284                                         wM(jj)=dabs(jacVal(ii))+resvec(jj) ! absolute row sum of the jacobian 
00285 
00286 
00287                                 end do
00288 
00289 
00290                         end do
00291 
00292 
00293 
00294                         ! compute ilu(0)
00295                         !call basicilu0CRS(luval,JacColInd(1:valsize),JacRowPtr,globalcurrentsetlength,valSize,uptr)
00296                         call basicilu0CRS(jacval(1:valsize),JacColInd(1:valsize),JacRowPtr,globalcurrentsetlength,valSize,uptr)
00297 
00298                         !luval=JacVal(1:valsize)
00299 
00300                         !discOrder=myModel%physics(currentPhysics)%discOrder
00301 
00302                         JacVal(1:valsize)=JacVal(1:valsize)+epsilon(0.d0)
00303                         JacVal(1:valsize)=JacVal(1:valsize)-epsilon(0.d0)
00304 
00305                 end if
00306 
00307                 1 continue
00308 
00309                 !print*,'==============',valsize
00310                 !pause
00311                 
00312         tempy=y
00313 
00314 
00315 
00316 
00317         if (gsSolver.eqv.no) print*,'newton iter=',newtonIter
00318 
00319         !inexactLinearTolerance=1e-5
00320 
00321         !if (gsSolver.eqv.no) print*,'linear',inexactLinearTolerance
00322 
00323         
00324         !eta=1e-2
00325         
00326         !inexactLinearTolerance= normNewFcurrent*eta
00327 
00328 
00329         !print*,inexactLinearTolerance,normNewFcurrent,eta
00330         !pause
00331 
00332         deltax=0.d0
00333 
00334         !inexactLinearTolerance=eta
00335 
00336 
00337         !wM=jacval(uptr) !should be the diagonals
00338 
00339         wM=1.d0
00340 
00341         resVec=wM
00342 
00343         !inexactLinearTolerance=1e-5
00344 
00345         print*,'need local set, global set'
00346         !print*,maxval(setvec)
00347         stop
00348 
00349         !call bicgstabRP(-F0,deltax,inexactLinearTolerance,1000,resVec,luval,JacColInd,JacRowPtr,uptr,valsize,LinearIterations)
00350         call gmres(-F0(setvec),deltax,500,40,inexactLinearTolerance,resVec,JacVal(1:valsize),jaccolInd(1:valsize),jacrowPtr,uptr,valsize,LinearIterations,preconditioner)
00351 
00352 
00353         if (newtoniter.eq.2) then
00354         
00355         
00356         !do j=1,globalcurrentsetlength
00357 
00358         !       print*,setvec(j),deltax(j)
00359 
00360         !end do
00361         !stop
00362 
00363 
00364 
00365         end if
00366 
00367 
00368 
00369 
00370         !print*,deltax
00371         !pause
00372 
00373 
00374         !stop
00375         !wM=1.d0
00376         
00377         !print*,tempy(133),y(133),deltax(133)
00378 
00379         !totalLinearIterations=totalLinearIterations+LinearIterations
00380 
00381     lambda=1.d0
00382     lambda2=lambda
00383         
00384         !F0(setVec)=discrete(0)
00385 
00386         !if (newtonIter.eq.2) then
00387 
00388         !print*,F0
00389         !print*,'----'
00390         !stop
00391         !end if
00392 
00393         resVec=-F0(setVec)-matVec(deltax)
00394 
00395         normRes=dsqrt(dot_product(resVec,resVec))
00396         !normRes=wNorm(resVec,wM)
00397 
00398 
00399         tempy(setVec)=tempy(setVec)+lambda*deltax
00400         !y(setVec)=y(setVec)+deltax
00401 
00402         !print*,tempy(setVec(133)) !+lambda*deltax(133)
00403         !stop
00404 
00405 
00406         !here, we should move to the fullset to calculate F!!!
00407         myModel%physics(currentPhysics)%domain(currentDomain)%dofMatrix(:,:,:,flag)=1
00408 
00409         call currentSet
00410 
00411         print*,size(setvec)
00412 
00413 
00414         F2(setVec)=discrete(1)
00415 
00416 
00417         !print*,tempy(133),y(133),deltax(133)
00418         !stop
00419 
00420 
00421         !fZero=0.5d0*dot_product(F0,F0)
00422         !fzero=0.5d0*dot_product(invscaler(F0(setVec),wM),invScaler(F0(setVec),wM))
00423     
00424         !fTwo=0.5d0*dot_product(F2,F2)
00425         !fTwo=0.5d0*dot_product(invscaler(F2(setVec),wM),invScaler(F2(setVec),wM))
00426 
00427         !slope=-dot_product(F0(setVec),F0(setVec)) -dot_product(F0(setVec),resVec)
00428         !slope=-dot_product(invscaler(F0(setVec),wM),invScaler(F0(setVec),wM)) -dot_product(invscaler(F0(setVec),wM),invScaler(resVec,wM))
00429 
00430         normFcurrent=dsqrt(dot_product(F0,F0))
00431         !normFcurrent=wNorm(F0,wM)
00432 
00433         print*,'------'
00434         print*,normFcurrent
00435         !pause
00436 
00437         normNewFcurrent=dsqrt(dot_product(F2,F2))
00438         !normNewFcurrent=wNorm(F2,wM)
00439         print*,normNewFcurrent
00440         !pause
00441         
00442         totalLinearIterations=totalLinearIterations+LinearIterations
00443 
00444 
00445         if (doBacktrack.eqv.on) then
00446                 if (slope.ge.0.d0) then
00447 
00448                         !print*,'positive slope'
00449                         !print*,'line search fails'
00450                         
00451                         !print*,slope,dot_product(F0,F0),dot_product(F0,resVec)
00452                         !stop
00453                         !exit
00454 
00455                 end if
00456         end if
00457         
00458         !print*,normNewFcurrent,normFcurrent
00459         !stop
00460 
00461         
00462         if (doBacktrack.eqv.off) then
00463 
00464                 !y(setVec)=y(setVec)+deltax
00465                 y=tempy
00466                 !print*,normNewFcurrent,normFcurrent
00467                 !stop
00468 
00469         else
00470 
00471 backtrack:      do j=1,maxBacktrackCount
00472                                 
00473 
00474                                 
00475                         testLHS=normNewFcurrent
00476                         testRHS=(1.d0-tIn*(1.d0-eta))*normFcurrent
00477 
00478                         !print*,testLHS,testRHS
00479                         !stop
00480 
00481                         if (testLHS.le.testRHS) then
00482 
00483                                 normFprev=normFcurrent
00484                                 F0(setVec)=F0(setVec)+matVec(deltaX) ! use F0 as temporary vector so dot_product is simple
00485                                 
00486                                 normRes=dsqrt(dot_product(F0(setVec),F0(setVec)))
00487                                 !normRes=wNorm(F0,wM)
00488 
00489                                 y(setVec)=y(setVec)+deltax
00490 
00491                                 F2(setVec)=discrete(0)
00492 
00493                                 normNewFcurrent=dsqrt(dot_product(F2,F2))
00494                                 !normNewFcurrent=wNorm(F2,wM)
00495 
00496                                 !print*,normNewFcurrent,normFcurrent
00497                                 !stop
00498 
00499                                 exit
00500 
00501                         else
00502 
00503 
00504                         if (j.eq.1)  then
00505                                 
00506                                 lambda=quadratic(slope,fZero,fTwo)
00507 
00508                         else
00509 
00510                                 !lambda=cubic(slope,fZero,fOne,fTwo,lambda1,lambda2)
00511                                 lambda=quadratic(slope,fZero,fTwo)
00512                          
00513                         end if
00514 
00515 
00516 
00517                         !print*,lambda
00518                         !stop
00519 
00520 
00521                         if (lambda.lt.thetaINMin) then
00522                                 lambda=thetaINMin
00523                         elseif (lambda.gt.thetaINMax) then 
00524                                 lambda=thetaINMax
00525                         end if
00526 
00527 
00528                                 deltaX=lambda*deltaX
00529                                         
00530                                 eta=1.d0-lambda*(1.d0-eta)
00531 
00532                                 tempy=y
00533 
00534                                 tempy(setVec)=tempy(setVEc)+deltax
00535 
00536                                 F2(setVec)=discrete(1)
00537 
00538                                 normNewFcurrent=dsqrt(dot_product(F2,F2))
00539                                 normNewFcurrent=wNorm(F2(setVec),wM)
00540 
00541                                 fTwo=0.5d0*normNewFcurrent**2.d0
00542                                 fTwo=0.5d0*dot_product(invscaler(F2(setVec),wM),invScaler(F2(setVec),wM))
00543 
00544 
00545                                 if (j.eq.maxBacktrackCount) then
00546 
00547                                         !if (gsSolver.eqv.no) print*,'not enough backtracks'
00548                                         !print*,'not enough backtracks'
00549                                         exit
00550 
00551                                 else
00552 
00553                                         cycle backtrack
00554 
00555                                 end if
00556 
00557                         end if
00558                         
00559 
00560 
00561 
00562 
00563                 end do backtrack
00564 
00565         end if
00566 
00567 
00568         !print*,normNewFcurrent,normRes,normFcurrent
00569         !stop
00570 
00571 
00572         wM=0.d0
00573 
00574 
00575         !do j=1,globalcurrentsetlength
00576         do j=1,mylength
00577 
00578         
00579 
00580                 wM(j)=1.d0/(1e-3*dabs(y(setvec(j)))+1e-8)
00581 
00582         end do
00583 
00584         wM=scaler(deltax,wM)
00585 
00586 
00587         !stepWNorm=dsqrt(dot_product(wM,wM))/globalcurrentsetlength
00588         stepWNorm=dsqrt(dot_product(wM,wM))/mylength
00589 
00590         
00591 
00592         if (stepWNorm.lt.1) then
00593 
00594                 stepLengthConverged=1
00595 
00596         else
00597 
00598                 stepLengthConverged=0
00599 
00600         end if
00601 
00602 
00603         if (gsSolver.eqv.no) print*,'stepWNorm=',stepWNorm
00604 
00605 
00606 
00607 
00608         select case(choiceIN)
00609 
00610                 case(1)
00611 
00612                         eta=dabs( normNewFcurrent  -  normRes)/normFcurrent
00613 
00614 
00615 
00616                         compareIN=etaPrev**((1.d0+dsqrt(5.d0))/2.d0)
00617 
00618                         if (compareIN.gt.0.1d0) then
00619                 
00620                                 eta=max(eta,compareIN)
00621 
00622                         end if
00623 
00624                         
00625                 case(2)
00626 
00627                         eta=gammaIN*(normNewFcurrent/normFcurrent)**alphaIN
00628 
00629                         compareIN=gammaIN*etaPrev**alphaIN
00630                 
00631                         if (compareIN.gt.0.1d0) then
00632                                         
00633                                 eta=max(eta,compareIN)
00634 
00635                         end if
00636 
00637 
00638                 end select
00639 
00640 
00641                 !if (gsSolver.eqv.no) print*,eta
00642                 !stop
00643 
00644         !print*,eta
00645         !pause
00646 
00647                 
00648         !eta=eta*1e1
00649 
00650 
00651         eta=min(eta,etaMax)
00652 
00653         
00654         etaPrev=eta
00655 
00656         if (gsSolver.eqv.no) print*,'current eta is=',eta
00657 
00658 
00659 
00660         if (gsSolver.eqv.no) print*,'normNewFcurrent=',normNewFcurrent,' <-------------'
00661         if (gsSolver.eqv.no) print*,' against conv tol=',newtonConvergenceTolerance
00662 
00663 
00664         write(2,'(F20.12,1X)')normNewFcurrent
00665 
00666 
00667         if ((normNewFcurrent.le.newtonConvergenceTolerance).and.(stepLengthConverged.eq.1)) then
00668 
00669                 if (gsSolver.eqv.no) print*,'newton converged'
00670                 if (gsSolver.eqv.no) print*,'LinearIterations=',totalLinearIterations
00671 
00672                 totalNonLinearIterations=newtonIter
00673 
00674                 newtonSuccess=1
00675                 exit
00676 
00677         elseif (normNewFcurrent.le.newtonConvergenceTolerance) then
00678 
00679                 if (gsSolver.eqv.no) print*,'step length NOT converged'
00680 
00681         elseif (stepLengthConverged.eq.1) then
00682 
00683                 if (gsSolver.eqv.no) print*,'F NOT converged'
00684 
00685 
00686         end if
00687 
00688 
00689         F0=F2
00690         
00691         !if (newtonIter.eq.1) then
00692 
00693         !print*,F0
00694         !print*,'----'
00695 
00696         !end if
00697         !stop
00698 
00699         
00700         if (newtonIter.eq.maxNewtonIter) then
00701 
00702                 print*,'newton failed in',maxNewtonIter,' iterations'
00703                 print*,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
00704 
00705         end if
00706 
00707 
00708 
00709 
00710         
00711         ! if step is accepted
00712         
00713         !y=tempy
00714 
00715 
00716         if (resSet.eqv.yes) then
00717 
00718                 !print*,f2
00719                 !print*,dsqrt(dot_product(F2,F2)),normNewFcurrent
00720 
00721 
00722                 open(9,file='F2.dat',status='unknown')
00723                 write(9,'(E20.12,1X)')F2
00724                 close(9)
00725 
00726 
00727                 normNewFcurrentRMS=normNewFcurrent/dsqrt(globalcurrentsetlength*1.d0)
00728 
00729                 !print*,setMin
00730 
00731                 setMin=min(setMin,normNewFcurrentRMS)
00732 
00733                 !print*,setMin,normNewFcurrentRMS
00734                 !print*,normNewFcurrentRMS
00735 
00736 
00737                 do j=1,1
00738 
00739                 call setRule(F2,size(F2),eta,normNewFcurrentRMS)
00740                 !call setRule(F2,size(F2),newtonTolerance,normNewFcurrentRMS)
00741                 !call setRule(F2,size(F2),normNewFcurrent,setMin)
00742 
00743                 !call set
00744 
00745 
00746                 print*,globalcurrentsetlength
00747 
00748                 call currentSet
00749 
00750                 print*,globalcurrentsetlength
00751 
00752                 mylength=globalcurrentsetlength
00753 
00754                 F2(setVec)=discrete(1)
00755                 
00756                 print*,normNewFcurrentRMS
00757                 print*,dsqrt(dot_product(F2,F2))/dsqrt(size(F2)*1.d0)
00758                 print*,dsqrt(dot_product(F2(setVec),F2(setVec)))/dsqrt(globalcurrentsetlength*1.d0)
00759 
00760                 
00761                 end do
00762 
00763                 !pause
00764 
00765 
00766         if (allocated(deltax)) deallocate(deltax)
00767         if (allocated(resVec)) deallocate(resVec)
00768         if (allocated(wM)) deallocate(wM)
00769                 
00770         allocate(deltax(globalcurrentsetlength))
00771         allocate(resVec(globalcurrentsetlength))
00772         allocate(wM(globalcurrentsetlength))
00773 
00774 
00775         end if
00776 
00777 
00778 
00779 
00780 
00781         !if (newtoniter.eq.3) stop
00782 
00783 
00784 end do
00785 
00786 close(2)
00787 
00788 end subroutine newton