FiPy or Clawpack for advection-diffusion like equation in non-rectangular domain

[cgadal]

Dear all,

I am trying to solve the following equation:

$$\frac{\partial \phi}{\partial t} + \nabla\cdot(\phi \boldsymbol{u}) - \frac{\partial}{\partial z} \left[\mathcal{F} \phi(1-\phi)\right] = \nabla \cdot (\mathcal{D} \nabla\phi),$$

where:

• $\phi$ is a scalar, $\boldsymbol{u}$ is a vector field advecting the scalar
• $\mathcal{D}$ and $\mathcal{F}$ are coefficients that can be taken constant, but that will probably become functions of $\phi$ but also $x$ and $z$.

The domain is bounded by two parabolas (red and brown) :

and the black lines represent streamlines of the advection field $\boldsymbol{u}$.

The boundary conditions are non-flux through the boundaries. As the vector field also respects this condition, it becomes:

$$ ( \mathcal{D} \nabla\phi ) \cdot \boldsymbol{n} + \mathcal{F} \phi(1-\phi) n_{z} = 0,$$

where $\boldsymbol{n} = (n_{x}, n_{z})$ is a vector normal to the boundary.

From what I have read in FiPy documentation:

• The first, second and RHS terms in the equation can be discredited using, in the right order: TransientTerm, PowerLawConvectionTerm, DiffusionTerm.

• the non-rectangular domain can be dealt with by generating a mesh using Gmsh

Remaining questions:

• I am unsure about how to discretize the third term of the equation. It does not seem to fall into any of the typical categories, so should I deal with it as a source term ?
• I am unsure about how to represent those boundary conditions, but it seems to me they match the default bc used by Fipy?
• If I set the diffusion coefficient $\mathcal{D}$ to 0, the system becomes hyperbolic. What can I expect from FiPy?
• Should I expect issues when specifying more complex coefficients $ \mathcal{D}$ and $\mathcal{F}$ (i.e function of space, or $\phi$) ?

More generally, would ClawPack be more suited for this? I initially chose FiPy because the documentation is much clearer, and also because I found no example of using custom meshes for non-rectangular domain in ClawPack documentation and examples.

[guyer]

Remaining questions:

• I am unsure about how to discretize the third term of the equation. It does not seem to fall into any of the typical categories, so should I deal with it as a source term ?

I would be inclined to approach the third term by writing

$$ -\frac{\partial}{\partial x}\left[\hat{\imath} F_x \phi \left(1 - \phi\right)\right] -\frac{\partial}{\partial z}\left[\hat{k} F_z \phi \left(1 - \phi\right)\right] = -\nabla\cdot\left[\vec{F} \phi \left(1 - \phi\right)\right] $$

and then representing that as

F = ...
Fvec = [0, F]
eq = ... + ConvectionTerm(coeff=Fvec * (1 - phi), var=phi) == ...

• I am unsure about how to represent those boundary conditions, but it seems to me they match the default bc used by Fipy?

Correct. The default BC for FiPy are no flux.

• If I set the diffusion coefficient

  D
  

  to 0, the system becomes hyperbolic. What can I expect from FiPy?

There are better tools for purely hyperbolic equations (such as Clawpack), but I have successfully used FiPy to solve Beer's law. I think you should be fine.

• Should I expect issues when specifying more complex coefficients $ \mathcal{D}$ and

  F
  

  (i.e function of space, or
  ϕ
  ) ?

More complex coefficients should be fine. You'll probably need to sweep each time step to convergence.

More generally, would ClawPack be more suited for this? I initially chose FiPy because the documentation is much clearer, and also because I found no example of using custom meshes for non-rectangular domain in ClawPack documentation and examples.

I have not personally used Clawpack, but I think FiPy should be able to handle what you're trying to do. If you're going to be in an extreme range of fluid dynamics, such as involving shocks, Clawpack is probably better. Maybe @wd15 can weigh in.

[cgadal]

Thank you for this (amazingly) quick answer! I will try this and tell you when I have succeeded or failed. I will start with reasonably high diffusion coefficients compared to convection to avoid shocks for starting.

[cgadal]

Script

After about a day of fiddling around, I wrote a code that gave promising first results.
The equation is defined as:

# ###### load mesh
mesh = Gmsh2D(os.path.join(f"{name}.geo_unrolled"), communicator=serialComm)
yz_meshface = np.array(mesh.faceCenters)
yz_meshcell = np.array(mesh.cellCenters)

# ###### define variables
phi = CellVariable(name=r"$\phi$", mesh=mesh, value=0.8, hasOld=sweep)
Fvec = [0, F]

velocity = FaceVariable(name="vel", mesh=mesh, rank=1)
yz_meshface = np.array(mesh.faceCenters)

velocity[0, :] = U * MODEL_DFLT.compute_v(yz_meshface[0, :], yz_meshface[1, :])
velocity[1, :] = U * MODEL_DFLT.compute_w(yz_meshface[0, :], yz_meshface[1, :])

# ###### define equation
eq = TransientTerm(var=phi) + VanLeerConvectionTerm(coeff=Fvec * (1 - phi), var=phi) + VanLeerConvectionTerm(
    coeff=velocity, var=phi
) == DiffusionTerm(coeff=D, var=phi)

and the temporal loop:

while res >= stop_condition:
    print(step)
    if sweep:
        phi.updateOld()
        for sweep in range(sweep_number):
            res_sweep = eq.sweep(var=phi, dt=timestep)
            print(f"res_sweep = {res_sweep:1.2e}")
    else:
        eq.solve(var=phi, dt=timestep)
    res = np.abs(phi.value - phi_old).max()
    print(f"res_dt = {res:1.2e}")
    phi_old = np.copy(phi.value)
    step += 1

Results

I obtain convergence for values of the diffusion coefficient D (here D = 0.05, F = -1, U = 3):

However, when D becomes too low, 2 unstable behaviours are observed:

1. if D becomes too small compared to F, waves are generated at the bottom boundary and propagate towards the top while growing.
2. if D becomes too small compared to U, some mass builds up at the surface on the left and right sides.

Both effects can be seen in this image:

Note that both are not physical as $\phi$ gets larger than 1.

What I have tried

• sweeps or smaller time steps: Strangely, 1. can be helped by using sweeps (and to some extent by decreasing the time step), allowing smaller diffusion coefficients (but it always ends up unstable) to be used, but this does not affect 2.

• Convection schemes: I have used both PowerLawConvectionScheme and VanLeerConvectionTerm, hoping that the higher order of the second would help, but I have seen no visible effect switching between the two.

Would you have, by any chance, any ideas about all of this?

[guyer]

I'm happy to see you've made as much progress as you have. My suggestions would be to tinker with sweeping and to try the VanLeerConvectionTerm, but you've already done both.

What happens with res_sweep? Does it appear to be converged?

You are defining velocity on mesh.faceCenters, which I think is the right answer, but you might try defining it on .cellCenters. The ConvectionTerm will interpolate it to faces. FWIW, I don't think this is better, but it's something to try. The displacement of velocity and pressure fields between faces and cells can be a source of problems in Navier-Stokes solutions. That's not explicitly what you have here, but maybe it's related.

@wd15 has much more experience with flow problems and may have other ideas.

[cgadal]

(just confirming that 1. can be entirely solved by adjusting the timesteps wrt the characteristics of the strictly hyperbolic system)

About 2.:

What happens with res_sweep? Does it appear to be converged?

Yes

You are defining velocity on mesh.faceCenters, which I think is the right answer, but you might try defining it on .cellCenters.

This neither seems to help nor to make it worse.

I am now more and more convinced that this is due to the discretization of the velocity field close to the boundaries, which is supposed to be divergence-free (and the analytical expression is), but its numerical estimation close to the boundaries might not be, resulting from curved boundaries not being perfectly mapped by the triangular mesh, hence cutting streamlines.

What I would like to try:

• maybe solving directly $\boldsymbol{u}\nabla\phi$ instead of $\nabla \cdot (\boldsymbol{u}\phi)$ ? I see in the API there are things called AdvectionScheme.

• Not enforcing a 0 flux boundary condition on the boundaries, but rather only $( \mathcal{D} \nabla\phi ) \cdot \boldsymbol{n} + \mathcal{F} \phi(1-\phi) n_{z} = 0$ and see if there are any differences? If this is due to the velocity problem mentioned above, doing this should not build up mass anymore but slowly leak some mass across the boundary.

However, while I have understood how to set both fluxes to 0 from the documentation, I have not yet understood how to set the sum of the two equal to 0.

[guyer]

(just confirming that 1. can be entirely solved by adjusting the timesteps wrt the characteristics of the strictly hyperbolic system)

👍

About 2.:

What happens with res_sweep? Does it appear to be converged?

Yes

Good

You are defining velocity on mesh.faceCenters, which I think is the right answer, but you might try defining it on .cellCenters.

This neither seems to help nor to make it worse.

I'm not surprised, but thanks for checking.

I am now more and more convinced that this is due to the discretization of the velocity field close to the boundaries, which is supposed to be divergence-free (and the analytical expression is), but its numerical estimation close to the boundaries might not be, resulting from curved boundaries not being perfectly mapped by the triangular mesh, hence cutting streamlines.

That sounds plausible to me. Possibly related, I think FiPy doesn't handle (or doesn't handle well) flux components tangent to the boundaries. At one point, Daniel was working on a Riemann solver that I think was supposed to help with this, but that never went anywhere.

What I would like to try:

• maybe solving directly
  ...
  ? I see in the API there are things called AdvectionScheme.

I'm not sure AdvectionTerm has been (or can be) used for anything other than level-set problems. That's in @wd15's department.

• Not enforcing a 0 flux boundary condition on the boundaries, but rather only
  ...
  and see if there are any differences? If this is due to the velocity problem mentioned above, doing this should not build up mass anymore but slowly leak some mass across the boundary.

However, while I have understood how to set both fluxes to 0 from the documentation, I have not yet understood how to set the sum of the two equal to 0.

Take a look at the discussion on applying Robin conditions.

[cgadal]

I realized that there was a small but existent typo in a term of the velocity field, inducing some streamlines to cut the top boundary. Once fixed, everything behaved as expected.

Some small noise might arise at some places, but if the grid, time steps, and sweeps are chosen well, they do not grow and invade the system.

Thanks again for all your help!

[guyer]

Great news!