[t.teschner]

I have some issues implementing the dual time stepping procedure and was hoping some experienced eyes can help me out. Suppose I want to solve the advection equation using a dual time stepping procedure, then I have

,

where and are the pseudo and real time, respectively. Using first order (for simplicity) for all the derivatives, I get

.

I define the right-hand side (RHS) as



and then can rewrite the equation in explicit time integration form as

.

Once , where is some small (convergence number), the pseudo time derivative vanishes and i obtain the solution at the next time step. The pseudo code for my example follows below

while()
   while( and )
   
      Do i=2 to n
         
      End

      Do i=2 to n
         
      End

      calculate residual
      set iteration=iteration+1
      

   End

   
   

End

As long as I have the above code works fine, but as soon as I use different time steps, say , then my solution is off (by a factor of 2, in terms of propagated distance). It is probably an easy fix, but at the moment I am not quite sure what I am doing wrong.

[FMDenaro]

I never worked with such a procedure, thus I have no experience about the advantage in doing that. However, I see that for dt=dtau you get the increment without
u^{n+1,m}

so I wonder what about the iteration over m...

[t.teschner]

well, as far as i understand, the method was introduced to allow usage of fast time marching algorithms (implicit schemes, multigrid etc.) which need not to be accurate but give a steady state solution fast. the time accurate solution is recovered with the second timestep which can follow its on stability criterion.
well, so far the theory, but all method derived from the artificial compressibility method in incompressible flows (which i am working on) suffer from the inherent pseudo time marching so the only way, really, to get an unsteady simulation, is to use the dual time stepping procedure.
so i though i start with a simple advection equation but yeah, if than vanishes. I suspect the problem to be somewhere in the update of the variables, though.

[sbaffini]

The DTS formulation refers to the fact that you have 2 time derivative terms, one of which is in pseudo-time, the other in actual time so, formally, by construction, it is used for unsteady computations. By extension, it is also possible for steady states, where the true time derivative term is assumed to be 0.

Now, what the problem is for such unsteady computations? That the explicit stability limit to advance them might require much lower time steps than accuracy would, possibly only because of few problematic cells.

The idea of the DTS is thus to introduce a pseudo-time derivative, discretized with a local time step (on a per cell basis) based on the local stability limit, and use such pseudo-time to advance explicilty the solution until it doesn't change anymore in pseudo-time, which means that the pseudo-time derivative has disappeared and the remaining terms (i.e., the original equation) are fully balanced.

Now, what should be clear is that, really, it only makes sense for implicit discretizations in time or implicit steady states. Otherwise, your physical time step, by construction, has to be, for each cell, equal to the lowest allowed pseudo-time step in the grid.

So, to start with, as your discretization is explicit, both in time and pseudo-time, both time steps have to respect the stability limits. For your simple case it means:





Do the both of them respect these constraints?

Second question, how do you determine convergence within pseudo-time? Are you sure you are not simply using all the allowed iterations without actually converging? (Also note that, typically, the condition used in the check is an AND, not an OR)

[t.teschner]

thanks for spotting that, in my code i have indeed a "and" statement, not or. and yes, i specify the timestep according to the stability criterion stated in your post. as i was saying before, when dealing with the artificial compressibility method, dual time stepping is the only way to achieve unsteady simulations, so we don't use it for the larger timestep reason as in compressible flow and hence an explicit method will work just fine for me.
i set my maximum iterations to 100 and found convergence for eps=1e-8 to be around 30 iterations, i played around with the maximum iteration number and eps but usually the solution converged in pseudo time before the maximum iteration count was achieved. the residual is calculated based on the L0 norm and i just take the maximum difference of u between m+1 and m (and normalise that difference so that the residual is 1 at the first iteration). if it helps, i put the matlab code below (i use k here instead of m to indicate the pseudo time level)

%% the program solves the advection equation with dual-time stepping
clc;
clear;
close all;

%% initial parameters

% problem specific parameters
a     = 1.0;   % wave speed
n     = 1001;  % number of elements in x
CFL_r = 0.8;   % CFL, real time
CFL_p = 0.4;   % CFL, pseudo time
t_end = 5.0;   % compute until

% convergence criterion in pseudo time
eps   = 1e-8;  % convergence parameter in pseudo time
imax  = 100;   % maximum iterations in pseudo time

% calculated parameters
dx    = 20/(n-1);   % spacing
dt_p  = CFL_p*dx/a; % delta t(au), pseudo time
dt_r  = CFL_r*dx/a; % delta t, real time

% solution arrays
un1k1 = zeros(n,1); % u^{n+1,k+1}
un1k  = zeros(n,1); % u^{n+1,k}
un    = zeros(n,1); % u^{n}
unm1  = zeros(n,1); % u^{n-1}
uana  = zeros(n,1); % u analytic
RHS   = zeros(n,1); % right-hand side vector

% space vector
x = linspace(-10,10,n); % x vector

% initialise solution
for i=1:n
    un1k1(i) = 0.5*exp(-x(i)^2);
    un1k(i)  = 0.5*exp(-x(i)^2);
    un(i)    = 0.5*exp(-x(i)^2);
    unm1(i)  = 0.5*exp(-x(i)^2);
end

% compute analytic solution for comparison
for i=1:n
    uana(i) = 0.5*exp(-(x(i)-t_end)^2);
end

%% time loop

% loop in real time (outer loop)
t_r = 0.0;
while(t_r<t_end)
    
    % copy values
    unm1 = un;
    un   = un1k;
   
    % loop in pseudo time (inner loop)
    res  = 1.0;
    iter = 1; 
    while(eps<res && iter<imax)
        
        % copy values
        un1k = un1k1;
        
        % loop over space
        for i=2:n            
            % assemble right-hand side
            temp = getRHS(un1k, un, unm1, dt_r, dx, a, i);
            RHS(i) = temp;            
        end
        
         % time integration in pseudo time
        for i=2:n            
            un1k1(i) = un1k(i)+dt_p*RHS(i);
        end
        
        % calculate residual
        res = abs(max(un1k1-un1k));
        
        if (iter == 1)
            if (res ~= 0.0)
                resnorm = res;
            else
                resnorm = 1.0;
            end
        end
        
        % normalise residuals
        res = res/resnorm;
        
        % print to screen
        if(mod(iter,100) == 0)
            fprintf('   iter: %d   residual: %12.5e\n', iter, res)
        end
        
        % increase iteration counter
        iter = iter + 1;
        
    end
    
    fprintf('time: %12.2e   iterations in pseudo time: %d\n',t_r,iter);
    t_r = t_r + dt_r;
end

%% plot solution
plot (x,un1k1,x,uana)
legend('numerical','analytical','location','northwest')
ylim([0 0.5])

and the function call

function [temp] = getRHS(un1k, un, unm1, dt_r, dx, a, i)
    
    temp = (-1.0)*(3*un1k(i)-4*un(i)+unm1(i))/dt_r - a*(un(i)-un(i-1))/dx;
    
end

with the code above you'll get the following output

[sbaffini]

I can'try your code now. But I suspect that the following 3 points are creating troubles:

1) res = abs(max(..., where it should be max(abs(...

2) you are not using the 1st order explicit euler in real time as claimed, but a 2nd order backward discretization. Still, you have it wrong, because you need to divide that by 2dt

3) as per point 2 above, your discretized physical unsteady term is actually defined at n+1, which means you would need to go implicit. I don't know if staying explicit here just means reverting to 1st order

[sbaffini]

Did the check in 2. It seems to solve your problem. Check in 1 does not seem to be relevant for this case.

[t.teschner]

it was about point 2 ... thanks, problem is solved now (forgot to devide by 2dt). i just used first order in the opening post for simplicity but 2nd order in the code. and there does not seem to be a problem with having a double explicit scheme, as stated before, for artificial compressibility based methods its the only way to introduce real unsteady behaviour so both pseudo and real time can be explicit. in this case i just use the advection equation to simplify the problem. but thanks again for spotting the bug