CCAPS (all flavours) - RTI bug (and avoid hard-coding of fixed eta for MG boundary)

[jothurgood]

CCAPS (all flavours) - RTI bug (and avoid hard-coding of fixed eta for MG boundary)

[JOThurgood]

Regarding the apparent bug:

• Noticed that for the mg eta bc at top set to 1/7 that the conditioning is much better, in the sense that typical number of V cycles is reduced from O(100) to O(20). Since ideally each V cycle should reduce an order of magnitude on the L2, this is closer to ideal.

• The better conditioning is noticed at two resolutions (both nx=16 and nx=66, with ny=4*nx in both).

• So better conditioning => ~5x faster code!

• However, at higher resolutions, the 1./7 ymax condition will begin to exhibit nonconvergence in MG (despite doing better up to that stage).

• Also noticed that for lower resolutions (nx=16, ny= 4*nx) that don't seem to get the crash. Might hint that artificial diffusion mitigates the problem?

Thoughts

• Maybe it just isn't really physical to have no slip at the top, especially in the absence of any sort of diffusion. Perhaps either ZG or slip (uy =/0) on top would avoid the issue. Also wonder if no slip is really appropriate for bottom.

[JOThurgood]

Also, regarding bug:

• With "wrong" eta fixing (=7) you get white on the rho plots (not on the colour table), but its okay in the "correct" runs (at least up until a crash)

• Don't know if the white corresponds to new or increased mass, or if it is a vacuum - check this. Do know the grid-integrated rho is about comparable in each simulation, though.

[JOThurgood]

• bug persists for other RTI (e.g., unperturbed interface with high enough resolution goes RTI unstable), if sufficiently high resolution.
• Again, "wrong" boundaries give it poor conditioning but it seems to be solvable (or at least, solvable for longer than otherwise)

[JOThurgood]

Don't see why I was using fixed eta bc for MG when, since it basically corresponds to rho (or 1/rho), zero gradient is probably more appropriate in the first instance? (e.g. in the godonov solver, bc on rho are taken as ZG if using noslip)

Experimented with the "fixed but wrong eta" versus zero gradient in the main CCAPS code

• ZG has much better V-cycle conditioning
• ZG shows no rho error (white in plot) at top
• otherwise very similar sln
• At higher nx (nx=64) still get NAN in multigrid L2

Also experimented with ZG on no-slip for mg eta bcs (as above), but with fixing one of the phi values.

• This avoids the crash (if either top or bottom phi is fixed, arbitrarily to zero) previously noticed at nx=64
• Results look good otherwise.
• Conditioning looks decent, V-cycles are still <=O(20)

Leads to questions

• does it work for higher nx / longer run times (i.e. is this a robust fix or a fluke)
• does it fix all flavours

[JOThurgood]

Tried the "fixed wrong" conditions (eta fixed, but set to =7 on top instead of expected 1/7, which is also likely wrong) for a higher resolution

nx = 64 -> 128

for t_end = 10

It completes for 64 (already known) but for 128 it fails at about t=3.66 (NAN in mg converge).

This suggests that eta fixed=7 = rho at the top wasn't magically making it robust against the crash, it just so happened to remain solvable at that resolution. There is nothing special about the value.

[JOThurgood]

does it work for higher nx / longer run times (i.e. is this a robust fix or a fluke)

set eta=zg and phi = fixed = 0 on bottom in mg arguments

nx = 64 confirmed working earlier

tried nx = 128, found

• It crashes a bit earlier than its "fixed to 7" equivalent (t=3.3 rather than 3.6)
• conditioning is a bit better

This is very disappointing, I thought it had proven to be a fix.

The same, but with fixing phi on the top, runs a bit longer but also ultimately fails.

[JOThurgood]

Maybe if you fix phi on a top for a mg BC you also want it fixed in the fluid bc ?

[JOThurgood]

See issue #101 - might be that the ymin = no slip + ymax = outflow helps

Might also need to change bottom from no slip to slip.

[JOThurgood]

Also - check if it is because nx/= ny in MG. Might be that the bottom solve isn't going anywhere (this is the first anisotropic case you've considered).

Things like setting up the problem with nx=ny (i.e. just simulate more periodic domains st Lx=Ly nx=ny), and also diagnostics from MG's bottom solve.

[JOThurgood]

see if .. ymin = no slip + ymax = outflow helps

Still get crashes.

Asking the MG to print out L2 on nonconvergence gives L2=NaN.
Asking to print out phi gives phi = NaN (everywhere, apparently).
Asking to print out eta gives finite array.

[JOThurgood]

Noticed for the setup where x = -2, 2 but "d" in initial conditions isn't /4 that you get some divergence at the lower boundary early on. Could be a hint about improper treatment down there.

[JOThurgood]

In addition to slip conditions, should see if linear extrapolations of density (as well as pressure) work

[JOThurgood]

Attempts to replicate the bug in a less expensive setup

If you set the same problem for a shorter tube

  subroutine rti1_control                                            
    x_min = -0.5_num
    x_max = 0.5_num
    y_min = -1.0_num
    y_max = 1.0_num
    nx = 128 
    ny = 2 * nx
    CFL = 0.5_num
    t_end = 10.0_num 
    nsteps = -1
    use_minmod = .false.    
    use_viscosity = .false.
    visc = 0.0_num
    bc_xmin = periodic
    bc_xmax = periodic
    bc_ymin = no_slip
    bc_ymax = outflow
    dumpfreq = -1  
    dt_snapshot = 0.10_num
    use_vardens = .true.
    grav_x = 0.0_num
    grav_y = -1.0_num
  end subroutine rti1_control

You'll get a NaN crash eventually. Lower res (nx=64) doesn't crash. Also the crash isn't when the droplet hits the bottom but rather the top - so it may or may not be the same issue.

An even shorter also crashes:

    x_min = -0.5_num
    x_max = 0.5_num
    y_min = -0.5_num
    y_max = 0.5_num
    nx = 128 
    ny = nx
    CFL = 0.5_num
    t_end = 10.0_num 
    nsteps = -1
    use_minmod = .false.        
    use_viscosity = .false.
    visc = 0.0_num
    bc_xmin = periodic
    bc_xmax = periodic
    bc_ymin = no_slip
    bc_ymax = outflow
    dumpfreq = -1  
    dt_snapshot = 0.10_num
    use_vardens = .true.
    grav_x = 0.0_num
    grav_y = -1.0_num    

edit - actually the above crashes for nx=64 too (but not 32)

[JOThurgood]

For the smallest above case, find that there is a vacuum formation in spots and eventually -ve density (that aren't located at boundaries but rather interior). The crash occurs the first time -ve density is produced in advect_dens and correspondingly passed to mg via eta.

The vacuum formation is avoided in coarser simulations (which do not crash).

Confirmed this behaviour is the case in all the RTI crashes.

[JOThurgood]

did quick test - varying cfl modifier in 0.1-0.9 range doesn't stop crashes.

[JOThurgood]

Gradient limiting might circumvent this behaviour?

[JOThurgood]

Gradient limiting might circumvent this behaviour?

it does help.

Close this issue for now.