Two Adaptive Timestepping questions/suggestions

[Franz Pichler]

Hello Everbody,
Hello dear Elmer Team,
i hope everybody had a good summer and that you guys missed me as much as i missed you. i had a bit of an Elmer braeak the last months but will start hacking around slowly now.

I just read through the ElmerSolver file and the adaptive timestepping therein. I am a big fan of thew "two small steps, one big , compare" scheme, that i used myself already.

Ihave two questions to the Elmer implementation now:
1)
When the step is accepted wouldn't it be better to take the solution from the second small step? because this comes from a finer time discretisation?
I think that would give some cheap performance improvement in transient simulations.

2)
The second thing is the way the error is calculated. it is calculated by the difference of the norms. Wouldn't the norm of the differences be a better measure? But maybe it's hard to calculate the norm of the difference in the Elmersolver module?

As always please don't blame me if i got everything wrong.

i am looking forward to your reply and wish a nice weekend.
best regards
Franz

[Franz Pichler]

OOops, SOrry.
Forget about my first point. Elmer does it the good way, of course! Sorry for being a wise guy.
Still the second would be interesting.

And as another suggestion
3)
what about introducing another limit. so that if the error is between the two limits the step is accepted but the timestep is not changed.
if the error is smaller than both limits the step is accepted and the timestep is increased.
and if the error is bigger than both limits than step is not accepted and timestep is decreased.
That way if the limits are chosen carefully often the simulation "finds" a good timestep by itself.

so
sorry for the mistake again
and have a nice one
Franz

[mzenker]

Hi,

the problem with the adaptive timestepping as it is now in Elmer is that it often makes the simulation time explode. As in some kinds of simulations, small timesteps are needed in order to achieve convergence at all, a suggestion from my side is to take convergence as additional criterion: Decrease timestep if solution did not converge.
I have implemented that one in a custom solver and use it frequently.

Matthias

[nguyent]

@mzenker: Could you share your code? The current adaptive timestepping continues to bisect the time until it reaches the limit in my cases (MEMS). I was just getting ready to try to write an adaptive timestepper dependent on convergence and seeing your code would really be great.

[mzenker]

Hi,

my solver does actually some other things as well, but I can share the relevant code pieces.
The idea is to check convergence after each timestep. If the calculation did not converge, the time and result is reset to previous, and the timestep is divided by two. If Resultoutput solver is used, a variable is set in order not to save the nonconverged result. As I don't use ElmerPost any more, I didn't take care of .ep output.
All this would ideally be done in main Elmersolver, but as I was not able to compile it at the time, I implemented the functionality in a custom solver.
I have made an extract of my solver which should only do the adaptive timestepping. I did not test it, so there may be some bugs and/or syntax errors in it. Here it is:

Code: Select all

SUBROUTINE AdaptiveSolver(Model, Solver, dt, Transient)

	USE DefUtils
	IMPLICIT NONE
	TYPE(Solver_t) :: Solver
	TYPE(Model_t) :: Model
	REAL(KIND=dp) :: dt
	LOGICAL :: Transient
	!----------------------------------------------------------------
	! Local variables
	!----------------------------------------------------------------
	TYPE(ValueList_t), POINTER :: SolverParams
    TYPE(Mesh_t), POINTER :: Mesh
    TYPE(Variable_t), POINTER :: Var
	REAL(KIND=dp), POINTER :: WrkPntr(:)
	INTEGER :: SaveResults = 1
	REAL(KIND=dp) :: MaxTimeStepSize = 0.0_dp, MinTimeStepSize, CurrentTimeStepSize=0.0_dp, MinSaveInterval
	REAL(KIND=dp) :: Time, Last_Time = 0.0_dp, Last_Save_Time = 0.0_dp, Duration
	REAL(KIND=dp) :: SaveResults_real
	LOGICAL :: Found, Visited = .FALSE.
	LOGICAL :: AdaptiveTimeStepping, Converged, StopSimulation = .FALSE.
	INTEGER :: istat
    TYPE(Solver_t), POINTER :: iSolver
    REAL(KIND=dp), ALLOCATABLE :: xx(:,:), prevxx(:,:,:)
	INTEGER :: i, j, k, n
	INTEGER :: AdaptiveKeepSmallest, KeepTimeStepSize = 0
	INTEGER :: SteadyConverged
	
	SAVE Last_Time
	SAVE CurrentTimeStepSize
	SAVE KeepTimeStepSize
	SAVE Last_Save_Time
	SAVE Visited
	SAVE xx, prevxx
	SAVE StopSimulation

	CALL Info( 'AdaptiveSolver', '*************************************************' )

	!Read in solver parameters
	!-------------------------
	SolverParams => GetSolverParams()
	IF (.NOT. ASSOCIATED(SolverParams)) CALL FATAL('AdaptiveSolver','No Solver section found')
	
	IF ( .NOT. Visited ) THEN
		! allocate space for previous result
		n = Model % NumberOfSolvers
		j = 0
		k = 0
		DO i=1,n
			iSolver => Model % Solvers(i)
			IF ( ASSOCIATED( iSolver % Variable  % Values ) ) THEN
				IF ( ASSOCIATED( iSolver % Variable % PrevValues ) ) THEN
					j = MAX( j, SIZE( iSolver % Variable % PrevValues,2 ) )
				END IF
				k = MAX( k, SIZE( iSolver % Variable % Values ) )
			END IF
		END DO
		ALLOCATE( xx(n,k), prevxx( n,k,j ) )
		! initialize variables, allocate memory
		Mesh => Model % Meshes
		DO WHILE( ASSOCIATED(Mesh) )
			IF ( Mesh % OutputActive ) EXIT 
			Mesh => Mesh % Next
		END DO
		CALL SetCurrentMesh( Model, Mesh )
		NULLIFY(WrkPntr)
		ALLOCATE(WrkPntr(1),STAT=istat)
		IF( istat /= 0 ) CALL Fatal('AdaptiveSolver','Memory allocation error') 	
		CALL VariableAdd( Model % Variables, Mesh, Solver, 'AdaptiveTimeStepSize', 1, WrkPntr )
		NULLIFY(WrkPntr)
		ALLOCATE(WrkPntr(1),STAT=istat)
		IF( istat /= 0 ) CALL Fatal('AdaptiveSolver','Memory allocation error 3') 	
		CALL VariableAdd( Model % Variables, Mesh, Solver, 'SaveResults', 1, WrkPntr )
	END IF
	
	Time = GetTime()
	
	IF ( StopSimulation ) THEN

		! Set SaveResults variable to zero
		Var => VariableGet( Model % Variables, 'SaveResults' )
		IF(.NOT. ASSOCIATED(Var)) THEN
			CALL FATAL( 'AdaptiveSolver', 'No memory allocated for SaveResults')       
		END IF
		! It seems impossible to save an Integer, so we convert it to Real
		Var % Values(1) = 0.0_dp
	
		! End of simulation was reached at previous timestep - now send stop signal
		WRITE(Message, *)  'End of Simulation! Sending stop signal to solver...'
		CALL ListAddConstReal(Model % Simulation,'Exit Condition',1.0_dp)
		
	ELSE
	
		Converged = .TRUE.
		SaveResults = 1
		DO i=1,Model % NumberOfSolvers
			iSolver => Model % Solvers(i)
			SteadyConverged = iSolver % Variable % SteadyConverged
			IF( ( Time > 0 ) .AND. ( SteadyConverged == 0 ) ) THEN
				Converged = .FALSE.
				CALL Info('AdaptiveSolver', 'No convergence for this timestep')
				EXIT
			END IF
		END DO
		
		AdaptiveTimeStepping = ListGetLogical(SolverParams,'Convergence Adaptive Timestepping',Found)
		
		IF ( AdaptiveTimeStepping ) THEN
			MinTimestepSize = ListGetConstReal( SolverParams, 'Convergence Adaptive Min Timestep', Found )
			IF ( .NOT. Found ) MinTimeStepSize = 1e-15
			AdaptiveKeepSmallest = GetInteger(SolverParams, 'Convergence Adaptive Keep Smallest', Found )
			IF ( .NOT. Found ) AdaptiveKeepSmallest = 2

			IF( Converged .OR. CurrentTimeStepSize - MinTimeStepSize < 1e-10 ) THEN
				! save current result
				DO i=1,Model % NumberOFSolvers
					iSolver => CurrentModel % Solvers(i)
					IF ( ASSOCIATED( iSolver % Variable % Values ) ) THEN
						n = SIZE( iSolver % Variable % Values )
						xx(i,1:n) = iSolver % Variable % Values
						IF ( ASSOCIATED( iSolver % Variable % PrevValues ) ) THEN
							DO j=1,SIZE( iSolver % Variable % PrevValues,2 )
								prevxx(i,1:n,j) = iSolver % Variable % PrevValues(:,j)
							END DO
						END IF
					END IF
				END DO
			ELSE
				! - reset to result for last convergence
				DO i=1,Model % NumberOFSolvers
					iSolver => CurrentModel % Solvers(i)
					IF ( ASSOCIATED( iSolver % Variable % Values ) ) THEN
						n = SIZE(iSolver % Variable % Values)
						iSolver % Variable % Values = xx(i,1:n)
						IF ( ASSOCIATED( iSolver % Variable % PrevValues ) ) THEN
							DO j=1,SIZE( iSolver % Variable % PrevValues,2 )
							iSolver % Variable % PrevValues(:,j) = prevxx(i,1:n,j)
							END DO
						END IF
					END IF
				END DO
				
				! - reset time to previous value
				Var  => VariableGet( Model % Solver % Mesh % Variables, 'Time' )
				IF ( ASSOCIATED( Var ) )  Var % Values(1)  = Time - CurrentTimeStepSize
				Time = GetTime()
				Write(Message, *) 'Adaptive Timestepping: Time reset to ', Time
				CALL Info('AdaptiveSolver', Message)
				SaveResults = 0
			END IF
		END IF ! AdaptiveTimeStepping

		! Get duration 
		Duration = ListGetConstReal( SolverParams, 'Duration',Found)
		IF ( .NOT. Found ) Duration = 0.0

		! Get time minimum saving interval 
		MinSaveInterval = ListGetConstReal(SolverParams, 'Min Data Saving Interval',Found)
		IF ( .NOT. Found ) MinSaveInterval = 0.0_dp
		
		IF( SaveResults /= 0 ) THEN
			IF( ( Time > 1e-10 ) .AND. ( Time - Last_Save_Time < MinSaveInterval ) ) THEN
				SaveResults = 0
			ELSE
				Last_Save_Time = Time
			END IF
		END IF
				
		! Check if duration is exceeded since Time was saved the last time
		! Duration is multiplied by a factor to avoid problems due to rounding errors
	    IF ( Time - Last_Time >= Duration - 1e-10 ) THEN
		! End of simulation after this timestep! Save one last time...
			SaveResults = 1
			StopSimulation = .TRUE.
		END IF
				
		! Get time step size
		MaxTimeStepSize = ListGetConstReal(SolverParams, 'Adaptive TimeStep Size',Found)
		IF ( .NOT. Found ) MaxTimeStepSize = 0.0
			
		! Determine timestep size
		IF ( .NOT. Visited ) THEN	
			CurrentTimeStepSize = MaxTimeStepSize
		ELSE IF ( AdaptiveTimeStepping ) THEN
			IF ( Converged ) THEN
				! check keepsmallest, increase timestep size if desired
				IF ( KeepTimeStepSize > 0 ) THEN
					KeepTimeStepSize = KeepTimeStepSize - 1
				ELSE
					! try to double the timestep until the timestep set by the user is reached
					CurrentTimeStepSize = CurrentTimeStepSize * 2
					IF ( CurrentTimeStepSize > MaxTimeStepSize ) THEN
						CurrentTimeStepSize = MaxTimeStepSize
					ELSE
						WRITE( Message, * ), 'Adaptive Timestepping: Time step set to ', CurrentTimeStepSize
						CALL Info('AdaptiveSolver', Message)
					END IF
				END IF
			ELSE
				! set timestep size to half if MinTimestepSize is not reached
				IF( CurrentTimeStepSize - MinTimeStepSize > 1e-10 ) THEN
					CurrentTimeStepSize = CurrentTimeStepSize / 2
					IF ( CurrentTimeStepSize < MinTimeStepSize ) CurrentTimeStepSize = MinTimeStepSize
					WRITE( Message, * ), 'Adaptive Timestepping: Time step set to ', CurrentTimeStepSize
					CALL Info('AdaptiveSolver', Message)
				ELSE
					CALL Warn('AdaptiveSolver','Cannot set smaller time step!')
					CALL Warn('AdaptiveSolver','Minimum time step might be too big to reach convergence.')
					CALL Warn('AdaptiveSolver','Continuing anyway...')
				END IF
				! reset KeepTimeStepSize
				KeepTimeStepSize = AdaptiveKeepSmallest
			END IF
		END IF
			
		! save TimeStepSize
		Var => VariableGet( Model % Variables, 'AdaptiveTimeStepSize' )
		IF(.NOT. ASSOCIATED(Var)) CALL FATAL( 'AdaptiveSolver', 'No memory allocated for AdaptiveTimeStepSize')       
		Var % Values(1) = CurrentTimeStepSize
					
		! Save SaveResults variable
		Var => VariableGet( Model % Variables, 'SaveResults' )
		IF(.NOT. ASSOCIATED(Var)) CALL FATAL( 'AdaptiveSolver', 'No memory allocated for SaveResults')       
		! It seems impossible to save an Integer, so we convert it to Real
		SaveResults_real = REAL(SaveResults)
		Var % Values(1) = SaveResults_real

		IF ( .NOT. StopSimulation ) THEN
			WRITE(Message, *)  ' SaveResults', SaveResults
		ELSE
			WRITE(Message, *)  'This was the last timestep! Preparing end of Simulation...'
		END IF
	END IF


	CALL Info( 'AdaptiveSolver', Message, level=4)
	CALL Info( 'AdaptiveSolver', '*************************************************' )
	
	IF ( .NOT. Visited ) Visited = .TRUE.
	
END SUBROUTINE AdaptiveSolver

In ResultoutputSolver, I have added a piece of code which gets the SaveResults variable and returns without saving if it is zero:

Code: Select all

  ! Check SaveResults variable from Adaptive solver, don't save if it is equal to 0
  IF( ListGetLogical(Params,'Check SaveResults Variable',Found) ) THEN
     Var => VariableGet( Model % Variables, 'SaveResults' )
     IF(.NOT. ASSOCIATED(Var)) THEN
        CALL Warn( 'ResultOutputSolver', 'SaveResults not found' )
        SaveResults = 1.0_dp
     ELSE 
        SaveResults = Var % Values(1)
     END IF
     IF( SaveResults < 0.1 ) THEN
       CALL Info( 'ResultOutputSolver', 'Nothing to save...')
       CALL Info( 'ResultOutputSolver', '-------------------------------------')
       RETURN
    END IF
  END IF  

Now in the sif file, the following is needed:

Code: Select all

Simulation
...
  Timestep intervals = 10000
  Timestep Size = Equals AdaptiveTimeStepSize
...
End
...
Solver 1
  Equation = Adaptive Solver
  Convergence Adaptive Timestepping = Logical True
  Convergence Adaptive Min Timestep = Real 0.1
  Convergence Adaptive Keep Smallest = Integer 5
  Duration = Real 1200
  Variable = -dofs 1 -global AdaptiveTimeStepSize
  Procedure = "AdaptiveSolver" "AdaptiveSolver"
  Adaptive TimeStep Size = 5
  Min Data Saving Interval = Real 1
  Exec Solver = After Timestep
End

Solver 2
  Equation = Result Output
  Procedure = "ResultOutputSolve" "ResultOutputSolver"
...
  Check SaveResults Variable = Logical True
  Exec Solver = After Timestep
End

I think that's all, but I have not tested it. My actual implementation which works is almost identical.
I hope it will be useful.

Question to the Elmer team: What do you think about implementing this functionality as an extension to the existing adaptive timestepping in Elmersolver?

Matthias

EDIT (7 Nov 2012): Removed some confusing remainders from full solver, added initialization for "CurrentTimeStepSize"

[nguyent]

Wow, thanks! It would be excellent if this functionality were integrated into ElmerSolver... in the mean time, I did a compile check and got the following:

Code: Select all

      IF ( IncState .OR. ( .NOT. Visited ) ) THEN   
                   1
Error: Symbol 'incstate' at (1) has no IMPLICIT type
mzenkersolver.f90:229.53:

         WRITE(Message, *)  'Next Step:  Loop ', Loop, ' State ',State, ' SaveR
                                                     1
Error: Symbol 'loop' at (1) has no IMPLICIT type
mzenkersolver.f90:229.70:

         WRITE(Message, *)  'Next Step:  Loop ', Loop, ' State ',State, ' SaveR
                                                                      1
Error: Symbol 'state' at (1) has no IMPLICIT type

Thank you very much for sharing. I am starting to look at your code to try to understand it. Are IncState, Loop, and State variables part of the adaptive timestepper or part of your other code?

[mzenker]

Hi,

IncState, State and Loop are part of the other code. You can replace

Code: Select all

IF ( IncState .OR. ( .NOT. Visited ) ) THEN

by

Code: Select all

IF ( .NOT. Visited ) THEN

and remove Loop and State from the WRITE statements (and also from the comments).

Matthias

[Franz Pichler]

Hi there,

i would like to stress my second question once again, hoping to start a discussion about that.
The way the error is calculated.
Right now Elmer is taking the norms of the solutions and taking the difference between them as the error measure.

Wouldn't it be better to use the norm of the difference?
best regards Franz

[Franz Pichler]

Hi, this one is for mzenker,

i had a look at your solver from abovew. looks really nice.
i think i will use it too, so thank you for sharing it.

just a short remark:
in line 114 you use CurrentTimeStepSize uninitialized in the first run, i think.
could that be true?
mybe your full solver hasn't gotr that problem.
best regards
Franz

[mzenker]

Hi Franz,

you are right, I have corrected that in the code. The missing initialization is also in my full solver and has not been a problem so far, maybe Fortran automatically initializes to zero. But it is safer not to rely to that, of course...
Thank you for pointing that out!

Matthias