Compute Equilibrium, iteratively
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(t_functional_mhd3d), | intent(inout) | :: | sf |
SUBROUTINE MinimizeMHD3D_descent(sf) ! MODULES USE MODgvec_MHD3D_Vars USE MODgvec_MHD3D_EvalFunc USE MODgvec_Analyze, ONLY:analyze USE MODgvec_Restart, ONLY:WriteState USE MODgvec_MHD3D_visu, ONLY:WriteSFLoutfile IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- ! OUTPUT VARIABLES CLASS(t_functional_mhd3d), INTENT(INOUT) :: sf !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES INTEGER :: iter,nStepDecreased,nSkip_Jac,nSkip_dw INTEGER :: JacCheck,lastoutputIter,StartTimeArray(8) REAL(wp) :: dt,deltaW,absTol INTEGER,PARAMETER :: ndamp=10 REAL(wp) :: tau(1:ndamp), tau_bar REAL(wp) :: min_dt_out,max_dt_out,min_dw_out,max_dw_out,sum_dW_out,t_pseudo,Fnorm(3),Vnorm(3),Fnorm0(3),Fnorm_old(3),W_MHD3D_0 INTEGER :: logUnit !globally needed for logging INTEGER :: logiter_ramp,logscreen LOGICAL :: restart_iter LOGICAL :: first_iter !=================================================================================================================================== SWRITE(UNIT_stdOut,'(A)') "MINIMIZE MHD3D FUNCTIONAL..." abstol=minimize_tol dt=start_dt nstepDecreased=0 nSkip_Jac=0 t_pseudo=0 lastOutputIter=0 iter=0 Vnorm=0.0_wp logiter_ramp=1 logscreen=1 first_iter=.TRUE. restart_iter=.FALSE. CALL U(-3)%set_to(U(0)) !initial state, should remain unchanged DO WHILE(iter.LT.maxIter) IF((first_iter).OR.(restart_iter))THEN JacCheck=1 !abort if detJ<0 CALL EvalAux( U(0),JacCheck) U(0)%W_MHD3D=EvalEnergy(U(0),.FALSE.,JacCheck) W_MHD3D_0 = U(0)%W_MHD3D CALL EvalForce( U(0),.FALSE.,JacCheck,F(0)) Fnorm0=SQRT(F(0)%norm_2()) Fnorm=Fnorm0 Fnorm_old=1.1_wp*Fnorm0 CALL U(-1)%set_to(U(0)) !last state CALL U(-2)%set_to(U(0)) !state at last logging interval !for hirshman method IF(MinimizerType.EQ.10)THEN CALL V(-1)%set_to(0.0_wp) CALL V( 0)%set_to(0.0_wp) tau(1:ndamp)=0.15_wp/dt tau_bar = 0.075_wp END IF min_dt_out=1.0e+30_wp max_dt_out=0.0_wp min_dW_out=1.0e+30_wp max_dW_out=-1.0e+30_wp sum_dW_out=0.0_wp nSkip_dW =0 IF(restart_iter) restart_iter=.FALSE. IF(first_iter)THEN CALL StartLogging() first_iter=.FALSE. END IF END IF !before first iteration or after restart Jac<0 !COMPUTE NEW SOLUTION P(1) as a prediction SELECT CASE(MinimizerType) CASE(0) !gradient descent, previously used for minimizerType=0 CALL P(1)%AXBY(1.0_wp,U(0),dt,F(0)) !overwrites P(1), predicts solution U(1) CASE(10) !hirshman method !tau is damping parameter tau(1:ndamp-1) = tau(2:ndamp) !save old tau(ndamp) = MIN(0.15_wp,ABS(LOG(SUM(Fnorm**2)/SUM(Fnorm_old**2))))/dt !ln(|F_n|^2/|F_{n-1}|^2), Fnorm=|F_X1|,|F_X2|,|F_LA| tau_bar = 0.5_wp*dt*SUM(tau)/REAL(ndamp,wp) !=1/2 * tauavg CALL V(1)%AXBY(((1.0_wp-tau_bar)/(1.0_wp+tau_bar)),V(0),(dt/(1.0_wp+tau_bar)),F(0)) !velocity V(1) CALL P(1)%AXBY(1.0_wp,U(0),dt,V(1)) !overwrites P(1), predicst solution U(1) Vnorm=SQRT(V(1)%norm_2()) END SELECT JacCheck=2 !no abort,if detJ<0, JacCheck=-1 P(1)%W_MHD3D=EvalEnergy(P(1),.TRUE.,JacCheck) IF(JacCheck.EQ.-1)THEN dt=0.9_wp*dt nstepDecreased=nStepDecreased+1 nSkip_Jac=nSkip_Jac+1 restart_iter=.TRUE. CALL U(0)%set_to(U(-3)) !reset to initial state SWRITE(UNIT_stdOut,'(8X,I8,A,E11.4,A)')iter,'...detJac<0, decrease stepsize to dt=',dt, ' and RESTART simulation!!!!!!!' ELSE !detJ>0 deltaW=P(1)%W_MHD3D-U(0)%W_MHD3D!should be <=0, IF(deltaW.LE.dW_allowed*W_MHD3D_0)THEN !valid step /hirshman method accept W increase! IF(ALL(Fnorm.LE.abstol))THEN CALL Logging(.FALSE.) SWRITE(UNIT_stdOut,'(4x,A)')'==>Iteration finished, |force| in relative tolerance' EXIT !DO LOOP END IF iter=iter+1 t_pseudo=t_pseudo+dt ! for simple gradient & hirshman CALL U(-1)%set_to(U(0)) CALL U(0)%set_to(P(1)) ! for hirshman method IF(MinimizerType.EQ.10)THEN CALL V(-1)%set_to(V(0)) CALL V(0)%set_to(V(1)) END IF CALL EvalForce(P(1),.FALSE.,JacCheck,F(0)) !evalAux was already called on P(1)=U(0), so that its set false here. Fnorm_old=Fnorm Fnorm=SQRT(F(0)%norm_2()) nstepDecreased=0 min_dt_out=MIN(min_dt_out,dt) max_dt_out=MAX(max_dt_out,dt) min_dW_out=MIN(min_dW_out,deltaW) max_dW_out=MAX(max_dW_out,deltaW) sum_dW_out=sum_dW_out+deltaW IF(MOD(iter,logIter_ramp).EQ.0)THEN CALL Logging(.NOT.((logIter_ramp.GE.logIter).AND.(MOD(logscreen,nLogScreen).EQ.0))) IF(.NOT.(logIter_ramp.LT.logIter))THEN !only reset for logIter logscreen=logscreen+1 min_dt_out=1.0e+30_wp max_dt_out=0.0_wp min_dW_out=1.0e+30_wp max_dW_out=-1.0e+30_wp sum_dW_out=0.0_wp nSkip_dW =0 END IF logIter_ramp=MIN(logIter,logIter_ramp*2) END IF ELSE !not a valid step, decrease timestep and skip P(1) dt=0.9_wp*dt nstepDecreased=nStepDecreased+1 nSkip_dW=nSkip_dW+1 !CALL U(0)%set_to(U(-2)) restart_iter=.TRUE. SWRITE(UNIT_stdOut,'(8X,I8,A,E8.1,A,E11.4)')iter,'...deltaW>',dW_allowed,'*W_MHD3D_0, skip step and decrease stepsize to dt=',dt END IF END IF !JacCheck IF(nStepDecreased.GT.130) THEN ! 0.9^130 ~10^-6 SWRITE(UNIT_stdOut,'(A,E21.11)')'Iteration stopped since timestep has been decreased by 0.9^130: ', dt SWRITE(UNIT_stdOut,fmt_sep) RETURN END IF IF((MOD(iter,outputIter).EQ.0).AND.(lastoutputIter.NE.iter))THEN __PERFON('output') SWRITE(UNIT_stdOut,'(A)')'########################## OUTPUT ##################################' CALL Analyze(iter) CALL WriteState(U(0),iter) SWRITE(UNIT_stdOut,'(A)')'#####################################################################' lastOutputIter=iter __PERFOFF('output') END IF END DO !iter IF(iter.GE.MaxIter)THEN SWRITE(UNIT_stdOut,'(A,E21.11)')"maximum iteration count exceeded" END IF SWRITE(UNIT_stdOut,'(A)') "... DONE." SWRITE(UNIT_stdOut,fmt_sep) IF(lastoutputIter.NE.iter)THEN CALL Analyze(MIN(iter,MaxIter)) CALL WriteState(U(0),MIN(iter,MaxIter)) END IF CALL FinishLogging() CALL writeSFLoutfile(U(0),MIN(iter,MaxIter)) CONTAINS !================================================================================================================================= !> all screen and logfile tasks, can use all variables from subroutine above !! !================================================================================================================================= SUBROUTINE StartLogging() USE MODgvec_Globals, ONLY: GETFREEUNIT USE MODgvec_Output_Vars, ONLY: ProjectName,outputLevel USE MODgvec_MHD3D_visu, ONLY: checkAxis IMPLICIT NONE !--------------------------------------------------------------------------------------------------------------------------------- CHARACTER(LEN=255) :: fileString INTEGER :: TimeArray(8),iLogDat REAL(wp) :: AxisPos(2,2),W_MHD3D INTEGER,PARAMETER :: nLogDat=20 REAL(wp) :: LogDat(1:nLogDat) !================================================================================================================================= IF(.NOT.MPIroot) RETURN __PERFON('log_output') W_MHD3D=U(0)%W_MHD3D CALL DATE_AND_TIME(values=TimeArray) ! get System time WRITE(UNIT_stdOut,'(A,E11.4,A)')'%%%%%%%%%% START ITERATION, dt= ',dt, ' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%' WRITE(UNIT_stdOut,'(A,I4.2,"-",I2.2,"-",I2.2,1X,I2.2,":",I2.2,":",I2.2)') & '%%% Sys date : ',timeArray(1:3),timeArray(5:7) WRITE(UNIT_stdOut,'(A,3E21.14)') & '%%% dU = |Force|= ',Fnorm(1:3) IF(MinimizerType.EQ.10) THEN WRITE(UNIT_stdOut,'(A,E11.4,A,3E11.4)') & '%%% accel.GD: tau= ',tau_bar,' |vel|= ',Vnorm(1:3) END IF WRITE(UNIT_stdOut,'(40(" -"))') !------------------------------------ StartTimeArray=TimeArray !save first time stamp logUnit=GETFREEUNIT() WRITE(FileString,'("logMinimizer_",A,"_",I4.4,".csv")')TRIM(ProjectName),outputLevel OPEN(UNIT = logUnit ,& FILE = TRIM(FileString) ,& STATUS = 'REPLACE' ,& ACCESS = 'SEQUENTIAL' ) !header iLogDat=0 WRITE(logUnit,'(A)',ADVANCE="NO")'"#iterations","runtime(s)","min_dt","max_dt"' WRITE(logUnit,'(A)',ADVANCE="NO")',"W_MHD3D","min_dW","max_dW","sum_dW"' WRITE(logUnit,'(A)',ADVANCE="NO")',"normF_X1","normF_X2","normF_LA"' LogDat(ilogDat+1:iLogDat+11)=(/0.0_wp,0.0_wp,dt,dt,W_MHD3D,0.0_wp,0.0_wp,0.0_wp,Fnorm(1:3)/) iLogDat=11 IF(MinimizerType.EQ.10) THEN WRITE(logUnit,'(A)',ADVANCE="NO")',"tau","normV_X1","normV_X2","normV_LA"' LogDat(ilogDat+1:iLogDat+4)=(/tau_bar,Vnorm(1:3)/) iLogDat=iLogDat+4 END IF IF(doCheckDistance) THEN WRITE(logUnit,'(A)',ADVANCE="NO")',"max_Dist","avg_Dist"' LogDat(iLogDat+1:iLogDat+2)=(/0.0_wp,0.0_wp/) iLogDat=iLogDat+2 END IF!doCheckDistance IF(doCheckAxis) THEN WRITE(logUnit,'(A)',ADVANCE="NO")',"X1_axis_0","X2_axis_0","X1_axis_1","X2_axis_1"' CALL CheckAxis(U(0),2,AxisPos) LogDat(iLogDat+1:iLogDat+4)=RESHAPE(AxisPos,(/4/)) iLogDat=iLogDat+4 END IF!doCheckAxis WRITE(logUnit,'(A)')' ' !first data line WRITE(logUnit,'(*(e23.15,:,","))') logDat(1:iLogDat) __PERFOFF('log_output') END SUBROUTINE StartLogging !================================================================================================================================= !> all screen and logfile tasks, can use all variables from subroutine above !! !================================================================================================================================= SUBROUTINE Logging(quiet) USE MODgvec_MHD3D_visu, ONLY: checkDistance USE MODgvec_MHD3D_visu, ONLY: checkAxis IMPLICIT NONE LOGICAL, INTENT(IN) :: quiet !! True: no screen output !--------------------------------------------------------------------------------------------------------------------------------- INTEGER :: TimeArray(8),runtime_ms,iLogDat REAL(wp) :: AxisPos(2,2),maxDist,avgDist,W_MHD3D INTEGER,PARAMETER :: nLogDat=20 REAL(wp) :: LogDat(1:nLogDat) !================================================================================================================================= IF(.NOT.MPIroot) RETURN __PERFON('log_output') CALL DATE_AND_TIME(values=TimeArray) ! get System time W_MHD3D=U(0)%W_MHD3D IF(.NOT.quiet)THEN WRITE(UNIT_stdOut,'(80("%"))') WRITE(UNIT_stdOut,'(A,I4.2,"-",I2.2,"-",I2.2,1X,I2.2,":",I2.2,":",I2.2)') & '%%% Sys date : ',timeArray(1:3),timeArray(5:7) WRITE(UNIT_stdOut,'(A,I8,A,2I8,A,E11.4,A,2E11.4,A,E21.14,A,3E12.4)') & '%%% #ITERATIONS= ',iter,', #skippedIter (Jac/dW)= ',nSkip_Jac,nSkip_dW, & '\n%%% t_pseudo= ',t_pseudo,', min/max dt= ',min_dt_out,max_dt_out, & '\n%%% W_MHD3D= ',W_MHD3D,', min/max/sum deltaW= ' , min_dW_out,max_dW_out,sum_dW_out WRITE(UNIT_stdOut,'(A,3E21.14)') & '%%% dU = |Force|= ',Fnorm(1:3) !------------------------------------ END IF!.NOT.quiet iLogDat=0 runtime_ms=MAX(0,SUM((timeArray(5:8)-StartTimearray(5:8))*(/360000,6000,100,1/))) LogDat(ilogDat+1:iLogDat+11)=(/REAL(iter,wp),REAL(runtime_ms,wp)/100.0_wp, & min_dt_out,max_dt_out,W_MHD3D,min_dW_out,max_dW_out,sum_dW_out, & Fnorm(1:3)/) iLogDat=11 IF(MinimizerType.EQ.10) THEN IF(.NOT.quiet)THEN WRITE(UNIT_stdOut,'(A,E11.4,A,3E11.4)') & '%%% accel.GD: tau= ',tau_bar,' |vel|= ',Vnorm(1:3) END IF!.NOT.quiet LogDat(ilogDat+1:iLogDat+4)=(/tau_bar,Vnorm(1:3)/) iLogDat=iLogDat+4 END IF IF(doCheckDistance) THEN CALL CheckDistance(U(0),U(-2),maxDist,avgDist) CALL U(-2)%set_to(U(0)) IF(.NOT.quiet)THEN WRITE(UNIT_stdOut,'(A,2E11.4)') & ' %%% Dist to last log (max/avg) : ',maxDist,avgDist END IF!.NOT.quiet LogDat(iLogDat+1:iLogDat+2)=(/maxDist,avgDist/) iLogDat=iLogDat+2 END IF!doCheckDistance IF(doCheckAxis) THEN CALL CheckAxis(U(0),2,AxisPos) IF(.NOT.quiet)THEN WRITE(UNIT_stdOut,'(2(A,2E22.14))') & '%%% axis position (X1,X2,zeta=0 ): ',AxisPos(1:2,1), & '\n%%% axis position (X1,X2,zeta=pi/nfp): ',AxisPos(1:2,2) END IF!.NOT.quiet LogDat(iLogDat+1:iLogDat+4)=RESHAPE(AxisPos,(/4/)) iLogDat=iLogDat+4 END IF !doCheckAxis IF(.NOT.quiet)THEN WRITE(UNIT_stdOut,'(40(" -"))') END IF!.NOT.quiet WRITE(logUnit,'(*(e23.15,:,","))') logDat(1:iLogDat) __PERFOFF('log_output') END SUBROUTINE Logging !================================================================================================================================= !> !! !================================================================================================================================= SUBROUTINE FinishLogging() IMPLICIT NONE !--------------------------------------------------------------------------------------------------------------------------------- CLOSE(logUnit) END SUBROUTINE FinishLogging END SUBROUTINE MinimizeMHD3D_descent