• From: Gothenburg, Sweden
  • Registered: 08-09-2009
  • Posts: 1,419

Topic: SUPG has broken time integration

I have been working on getting SUPG to properly use a primary field (dofdistributedprimaryfield).

I think it's vitally important that we are very clear on how time is discretized and what we are solving for.
After quite a lot of digging, I've convinced myself that the SUPG solver is trying to do

v = v_{n} * alpha + v_{n-1} * (1-alpha)    (VM_Intermediate)
p = p_{n}                                   (VM__Total)
a = (v_{n} - v_{n_1})/dt                  (VM_Velocity (velocity of velocities -> acc))

The iterations solve for a and dp. OK, straight forward, and, in my rewrite of the code, where I only store v_{n}, v_{n-1}, p_{n} (the other value modes are computed by the primary field, which now carries responsibility for time integration)

But entering the second time-step (of scctest01.in) I see that the SUPG doesn't save the v_{n} value, but the intermediate, v, which it now treats as the previous value. So, in sorts, you get an intermediate of intermediates. And well, frankly, It's almost impossible to follow the code paths in supg.C.
Based on the values, and side-by-side comparisons, I think the second iteration becomes

v = v_{n} * alpha + (v_{n-1} * alpha + v_{n-2} * (1-alpha)) * (1-alpha)

This doensn't just affect convergence, it affects the solution.

This all stems from not having a clear enough distinction between VM_Total and VM_Intermediate, which, consdering VM_Intermediate wasn't even available when SUPG was written, shouldn't be that surprising.

I have a working SUPG solver, which, as far as I can see, does a consistent (correct) time integration. It also simplifies the code significantly, and is what is needed to prepare for the changes to BC/IC/dof handling as we have discussed before.
It's also significantly simpler, with a 200-300 lines less code, with another 100 that can be removed after a few other things have been fixed.

  • From: Gothenburg, Sweden
  • Registered: 08-09-2009
  • Posts: 1,419

Re: SUPG has broken time integration

I pushed my version on a new branch "fix_supg"

https://github.com/oofem/oofem/tree/fix_supg

I would appreciate someones input on this.

  • Registered: 25-06-2013
  • Posts: 241

Re: SUPG has broken time integration

Hi Mikael,
the supg solver is based on this paper: https://www.sciencedirect.com/science/a … 608701534, see pages 15-19.
It seems to me that there is no calculation of the first of your equations, i.e., VM_Intermediate. Can you point me to this in the code?

  • From: Gothenburg, Sweden
  • Registered: 08-09-2009
  • Posts: 1,419

Re: SUPG has broken time integration

Disclaimer: I found a serious bug (which, despite my error of having minus instead of plus, didn't actually change the results that much).
It's fixed, and pushed, already.

Hi Martin, (unfortunately, I'll have to borrow the physical book from my uni library, and I'm on vacation and won't have access to it for some time).

The idea is that the primaryfield takes is responsible of the time integration (and determines on what values needs to be stored and how to compute the respective value modes), so the value mode calculations are baked into:
DofDistributedPrimaryField::giveUnknownValue
https://github.com/oofem/oofem/blob/fix … ield.C#L60

So, basically, the idea is that v_{n-1} is known from before (or from IC), and after each iteration, update  v_{n} = v_{n-1} + a*dt in the field.
https://github.com/oofem/oofem/blob/fix … upg.C#L380
Then, the elements can request VM_Intermediate (perhaps I should have called this VM_TotalIntrinsic?) which the field can compute:
v = v_{n} * alpha + v_{n-1} * (1-alpha)
for which the equilibrium is computed for.
What i believe was the error in the old code stems from the fact that VM_Total is used for for v and v_{n}, making it easy to get it wrong.

  • From: Gothenburg, Sweden
  • Registered: 08-09-2009
  • Posts: 1,419

Re: SUPG has broken time integration

Hi Martin, I'll reply here just to keep a public record. I checked the text, and I didn't see anything that specified something obscure w.r.t. time integration, so I still think I'm correct.

Summary:
The problem with the old code is that "prevSolutionVector" doesn't contain v_{n-1}, it contains the previous v.

I'm going through the scctest01.in test for 2 time-steps (which is all that's needed to see the issue). alpha is 0.5 here.
I emphasize that v != v_1


Initial condition is:
v_0 = 0

First step, time increment = 1.754642e-03  final acceleration value after the NR solver loop is

a = [-8.876e+00   2.344e+00   1.443e+01   ...]
v = [-7.787e-03   2.057e-03   1.266e-02   ...]  ( which is equal to alpha*(a*dt) + (1-alpha)*v0 )

so far so good, the result is thus that v1 should be equal to:

v1 = v0 + a * dt = (v - (1-alpha)*v0)/alpha = [-1.557e-02   4.113e-03   2.532e-02 ...]

but, when we proceed to step 2, we still have "prevSolutionVector" contains the old intrisic value, not v_1, like it should

prevSolutionVector = [-7.787e-03   2.057e-03   1.266e-02   ...]

Currently, the scctest01 test gets me

186: Check failed in scctest01.out: tstep 20, element 3, gpnum 1, ist 43, component 1:
186: value is 2.04227449e-01, but should be 1.95192016e-01 ( error is 9.035433e-03 but tolerance is 1.000000e-06 )

so the difference is not drastic.