Skip to content
Snippets Groups Projects

Adaptive time integration

Last edited by Christophe Geuzaine

Efficient transient simulations can be done with TimeLoopAdaptive. It varies the time step size according the local truncation error (LTE). Therefore a fine time resolution is only applied when needed in order to achieve an accurate result.

TimeLoopAdaptive works with the so-called predictor-corrector method. This means that before calculating the actual time step the result is predicted from the past solutions by polynomial extrapolation. This has the benefit of better convergence in case of non-linear problems because the initial guess for the iterative loop is already close to the result.

Local Truncation Error (LTE)

The LTE is the error made in one time step by the numerical integration:

\epsilon_{LTE} = ( y(t_{n+1}) - y(t_n) ) - ( y_{n+1} - y_n )

where y(t_{n+1}) is the true solution and y_{n+1} is the calculated solution of the actual time step. The LTE is internally calculated and the time step size is adjusted to keep the LTE within defined limits. Those limits are specified by a relative and absolute tolerance (reltol \epsilon_{rel} and abstol \epsilon_{abs}) for the quantities of interest (e.g. all DOFs of a system or results from a post-operation like temperature, current, etc.).

For each quantity of interest x_i an error ratio q_i is calculated:

q_i = \frac{|LTE_i|}{\epsilon_{abs} + \epsilon_{rel}|x_i|}

Error Norm

As the error ratio is calculated for each quantity of interest (e.g. one q per DOF of a given system) we need a norm in order to assess the overall error ratio. There are several norm types available for TimeLoopAdaptive:

L1 Norm

\|q\|_{L1} = \sum_{i=1}^n |q_i|

Mean L1 Norm

\|q\|_{Mean L1} = \frac 1n \sum_{i=1}^n |q_i|

L2 Norm

\|q\|_{L2} = \sqrt{\sum_{i=1}^n \left( q_i \right)^2}

Mean L2 Norm

\|q\|_{Mean L2} = \sqrt{\frac 1n \sum_{i=1}^n \left( q_i \right)^2}

Linfinity Norm

\|q\|_\infty  = \max \left( |q_i| \right)

If \|q\| \le 1 then the actual time step is valid. In the other case the calculated time step is rejected and must be calculated again with a reduced step size. For the next step the step size is adjusted to aim for \|q\| = 0.8.

For the mathematical background, please have a look at: http://qucs.sourceforge.net/tech/node23.html

Usage

A typical usage of TimeLoopAdaptive is shown in the following piece of code:

Resolution {
  { Name SimTimeLoopAdaptive;
    System {
      { Name SystemT; NameOfFormulation FormulationT; }
    }
    Operation { 
      InitSolution[SystemT]; 
      SaveSolution[SystemT];

      TimeLoopAdaptive[ tStart, tStop, dtInit, dtMin, dtMax, Method, List[Breakpoints], 
        System { { SystemT, Treltol, Tabstol, LinfNorm } } ]
      {         
        // This first part is executed for each time step
        Generate[SystemT]; 
        Solve[SystemT]; 
      }
      {
        // This second part is executed only if the time step is accepted,
        // i.e. the LTE is within the specified tolerances
        SaveSolution[SystemT];
      } 
    }
  }
}

The parameters are:

  • tStart: Start time

  • tStop: Stop time

  • dtInit: Initial time step size

  • dtMin: Minimum time step size

  • dtMax: Maximum time step size

  • Method: Integration method: Euler, Trapezoidal, Gear_2, Gear_3, ..., Gear_6;

    Gear's method correspond to the backward differentiation formulas (BDF) of order 2, ..., 6. If you like you can use the synonyms BDF_2, ..., BDF_6 instead.

    Important note: For Trapezoidal there are the same restrictions as for TimeLoopTheta! See the warning in the reference manual for TimeLoopTheta. I.e. it's not a valid method when there are time varying parts in the Galerkin statements using the standard equations (e.g. when we have a time dependent source term there or for non-linear problems). For a workaround have a look in the reference manual part of TimeLoopTheta. All other methods (Euler, Gear_2, ..., Gear_6) do not have this restriction.

  • Breakpoints: Here you can specify time points to be met during time stepping, e.g. when you expect an abrupt change of the solution due to a piece wise linear source definition.

  • System: Keyword System indicates that the quantity of interest are all DOFs of the given system. In addition or as an alternative you can specify post-operation results as quantities of interest by the keyword PostOperation (see example below).

  • System name: Here the name of the system of interest is declared. In case of a post-operation you specify here the name of the post-operation.

  • Relative tolerance: Relative tolerance for the quantity of interest.

  • Absolute tolerance: Absolute tolerance for the quantity of interest.

  • Error norm: Type of the error norm. Possible choices are: L1Norm, MeanL1Norm, L2Norm, MeanL2Norm, LinfNorm (→ see definitions above).

  • 1st set of operations:
    This set of operations are executed for each time step. Usually the solution is calculated here.

  • 2nd set of operations: This operations are executed only when the time step is valid, i.e. LTE is in the specified tolerance.

Quantities of Interest

Quantities of interest can be all degrees of freedom of one or several systems. Additionally or alternatively the quantities of interest can be the results of one or several post-operation(s). Each of them can have their own tolerances. Here is an example:

TimeLoopAdaptive[ tStart, tStop, dtInit, dtMin, dtMax, Method, List[Breakpoints], 
  System        { { MySystemA,   Areltol, Aabstol, LinfNorm   }
                  { MySystemB,   Breltol, Babstol, L1Norm     } }
  PostOperation { { Voltage,     Vreltol, Vabstol, LinfNorm   }
                  { Temperature, Treltol, Tabstol, MeanL2Norm } } ]
{         
  ...
}
{
  ...
}

Comparison of system versus post-operation

  • For non-linear or coupled problems it's recommended to specify the system(s) because then for each new time step the solution is predicted (predictor-corrector method) before the actual calculation is done. This yield in better convergence.
  • When a system is specified then all DOFs of this system are assessed, but they do not necessarily correspond to the values of the quantities of the formulation directly (e.g. in case of higher order basis functions).
  • In case of several quantities in a formulation they can be assessed individually by a specified post-operation.
  • With a post-operation it's also possible to assess different regions individually.
  • Also derived quantities can be assessed with a specified post-operation.
  • Performing a post-operation each time step can be time consuming for large problems.

Obtaining Results with Equidistant Time Steps

Sometimes it might be advantageous to have a result with equidistant time steps, e.g. for creating a movie with constant speed. There are two ways to do this:

1) During simulation with breakpoints

Function {
...
  dt = 125e-9;  // Time step size for breakpoints [s]
  Breakpoints = {(tStart+dt):tStop:dt};
}

Resolution {
...
      TimeLoopAdaptive[ tStart, tStop, dtInit, dtMin, dtMax, Method, List[Breakpoints], 
        System { { SystemT, Treltol, Tabstol, LinfNorm } } ]
      {         
        // This part is executed for each time step
        Generate[SystemT]; 
        Solve[SystemT]; 
      }
      {
        // Save the solution only at the equidistant breakpoints:
        // If we are at a breakpoint $Breakpoint gives its numer.
        // Otherwise $Breakpoint = -1
        If[ $Breakpoint >= 0] { SaveSolution[SystemT]; }
      } 
...
}

2) After simulation in the post-operation

There is a new feature, called ResampleTime. It interpolates the result at specified time points via a cubic spline interpolation. This feature is independent of the time loop type.

PostOperation {      
  { Name T_resampled; NameOfPostProcessing The; ResampleTime[tStart, tStop, dtResample];
    Operation {
      Print[ T, OnElementsOf Vol_The , File "Result_T_map.pos"];
    }
  }  
}

General Hints

  • If you specify several systems and/or post-operations then usually only the largest \|q\| is printed. If the verbosity level is set higher than 5 then \|q\| is printed for each system or post-operation.
  • If you start GetDP within Gmsh then time, time step size and the largest \|q\| (LTEmaxErrorRatio) can be plotted nicely.
  • When TimeLoopAdaptive aborts the first time step because the minimum step size is reached and LTEmaxErrorRatio is still quite high then there is usually a jump in a quantity of interest. E.g. this is usually the case when the initial solution is not really a solution of the system or if it is not completely initialized (This can easily happen e.g. when lumped element circuits are used). You can overcome this problem by using TimeLoopTheta for just the first time step(s) and applying TimeLoopAdaptive afterwards, as shown in this example:
Resolution {
  { Name SimTimeLoopAdaptive;
    System {
      { Name SystemT; NameOfFormulation FormulationT; }
    }
    Operation { 
      InitSolution[SystemT]; 
      SaveSolution[SystemT];
      TimeLoopTheta { Time0 tStart; TimeMax tStart+dtInit; DTime dtInit; Theta 1.0;
        Operation { 
          Generate[SystemT];
          Solve[SystemT]; 
          SaveSolution[SystemT];
        } 
      }
      TimeLoopAdaptive[ tStart+dtInit, tStop, dtInit, dtMin, dtMax, Method, List[Breakpoints], 
        System { { SystemT, Treltol, Tabstol, LinfNorm } } ]
      {         
        // This part is executed for each time step
        Generate[SystemT]; 
        Solve[SystemT]; 
      }
      {
        // This part is executed only if the time step is accepted,
        // i.e. the LTE is within the specified tolerances
        SaveSolution[SystemT];
      } 
    }
  }
}

Complete Example

Here a complete example is given for a thermal problem. We have an object consisting of two materials. On top a uniform heat source is applied.

Geometry: Thermal_Sim.geo

Problem: Thermal_Sim.pro

Gmsh_Simulation

*Initially written by @AsamMich.

Comments

    Please register or sign in to add a comment.