The adaptative VSVO Method

This section will present an algorithm created to adapt both the order and the time step to efficiently and more precisely solve EDPs using a backward Euler FEM code. This algorithm is presented by V. P. DeCaria in the dissertation Variable Stepsize, Variable Order Methods for Partial Differential Equations.

1. Algorithm Variable Stepsize, variable Order 1 and 2

First, let’s give the method for the initial problem :

\[ y'(t) = f(t,y(t))\text{, for }t>0\text{ and }y(0)=0\]

We write \(\Delta t_n\) as the \(n^{th}\) time stepsize. We have also \(t^{n+1}=t^n + \Delta t_n\) and \(y^n\) an approximation for \(y(t_n)\). \(TOL\) is the tolerance chosen by the user on the allowable error by step.

1.1. Step 1 : Backward Euler

\[\frac{y_{(1)}^{n+1}-y^n}{\Delta t_n}=f(t_{n+1},y^{n+1}_{(1)})\]

1.1.1. Implementation in FreeFEM

We start a While that will stop when the current time \(T\) reaches the final time \(T\) :

Initialization of the Loop and Step 1 in the FreeFEM file
while (t<=T){

     // previous solution are saved in case of
     // returning to previous time step
     ysave = yppp;
     yppp = ypp;
     ypp = uold;
     uold = u; //equivalent to u^{n-1} = u^n

     // STEP 1 : backward euler
     // Solve
     mqs;

1.2. Step 2 : Time Filter

\[y_{(2)}^{n+1}=y_{(1)}^{n+1}-\frac{\omega_n}{2\omega_n +1}(y_{(1)}^{n+1}-\left(1+\omega_n)y^n+\omega_n y^{n-1}\right)\]

with \(\omega_n=\frac{\Delta t_n}{\Delta t_{n-1}}\) the filtering weight.

1.2.1. Implementation in FreeFEM

Step 2 in the FreeFEM file
while (t<=T){
    //[...]

    // STEP 2 : time filtering
     wn = dt/dtp;
     wp = dtp/dtpp;

     y1=u; // y1(n+1)
     y2=y1 - wn/(2*wn+1)*(y1-(1+wn)*uold+wn*ypp);  // y2(n+1)

1.3. Step 3: Estimate Error in \(y_{(1)}^{n+1}\) and \(y_{(2)}^{n+1}\)

We compute an estimation of the local errors in this step. \(Est_1\) provides an estimation of the local error of the first order approximation and \(Est_2\) of the second order approximation.

\[EST_1=y_{(2)}^{n+1}-y_{(1)}^{n+1}\]
\[EST_2=\frac{\omega_{n-1}\omega_n(1+\omega_n)}{1+2\omega_n+\omega_{n-1}(1+4\omega_n+3\omega_n^2}\left( y_{(2)}^{n+1}\\ - \dfrac{(1+\omega_n)(1+\omega_{n-1}(1+\omega_n))}{1+\omega_{n-1}} y^n +\omega_n(1+\omega_{n-1}(1+\omega_n))y^{n-1}\\ -\dfrac{\omega_{n-1}^2\omega_n(1+\omega_n)}{1+\omega_{n-1}}y^{n-2}\right)\]

1.3.1. Implementation in FreeFEM

Step 3 in the FreeFEM file
while (t<=T){
    //[...]

    // STEP 3 : estimate error
     est1=y2-y1;
     est2=wp*wn*(1+wn)/(1+2*wn+wp*(1+4*wn+3*wn^2))*(
               y2
          -(1+wn)*(1+wp*(1+wn))/(1+wp)*uold
          +wn*(1-wp*(1+wn))*ypp
          -wp^2*wn*(1+wn)/(1+wp)*yppp
          );

     est1l2 = sqrt( int2d(Th)((est1)^2));
     est2l2 = sqrt( int2d(Th)((est2)^2));

    // safety net in order to not divide by 0 later
     if (est1l2==0) est1l2=1e-20;
     if (est2l2==0) est2l2=1e-20;

1.4. Step 4: Check if tolerance is satisfied

If \(||EST_1||<TOL\) or \(||EST_2||<TOL\), at least one approximation is correct, got to Step 5a. Otherwise the step is rejected, got to Step 5b

1.4.1. Implementation in FreeFEM

Step 4 in the FreeFEM file
while (t<=T){
    //[...]

    // STEP 4 : Check Tol
     if (( est1l2 > tol ) && (est2l2 > tol)){
          // STEP 5b : Neither approximation satisfy TOL

     }
     else if (( est1l2 <= tol ) || (est2l2 <= tol)){
          // STEP 5a : At least one approximation satisfy TOL
     }

1.5. Step 5a: At least one approximation satisfy TOL

If both approximation are acceptable, set

\[\Delta^{(1)}=0.9\Delta t_n\left(\frac{TOL}{||EST_1||}\right)^{1/2},\qquad\Delta^{(2)}=0.9\Delta t_n\left(\frac{TOL}{||EST_2||}\right)^{1/3}\]

And set

\[ i=arg\max_{i\in\{1,2\}}\Delta t^{(i)},\quad \Delta t_{n+1}=\Delta t^{(i)},\\ \quad t^{n+2}=t^{n+1}+\Delta t_{n+1},\quad y^{n+1}=y_{(i)}^{n+1}\]

(The approximation which have the largest estimated \(\Delta t\) is chosen to advance the solution.)

If only \(y^{(1)}\) (resp. \(y^{(2)}\)) satisfies \(TOL\), set \(\Delta t_{n+1}=\Delta t^{(1)}\) (resp. \(\Delta t^{(2)}\)), and \(y^{n+1}=y^{n+1}_(1)\) (resp. \(y^{n+1}=y^{n+1}_(2)\)). Proceed to Step 1 to calculate \(y^{n+2}\).

1.5.1. Implementation in FreeFEM

Step 5a in the FreeFEM file
while (t<=T){
    //[...]

    else if (( est1l2 <= tol ) || (est2l2 <= tol)){
          // STEP 5a : At least one approx is accepted

          dt1 = 0.9*dt*sqrt(tol/est1l2);
          dt2 = 0.9*dt*(tol/est2l2)^(1/3.);

          dtpp = dtp;
          dtp = dt;
          if (( est1l2 <= tol ) && (est2l2 <= tol)){
               dt = max(dt1, dt2);
          }
          else if( est1l2 <= tol ){
               dt = dt1;
          }
          else{
               dt = dt2;
          }

          if (dt == dt1){
               u = y1;
          }
          else{
               u = y2;
          }
          tp = t;

          t = t + dt;
     }

1.6. Step 5b: Neither Approximation satisfy TOL

Set

\[ \Delta^{(1)}=0.7\Delta t_n\left(\frac{TOL}{||EST_1||}\right)^{1/2},\qquad\Delta^{(2)}=0.7\Delta t_n\left(\frac{TOL}{||EST_2||}\right)^{1/3}\]

Set

\[ i=arg\max_{i\in\{1,2\}}\Delta t^{(i)},\quad \Delta t_{n}=\Delta t^{(i)}, \quad t^{n+1}=t^{n}+\Delta t_{n}\]

Return to Step 1 and try again.

1.6.1. Implementation in FreeFEM

Step 5b in the FreeFEM file
while (t<=T){
    //[...]

    if (( est1l2 > tol ) && (est2l2 > tol)){
        // STEP 5b : Neither approximation satisfy TOL

        dtsave = dtpp;
        dtpp = dtp;
        dtp = dt;

        dt1 = 0.7*dt*sqrt(tol/est1l2);
        dt2 = 0.7*dt*(tol/est2l2)^(1/3.);

        dt = max(dt1, dt2);

        t = tp+dt;

        dtp = dtpp;
        dtpp = dtsave;

        u = uold;
        uold = ypp;
        ypp = yppp;
        yppp = ysave;
    }
    //[...]
}

1.7. Additionnal Implementations

This was the basic implementation of the VSVO Algorithm, but there are several standard features that are missing.

1.7.1. Maximum stepsize

We impose a maximum time step in order to have reasonable stepsize and to prevent at a certain level to go to the Step 5b because the timestep was really large. This feature appears only in the Step 5a as this is the step where the stepsize can increase.

Implementation of the maximum stepsize in the FreeFEM file
real dtmax=1;

while (t<=T){
    //[...]

    else if (( est1l2 <= tol ) || (est2l2 <= tol)){
          // STEP 5a : At least one approx is accepted

          dt1 = 0.9*dt*sqrt(tol/est1l2);
          dt2 = 0.9*dt*(tol/est2l2)^(1/3.);

          dtpp = dtp;
          dtp = dt;
          if (( est1l2 <= tol ) && (est2l2 <= tol)){
               dt = max(dt1, dt2);
          }
          else if( est1l2 <= tol ){
               dt = dt1;
          }
          else{
               dt = dt2;
          }

          if (dt == dt1){
               u = y1;
          }
          else{
               u = y2;
          }
          // Check if dt>dtmax before t is changed
          if (dt>dtmax) dt=dtmax;

          tp = t;

          t = t + dt;
     }

1.7.2. Minimum stepsize

We impose a minimum time step in order to avoid the algorithm to be too long. If the timestep is always decreasing and the estimation is not acceptable each time, then it’s a reasonable choice to go to the next time and start the process once again. This feature appears only in the Step 5b as this is the step where the stepsize can decrease.

Implementation of the minimum stepsize in the FreeFEM file
real dtmax=1;

while (t<=T){
    //[...]

    if (( est1l2 > tol ) && (est2l2 > tol)){
          // STEP 5b : Neither approximation satisfy TOL

          dtsave = dtpp;
          dtpp = dtp;
          dtp = dt;

          dt1 = 0.7*dt*sqrt(tol/est1l2);
          dt2 = 0.7*dt*(tol/est2l2)^(1/3.);

          dt = max(dt1, dt2);

          // if dt<dtmin force it
          if (dt<dtmin){
               dt=dtmin;

               tp = t;

               t = t + dt;

          }
          // else try again last t with new dt
          else {
               t = tp+dt;

               dtp = dtpp;
               dtpp = dtsave;

               u = uold;
               uold = ypp;
               ypp = yppp;
               yppp = ysave;
          }
     }

1.7.3. Forced timestep

As the stepsize is variable, we cannot be sure that \(t\) will stop at certain values that we are interested in. That’s why being able to force \(t\) to go trough specific values is a must-have in this type of algorithm.

In order to do that, we creat a list with all the values we want \(t\) to go trought. Then before updating the new \(t\) in Step 5a and Step 5b, we check if the current \(t\) is in the list and if it is we go to the next one. Finally after applying the stepsize to obtain a new timestep, we check if the new timestep is not higher than the current value in the list and in this case we adapt the stepsize to reach precisely the forced value.

Implementation of the forced timesteps in the FreeFEM file
int nforced=3; //number of forced values
real[int] forcedtime = [1,20,22]
int force = 0; //index of forced values

while (t<=T){
     //[...]

     // STEP 4 : Check Tol
     if (( est1l2 > tol ) && (est2l2 > tol)){
          // STEP 5b : Neither approximation satisfy TOL

          //[...]

          // if dt<dtmin force it
          if (dt<dtmin){

               dt=dtmin;
               //check if t in forcedtime
               if (force<nforced) {
                    if (t==forcedtime[force]) force = force + 1;
               }

               tp = t;
               //check if t+dt!>forcedtime
               if (force<nforced) {
                    if ((t!=forcedtime[force]) && ((t+dt)>forcedtime[force])){
                         dt = forcedtime[force]-t;
                    }
               }

               t = t + dt;
          }
          // else try again last t with new dt
          else {
               //check if t+dt!>forcedtime
               if (force<nforced) {
                    if ((tp!=forcedtime[force]) && ((tp+dt)>forcedtime[force])){
                         dt = forcedtime[force]-tp;
                    }
               }

               t = tp+dt;

               //[...]
          }
     }
     else if (( est1l2 <= tol ) || (est2l2 <= tol)){
          // STEP 5a : At least one approx is accepted

          //check if t in forcedtime
          if (force<nforced) {
               if (t==forcedtime[force]) force = force + 1;
          }

          dt1 = 0.9*dt*sqrt(tol/est1l2);
          dt2 = 0.9*dt*(tol/est2l2)^(1/3.);

          dtpp = dtp;
          dtp = dt;
          if (( est1l2 <= tol ) && (est2l2 <= tol)){
               dt = max(dt1, dt2);
          }
          else if( est1l2 <= tol ){
               dt = dt1;
          }
          else{
               dt = dt2;
          }

          if (dt == dt1){
               u = y1;
          }
          else{
               u = y2;
          }

          tp = t;
          //check if t+dt!>forcedtime
          if (force<nforced) {
               if ((t!=forcedtime[force]) && ((t+dt)>forcedtime[force])){
                    dt = forcedtime[force]-t;
               }
          }

          t = t + dt;

     }

}

2. Implementation on an example

We will use this algorithm on a simple example of Maxwell Quasi Static problem, time-dependent with the A Formulation for a geometry of a torus surrounded by air in axisymmetric coordinates.

2.1. Run the Calculation

The command line to run this case is :

    ff-mpirun -np 4 mqs_vsvo.edp -wg

2.2. Data Files

The case data file are available in Github here :

2.3. Equation

The equation is described in AxiAtheta.adoc

2.4. Results

Let’s compare the evolution of stepsize compared to the evolution of the electrical tension \(U\) as a function of time. The maximum stepsize is \(\Delta t_{max}=1\) :

dt vsvo
Time evolution of the variable stepsize compared to the electric tension

We can observe that when the tension is constant, the stepsize goes quickly to the maximum stepsize, but when \(U\) undergoes important changes, the stepsize becomes very small to be more precise at those times.

We also can see that, when compared to the analytic solution, the model using the VSVO Algorithm is more precise than the model with constant stepsize :

Table 1. L2 relative error Norm between the Intensity calculated and the analytical solution

No VSVO

VSVO

L2 relative error Norm : \(1.72 \%\)

L2 relative error Norm : \(0.66 \%\)

3. References

  • Variable Stepsize, Variable Order Methods for Partial Differential Equations\\}_, DeCaria, V. P., University of Pittsburgh, ProQuest Dissertations Publishing, 2019.22615879. }