enforcing divergence free condition

[Makrand Ajay Khanwale]

Hello developers and users,

I am running into an interesting issue.  I am trying to use high-pass filter velocity in the linear forcing to Navier stokes like so A_ij (u_j), where (u_j) is the high-pass filter velocity which I compute in the main solver loop. However, after the filtering the velocity is not divergence free and it causes unphysical results (assymetry in Reynolds stresses). 

I am trying to use a LBVP solver to project the velocity in the divergence free space. The way I do it is as follows.

```
# Boundary tau field (for Chebyshev z-basis, scalar LBVP)
tau_phi = dist.Field(name='tau_phi')

# Two tau fields for z-boundary Neumann BCs
# tau_lower = dist.Field(name='tau_lower')  # for z=0
# tau_upper = dist.Field(name='tau_upper')  # for z=Lz

# Scalar potential field (Chebyshev z, Fourier x/y)
phi = dist.Field(name='phi', bases=(xbasis, ybasis, zbasis))
phi_rhs = dist.Field(name='phi_rhs', bases=(xbasis, ybasis, zbasis))

#laplacian_phi = d3.lap(phi) + lift(tau_phi)

# Build LBVP with tau lifting (no separate Neumann BC fields)
poisson = d3.LBVP([phi, tau_phi], namespace=locals())

# Laplacian with tau for BCs
poisson.add_equation("lap(phi) + tau_phi = phi_rhs")  
# Neumann BCs on z-walls
#poisson.add_equation("dz(phi)(z=0) = 0")
#poisson.add_equation("dz(phi)(z=Lz) = 0")

# Gauge condition for nullspace (fix phi constant)
poisson.add_equation("integ(phi) = 0")  

poisson_solver = poisson.build_solver()

def project_u_hp_cheb_z(u_hp_field, dealias=dealias):
   
    # 1) compute divergence field (Dedalus Field)
    div_u = d3.div(u_hp_field)        

    # 2) build mean (Dedalus Field scalar) and subtract (still a Field)
    #    mean_div_field is a Dedalus Field representing the domain average of div_u
    mean_div_field = d3.integ(div_u) / Volume   # Dedalus Field (scalar-like)

    # 3) fill phi_rhs safely by evaluating the full field (div - mean)
    #    evaluate the Field expression on the distributed grid (no hangs)
    phi_rhs.change_scales(dealias)              # ensure correct grid
    phi_rhs['g'][:] = (div_u - mean_div_field).evaluate()['g']

    # 4) solve LBVP (lap(phi) + tau_phi = phi_rhs). Tau enforces Neumann BCs.
    poisson_solver.solve()

    # 5) evaluate gradients (include lift of tau in z)
    #    change_scales for phi/tau to avoid aliasing / shape mismatch
    phi.change_scales(dealias)
    # If you made tau_phi with no bases, change_scales is harmless
    try:
        tau_phi.change_scales(dealias)
    except Exception:
        pass

    # Evaluate dx/dy/dz on the physical grid and include lift(tau_phi) in dz.
    dx_phi = dx(phi).evaluate()['g']
    dy_phi = dy(phi).evaluate()['g']
    # IMPORTANT: include lifted tau contribution in z-derivative
    dz_phi = dz(phi).evaluate()['g'] + lift(tau_phi).evaluate()['g']

    # 6) subtract gradient from u_hp_field (in-place)
    u_hp_field['g'][0] -= dx_phi
    u_hp_field['g'][1] -= dy_phi
    u_hp_field['g'][2] -= dz_phi

    # 7) diagnostics: compute RMS divergence using Dedalus integral (safe)
    div_after_field = d3.div(u_hp_field)                    # Field
    sumsq = d3.integ(div_after_field * div_after_field)     # Field scalar
    # Evaluate once and extract scalar
    sumsq_val = float(sumsq.evaluate()['g'].flatten()[0])
    rms_div = np.sqrt(sumsq_val / Volume)
    logger.info("proj check: div(u_hp) RMS=%g (min=%g max=%g)",
                rms_div,
                float(np.min(div_after_field.evaluate()['g'])),
                float(np.max(div_after_field.evaluate()['g'])))
```

However, I get factor singular errors with this.  I am attaching the full script if anyone wants to run it as a text file. 
I may be making a mistake in how the tau method is used and how the gauge constraints are enforced.  The boundary condition should be zero Neumann on the wall (z direction Chebyshev). 

Additionally, if there is a better way to project a velocity field onto a divergence-free space that can be used directly in D3, I would be happy to use that instead as well.
Any help and insights are much appreciated.

With sincere regards,

Makrand

[Keaton Burns]

Hi Makrand,

This is a bit hard to follow because of everything commented out, but yes its generally true that spectral filtering with Chebyshev does not retain the solenoidal property, and yes you should be able to do the filtering and then a divergence cleaning step. It would probably help if you wrote just a clean, isolated divergence cleaning LBVP and started with debugging that in isolation, and then combining it with the IVP.

Best,

-Keaton

--
You received this message because you are subscribed to the Google Groups "Dedalus Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dedalus-user...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/dedalus-users/38587cff-a89e-4925-adbb-414e4e87aa39n%40googlegroups.com.

[Makrand Ajay Khanwale]

Hello Keaton,

I created a minimal example of LBVP, which can be run in isolation. See the attached file (attaching as a .txt file). It still produces large values of divergence, and I am wondering if I messed up some tau setup or the setup of the LBVP.  The precise combination of taus and BCs still eludes my understanding in this scenario.

I would really appreciate your input on this setup.

With sincere regards,

Makrand

[Makrand Ajay Khanwale]

Hello everyone,

I have been trying multiple variations of the tau lifting and boundary conditions, the LBVP solver above still returns a field that has a large divergence.  Any insights and help on this is much appreciated.

With sincere regards,

Makrand