Several questions on defining (anisotropic, non-linear, ...) Coefficients

[cgadal]

Following #1086, I am solving:

$$\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$.

Using fiPy 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)

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, :])

diffcoeff = FaceVariable(name="D", mesh=mesh, value=D)
segcoeff = CellVariable(name="F", mesh=mesh, rank=1)
segcoeff[0, :] = 0
segcoeff[1, :] = F

# ###### define equation
eq = TransientTerm(var=phi) + VanLeerConvectionTerm(coeff=segcoeff * (1 - phi), var=phi) + VanLeerConvectionTerm(
    coeff=velocity, var=phi
) == DiffusionTerm(coeff=diffcoeff, 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

This leads to excellent results.

However, defining the entire first ConvectionTerm coefficient outside as:

segcoeff = CellVariable(name="F", mesh=mesh, rank=1)
segcoeff[0, :] = 0
segcoeff[1, :] = F* (1 - phi)  # <------------- Here

# ###### define equation
eq = TransientTerm(var=phi) + VanLeerConvectionTerm(coeff=segcoeff , var=phi) + VanLeerConvectionTerm(
    coeff=velocity, var=phi
) == DiffusionTerm(coeff=diffcoeff, var=phi)

induces completely different (and wrong) results.

1. Why is that?

2. I am moving towards complex non-linear expressions of the diffusion and convection coefficients that will be hard to define directly inside the creation of each term. How can I create a proper diffusion coefficient as a CellVariable or FaceVariable involving the dependent variables?

EDIT: looking through the examples talking about anisotropy, found that I could do something like:

seg_prefact = (u_z * d_av / (pi_d * d_av + p)) * (phi_m - phi) * (phi_m - lamb * phi)
# seg_prefact = phi_m - phi
seg_vec = Variable(value=[0, -1])
seg_coeff = seg_prefact * seg_vec

Is there a difference between using Variable and CellVariable here?

3. The documentation mentions that Higher order diffusion expressions can be expressed as DiffusionTerm(coeff=(Gamma1, Gamma2)). However, it also says in the FAQ that DiffusionTerm([5, 1]) represents anisotropic diffusion with 0 diagonal terms. Is the difference between the use of a tuple and a list?

[guyer]

This leads to excellent results.

Yay!

However, defining the entire first ConvectionTerm coefficient outside as:

segcoeff = CellVariable(name="F", mesh=mesh, rank=1)
segcoeff[0, :] = 0
segcoeff[1, :] = F* (1 - phi)  # <------------- Here

# ###### define equation
eq = TransientTerm(var=phi) + VanLeerConvectionTerm(coeff=segcoeff , var=phi) + VanLeerConvectionTerm(
    coeff=velocity, var=phi
) == DiffusionTerm(coeff=diffcoeff, var=phi)

induces completely different (and wrong) results.

1. Why is that?

When you write

segcoeff[1, :] = F* (1 - phi)

you are assigning the instantaneous value of that expression to segcoeff[1, :]; it won't update as phi changes.

2. I am moving towards complex non-linear expressions of the diffusion and convection coefficients that will be hard to define directly inside the creation of each term. How can I create a proper diffusion coefficient as a CellVariable or FaceVariable involving the dependent variables?

EDIT: looking through the examples talking about anisotropy, found that I could do something like:

seg_prefact = (u_z * d_av / (pi_d * d_av + p)) * (phi_m - phi) * (phi_m - lamb * phi)
# seg_prefact = phi_m - phi
seg_vec = Variable(value=[0, -1])
seg_coeff = seg_prefact * seg_vec

Exactly.

Assembling coefficients is illustrated in:

• examples.phase.anisotropy
• examples.phase.polyxtal

Is there a difference between using Variable and CellVariable here?

For fixed orientations, a Variable isn't necessary, as illustrated in examples.cahnHilliard.mesh2DCoupled

Variable and CellVariable work together. Variable is for single or multiple values that are not located in space. CellVariable is for a field of values evaluated at the center of the cells of the spatial mesh of your domain.

Multiplying a scalar CellVariable field by a vector Variable value results in a vector CellVariable field.

3. The documentation mentions that Higher order diffusion expressions can be expressed as DiffusionTerm(coeff=(Gamma1, Gamma2)). However, it also says in the FAQ that DiffusionTerm([5, 1]) represents anisotropic diffusion with 0 diagonal terms. Is the difference between the use of a tuple and a list?

The bit in the FAQ could perhaps be written more clearly.

DiffusionTerm([[[5, 0], [0, 1]]])

represents a 2nd order anisotropic diffusion term for a 2D field:

$$\nabla\cdot\left[\begin{pmatrix} 5 & 0 \\\ 0 & 1 \end{pmatrix}\nabla\phi\right]$$

DiffusionTerm([[5, 1]])

(a list containing a list containing two values) represents the same 2nd order anisotropic diffusion term.

DiffusionTerm([5, 1])

(a list containing two values) represents the 4th order isotropic diffusion term:

$$\nabla\cdot\left[5 \nabla\cdot\left(1 \nabla\phi\right)\right]$$

Although supported, we no longer recommend the use of higher order diffusion expressions. The same result can be obtained by factoring higher order PDEs into a coupling of 2nd order PDEs.

[guyer]

Great questions, BTW