CCAPS - No Slip + Grav issue

[jothurgood]

CCAPS - No Slip + Grav issue

[JOThurgood]

• If you set ZG on rho and v, you will avoid the issue of rho dropping below at top
• Of course, you can't set ZG as the bc on phi, otherwise you'll be messing up the correction
• E.g. if g = -1 and rho =1 , need grad_y p = -1. Zero gradient on phi is inconsistent with this.

[JOThurgood]

There is also a subtlety in that (I have the impression from literature) that often ZG on phi is a good / the best choice to pass to the multigrid solver. This could mean that you end up computing a phi around the edges that is polluted by the boundary somewhat (even if the interior L2 is good?) and then you'll not match rho g?

i.e. phi(:,iy=ny) itself is "wrong" when it comes out of MG and even say the backward difference calc of a gradient (i.e., not using any phi boundary data in the calculation)
(phi(nx/2,iy)-phi(nx/2,iy-1))/dy will not be what you want.

This also explains why the correction can be off at iy-1 - if you've used a centred difference then it has involved the bad phi(:,iy) data

Q - Why is there not an analogous issue at the bottom boundary ? There, you do get grad_p = rho g

[JOThurgood]

Q - Why is there not an analogous issue at the bottom boundary ? There, you do get grad_p = rho g

Maybe this hints that the MG is doing fine (despite using ZG on mg phi bc - which seems physically wrong - it could still converge to having the gradient in the domain ). Is it because only the top boundary exhibits the rho drop at the top cell, which in turn feeds through into MG via the 1/eta term, but the bottom remains rho = 1.

If in step 5 you pass eta = 1/1 instead of 1/rho to the mg, you get a much better gradient for phi at the iy = ny cells

[JOThurgood]

iy - altitude (ny=16)
grad_y psi - when tricked MG into thinking rho = 1
rho - actual rho in code at this stage

          iy                grad_y psi                        rho
           0 -0.99999999999996270        1.0000000000000000
           1 -0.99999999999996270        1.0000000000000000
           2 -0.99999999999995692        1.0000000000000000
           3 -0.99999999999993916        1.0000000000000000
           4 -0.99999999999992339        1.0000000000000000
           5 -0.99999999999990918        1.0000000000000000
           6 -0.99999999999989631        1.0000000000000000
           7 -0.99999999999988876        1.0000000000000000
           8 -0.99999999999988332        1.0000000000000000
           9 -0.99999999999988343        1.0000000000000000
          10 -0.99999999999988576        1.0000000000000000
          11 -0.99999999999989353        1.0000000000000000
          12 -0.99999999999990363        1.0000000000000000
          13 -0.99999999999991962        1.0000000000000000
          14 -0.99999999999993605        1.0000000000000000
          15 -0.99999999999995559        1.0000000000000000
          16 -0.99999999999997069       0.93750000000000000
          17 -0.99999999999997025        1.0000000000000000

This is much better.

Hopefully just moving to outflow (so that rho doesnt drop below 1) will fix it in of itself

[JOThurgood]

So when ymin = no_slip and ymax= zero gradient you get basically a divergence free flow except at the bottom boundary where, due to v = 0 on it, but v = - whatever nearby, you get a whopping divergence. This was there previously, but matched by an equal and opposite divergence at the top.

The MG converges to something bad with a huge L2 to try and eliminate this single source of divergence. It wont give you a phi-profile s.t. grad p = rho g

So ZG fixes the problem with rho "dripping" from the top, and so opens the door to a better correction there, but opens up this can of worms.

I wonder if, in the non-atmosphere subtracted codes, whether it would be worthwhile requiring the user to give grad p as part of the initial condition, in lieu of bootstrapping? Or, perhaps don't advect mass during the initial step ?

[JOThurgood]

when ymin = no_slip and ymax= zero gradient ... The MG converges to something bad with a huge L2 to try and eliminate this single source of divergence. It wont give you a phi-profile s.t. grad p = rho g

If, for the projection only, you set phi = 0 on the top boundary you'll get a linear phi profile out of MG. Then, if you also set phi = 0 in the hydro bcs, should get a reasonable correction at all cells, provided rho =1 remains everywhere (which ZG on velocity and rho otherwise should do?)

You still get the divergence now only at the bottom, but the MG manages to come up with a phi that has a linear gradient of -1 (of grad rho)

[JOThurgood]

Created a boundary type "outflow", which is currently only specified for ymax.

For the described test problem it correctly prescribes gradp to balance the atmosphere and then appears to remains stable. (tested for nx=16, ny=4nx, for t_end = 100 and div < 1e-8 by end - slowly growing vertical stripes - for higher resolutions nx=64, 128 magnitude of div is larger (but still small really) and everything is stable).

Number of V cycles on each step is good (2-5)

(hydro) Velocity

    ! outflow - zero gradient on v and rho, fixed phi

    if (bc_ymax == outflow) then

      if (present(arr_cc)) then
        arr_cc(:,ny+1) = arr_cc(:,ny)
        arr_cc(:,ny+2) = arr_cc(:,ny-1)
      endif

      if (present(arr_xface)) then
        arr_xface(:,ny+1) = arr_xface(:,ny)
        arr_xface(:,ny+2) = arr_xface(:,ny-1)
      endif

      if (present(arr_yface)) then
        arr_yface(:,ny+1) = arr_yface(:,ny-1)
        arr_yface(:,ny+2) = arr_yface(:,ny-2)
      endif

    endif

(hydro) Density

    if (bc_ymax == outflow) then

      if (present(arr_cc)) then
        arr_cc(:,ny+1) = arr_cc(:,ny)
        arr_cc(:,ny+2) = arr_cc(:,ny-1)
      endif

      if (present(arr_xface)) then
        arr_xface(:,ny+1) = arr_xface(:,ny)
        arr_xface(:,ny+2) = arr_xface(:,ny-1)
      endif

      if (present(arr_yface)) then
        arr_yface(:,ny+1) = arr_yface(:,ny-1)
        arr_yface(:,ny+2) = arr_yface(:,ny-2)
      endif

    endif

hydro Phi

    if (bc_ymax == outflow) then
      phi(:,ny+1) = phi(:,ny) + (phi(:,ny) - phi(:,ny-1))
      phi(:,ny+2) = phi(:,ny) + (phi(:,ny) - phi(:,ny-1))*2.0_num
    endif

multigrid phi and eta

    if (bc_ymax == outflow) then
      mg_bc_ymax = 'fixed'
      mg_etabc_ymax = 'zero_gradient'
    endif

Might be that there is a slight discrepancy between fixed eta in mg and then the extrapolation used for gradp calc in rest of code

[JOThurgood]

So in summary, this is basically because

• you want to enforce HSE in gravitational setups.

• Except in the background-subtracted model, this can only realistically be achieved by extrapolation.

• but usually at a solid wall rho and pressure would be chosen as zero gradient

I've refactored the boundaries and created no_slip slip no_slip_hse and slip_hse as distinct entities.

As for outflow, this is tricky. Usually on an outflow you'd want to enforce tangential pressure is zero (no shear on boundary), but this is somewhat at odds with achieving hse via linear extrap.

I've created outflow_hse, which prioritises hse over no shear (but treats velocity and density as zero grad). The idea is outflow could be implemented for non gravitational cases where the zero shear is respected. Additionally, for the atmosphere versions, this should be different.

Pushed to dev for CCAPS but not into the flavours, will be soon.