Missing acceleration when restart from existing solution

[stoykov]

Hello Elmer team,

In some of my old posts I described about a problem that in 2nd order time-depenent systems, the solution obtained by restarting from previous solution differs from solution obtained without restarting. I noticed this when I used it for elasticity problems, such as StressSolve, ElasticSolve and Smitc.

After some checks in the codes, I may conclude that this difference appears because when using an existence solution, the acceleration of that existence solution is not taken into account i.e. it is considered to be zero. If one looks in the Bossak method (eq. 5.7 from ElmerSolver Manual), one may see that the solution obtained at time i+1 is obtained by considering the solution, velocity and acceleration at time i. But when one wants to restart from solution at time i, the acceleration at time i is considered to zero.

I think that this is very important and I would like to suggest to the developers to consider this problem. One possible solution is to obtain the acceleration, in the cases when the user wants to restart from existing solution. Since, the displacement and the velocity are known, the acceleration might be obtained by replacing the displacement and the velocity in the equation of motion and solving an algebraic system of equations for the acceleration vector. Maybe this can be done in the SUBROUTINE LoadRestartFile at ModelDescription.src.

I would be grateful if you consider my note and update the source!

Best regards,
Stan

[mzenker]

Hi,

since the Elmer team has limited manpower, you might also try to implement a solution to your problem yourself. What you can do is to check if it is possible to modify an existing solver module (or write a new one) such that the desired effect is achieved. Solver modules can be very easily compiled using the Elmerf90 script. Compiling the core Elmersolver is also possible, but more difficult (at least under Windoze).

Matthias

[stoykov]

Yes, I am thinking to implement the modifications and write to the developers, so if they agree they might upload the source.

[rimple_sandhu]

I am facing the same issue. I am getting numerical damping and bias in frequency when restarting the transient simulation at each time step, when compared to the response obtained without stopping. Did some one found the solution to this issue?

[raback]

Hi

I think there is indeed a problem.

Now by default Elmer saves just the current displacement, and not any previous values (or velocity and acceleration terms). If the time-stepping scheme requires more than one previous value this unfortunately cannot be consistent with running everything in one sweep. So this means all 2nd order PDEs with time, and all higher order schemes for 1st order PDEs. In the latter case the code will always automatically start with 1st order scheme so having a small 1st timestep could be almost a sufficient remedy. For 2nd order PDEs this is not an option.

So, to have everything work consistently. One could do the following.

1) In the code decleare the previous values as variables that would then be automatically saved/loaded. The following code in MainUtils.F90 does this for one the velocity component of PrevValues for 2nd order schemes:

Code: Select all

      IF( ListGetLogical( Solver % Values,'Calculate Velocity',Found) ..
          ...
          Component => Solver % Variable % PrevValues(:,1)
          str = TRIM( Solver % Variable % Name ) // ' Velocity'
          CALL VariableAddVector( Solver % Mesh % Variables, Solver % Mesh, Solver, &
              str, Solver % Variable % Dofs, Component, Solver % Variable % Perm, &
              Secondary = .TRUE. )
      END IF

Something like this could work:

Code: Select all

      IF( ListGetLogical( Solver % Values,'Save Previous Values',Found) ..
          DO i=SIZE( Solver % Variable % PrevValues,2)
            Component => Solver % Variable % PrevValues(:,i)
            str = TRIM( Solver % Variable % Name ) // ' PrevValue'//TRIM(I2S(i))
            CALL VariableAddVector( Solver % Mesh % Variables, Solver % Mesh, Solver, &
              str, Solver % Variable % Dofs, Component, Solver % Variable % Perm, &
              Secondary = .TRUE. )
           END DO
      END IF

2) In addition you should then utilize the loaded values in the timestepping methods. So the "Solver % Order" should define the order always, and not the "Solver % DoneTime". Maybe adding some intelligence to MainUtils/AddEquationSolution could do the trick.

Ok, I'm actually almost there now. I could give it a short myself

-Peter

[raback]

Hi

Ok, this is now implemented and available in branch "devel". To use it say "Transient Restart = Logical True" at the Solvers where you need to apply the higher order routines in. Unfortunately the way I did it may result to bloated files because the previous values will be default now be saved at other solvers too.

There is a test case: ElasticBeamRestart, which demonstrates how exactly the same result is obtained using 100 or 50+50 timesteps.

Have fun!

-Peter

[rimple_sandhu]

Thanks peter. It works now. I just want to know what was changed in the source code. Is it that you now save displacements from several of previous iterations, or you use different time-descritization strategy for 2nd order PDEs. Also, with this implementation, I don't think we need to use 'Calculate Velocity = Logical True' in the simulation part of the sif file.

Nevertheless, thanks for all the effort.

-Rimple

[raback]

Hi

What was changed is that all previous stesps, or more generally all work vectors related to given time-stepping scheme are saved. These are actually named to be fields such that the standard saving automatically treats them.

-Peter