Solving a PDE problem considering the nonlinear diffusion coefficient of Richards equation

[essicolo]

Hello Fipy's community,

I'm working on a PDE problem, which was recently reformulated by Dr. Guyer on Stackoverflow. The problem is represented as:

$$ \frac{\partial \theta}{\partial t} + \nabla \cdot ( k(u) \nabla (u+z) ) - S = 0 $$

↳ I've edited this equation after @guyer's comment.

$$ \frac{\partial \theta}{\partial t} + \nabla \cdot ( k(u) \nabla u ) + \nabla \cdot ( k(u) \nabla z ) - S = 0 $$

$$ \frac{\partial \theta}{\partial t} + \nabla \cdot ( k(u) \nabla u ) + \nabla \cdot ( k(u) \hat{k} ) - S = 0 $$

$$ \frac{\partial \theta(u)}{\partial u} \times \frac{\partial u}{\partial t} + \nabla \cdot ( k(u) \nabla u ) + \nabla \cdot ( k(u) \hat{k} ) - S = 0 $$

Without the sink/source term and in steady-stade, we have:

$$ \nabla \cdot ( k(u) \nabla u ) + \nabla \cdot ( k(u) \hat{k} ) = 0 $$

Where:

• $\hat{k}$ is a unit vector in the $z$ direction (elevation).
• $\theta(u)$ is a non-linear function that is ruled out in steady state
• $k(u)$ is a non-linear function defined as:

$$ k(u) = k_{sat} \times (1 - (-a ~ u)^{n-1} \times (1 + (-a ~ u)^n)^{-m})^2 \times (1 + (-a ~ u)^n)^{-m ~ l} $$

With $k_{sat}$, $a$, $n$, $m$, and $l$ being material parameters.

Issue. After implementing Dr. Guyer's suggestions, my solution doesn't align with two analytical benchmarks. The expected pressure head should start at zero (the Dirichlet BC), then decrease and stabilize to a minimum value at a certain elevation from the water table. However, the FiPY solution appears influenced by the column height (Lz). The initial guess also heavily impacts the results. For instance, with an initial guess of 0.0, all u values (pressure head) remain at zero. Different solvers yield similar outcomes.

Request. I'm seeking suggestions to adjust my code in a way it fits approximately to my benchmarks. You can review it here.

Thank you for your assistance!

Best regards,

Essi

[guyer]

my solution doesn't align with two analytical benchmarks. The expected pressure head should start at zero (the Dirichlet BC), then decrease and stabilize to a minimum value at a certain elevation from the water table. However, the FiPY solution appears influenced by the column height (Lz).

Why wouldn't it be sensitive to the column height? What in the math would cause it to stabilize at some arbitrary elevation? I'm afraid $k(u)$ is too complicated for me to have any intuition about its value or about what effect it might have. I see that $k(u)$ has an inflection at intermediate depth, but I still don't know why that should mean that $u$ should saturate.

For instance, with an initial guess of 0.0, all u values (pressure head) remain at zero. Different solvers yield similar outcomes.

$u = 0$ is a solution to the steady state equation.

[guyer]

Better, thank you. Where is $q$ in your PDE? You use it to calculate the initial condition and then never again as far as I can tell.

Presently, at the right boundary, you have FiPy's natural boundary condition of zero flux (this is the natural boundary condition of the finite volume method). Note that, for your problem, flux is $-k(u)\nabla(u + z)$. Zero flux does not mean that $\hat{n}\cdot\nabla u$ is equal zero.

note: I think you have an error in your first equation at top. I think it should be

$$ \frac{\partial \theta}{\partial t} + \nabla \cdot ( k(u) \nabla (u+z) ) - S = 0 $$

[essicolo]

In fact, maybe the whole problem should be expressed in terms of total head $h$ ($h=u+z$)? This will make the code easier to express since, as you wrote, the flux is related to the gradient of $h$, not $u$ — it took me so much time to figure that out...

Action. I rearranged the code to solve for h instead. I'm wondering if, once expressed on $h$, is the + (k.faceValue * [[1]]).divergence term is still needed?

Results. It seems that eq.sweep(var=h, solver=solver) didn't alter the initial guess and all $h$ values are zeros.

[guyer]

Yes, solving for $h$ is probably better, in which case, no, you don't need (k.faceValue * [[1]]).divergence.

I'm not getting $h = 0$ from revision 13.

Please express your boundary conditions mathematically. What exactly is $q$?

[essicolo]

EDITED: You are right, I'm not getting $h=0$ either. Moreover, there was a sign error, in if np.max(np.abs(h.value - h.old.value)) > tolerance, which obviously should have been if np.max(np.abs(h.value - h.old.value)) < tolerance. Now h values are Inf.

What exactly is q?

I'll do my best, mathematically explained. $q$ in the context of my code is a flux going inward the domain, such as

$$ q = - \hat{n} \cdot k(u) \nabla h $$

where $u$ is a cell variable such as $u = h-z$, and $z$ is the dimension of the domain representing elevation.

The Neumann boundary condition expressed as $\nabla h$ is thus

$$ \nabla h = - \frac{q}{\hat{n} \cdot k(u) } $$

which I would express in my FiPY problem as h.faceGrad.constrain(- q/k.faceValue, where=mesh.facesRight).

Thanks again for your time!

[guyer]

OK, I assumed $q$ was something like that, but wanted to be sure. Just to confirm, it's definitely $q=-\hat{n}\cdot k(u)\nabla h$ and not $q=-\hat{n}\cdot k(u)\nabla u$?

Regardless, I've tried a couple of different formulations (it's actually more natural for finite volume methods to constrain the flux than the gradient[*]), but all cycle between a few (wrong) solutions. This behavior only reinforces two prejudices I've developed over the years: 1D formulations are always wrong and steady-state solutions are nonsense.

If I institute relaxation (artificial time evolution), I am able to obtain good agreement with your Kisch solution.
.

Specifically,

diff --git a/richards.py b/richards.py
index ee8ef7a..669cdb3 100644
--- a/richards.py
+++ b/richards.py
@@ -2,6 +2,8 @@ import numpy as np
 from scipy.optimize import fsolve
 import matplotlib.pyplot as plt
 import fipy as fp
+from petsc4py import PETSc
+
 PETSc.Options().setValue("snes_type", "newtonls")
 
 def vg_k(u, k_sat, a, n, m, l=0.5):
@@ -33,24 +35,23 @@ def richards_fem(q, k_sat, a, n, m, l=0.5, Lz=10, nz=100, max_iters=100, toleran
     # Variables
     z = fp.CellVariable(name="elevation", mesh=mesh, value=mesh.cellCenters[0])
     h = fp.CellVariable(name="total head", mesh=mesh, value=0.0, hasOld=True)
-    u = fp.CellVariable(name="matrix pressure", mesh=mesh, value=h.value-z.value, hasOld=True)    
-    k = fp.CellVariable(mesh=mesh, value=vg_k(u.value, k_sat, a, n, m, l))
+    u = h - z
+    u.name = "matrix pressure"
+    k = vg_k(u.faceValue, k_sat, a, n, m, l)
 
     # Boundary conditions
     h.constrain(0.0, mesh.facesLeft) # Dirichlet at bottom
-    h.faceGrad.constrain(- q/k.faceValue, where=mesh.facesRight)  # Neumann at top
+    h.faceGrad.constrain(q/k, where=mesh.facesRight)  # Neumann at top
 
     # Equation
-    eq = (fp.DiffusionTerm(coeff=k, var=h) == 0)
+    eq = (fp.TransientTerm(var=h) == fp.DiffusionTerm(coeff=k, var=h))
 
     # Solve with Picard iterations
     solver = fp.DefaultSolver(tolerance=1e-10, iterations=5000)
     for iter in range(max_iters):
-        u.value = h.value - z.value
-        k.value = vg_k(u.value, k_sat, a, n, m, l)
-        eq.sweep(var=h, solver=solver)
-        if np.max(np.abs(h.value - h.old.value)) > tolerance:
-            break
+        h.updateOld()
+        for sweep in range(3):
+            res = eq.sweep(dt=10000., var=h, solver=solver)
 
     return u.value, z.value
 
@@ -73,11 +74,12 @@ z_kisch = kisch(q=rain, u=u_kisch, u_q=u_q, k_sat=k_sat, a=a, n=n, m=m, l=l, z_b
 
 u_fem, z_fem = richards_fem(
     q=rain, k_sat=k_sat, a=a, n=n, m=m, l=l,
-    Lz=10, nz=100, max_iters=1000, tolerance=1e-8
+    Lz=10, nz=1000, max_iters=10000, tolerance=1e-8
 )

Notes:

• I have rewritten u and k to take advantage of FiPy's lazy evaluation. There is no need to manually set them after every solution.
• I had to flip the sign of your Neumann boundary condition.
• Time steps above about 10000 were unstable.
• Your if np.max(np.abs(h.value - h.old.value)) > tolerance: is backwards, but also isn't really helpful with relaxation and sweeping.

[*] To constrain flux, rather than gradient, make these additional changes:

diff --git a/richards.py b/richards.py
index 669cdb3..ab1d307 100644
--- a/richards.py
+++ b/richards.py
@@ -41,10 +41,11 @@ def richards_fem(q, k_sat, a, n, m, l=0.5, Lz=10, nz=100, max_iters=100, toleran
 
     # Boundary conditions
     h.constrain(0.0, mesh.facesLeft) # Dirichlet at bottom
-    h.faceGrad.constrain(q/k, where=mesh.facesRight)  # Neumann at top
 
     # Equation
-    eq = (fp.TransientTerm(var=h) == fp.DiffusionTerm(coeff=k, var=h))
+    eq = (fp.TransientTerm(var=h)
+          == fp.DiffusionTerm(coeff=k, var=h)
+          + (q * mesh.faceNormals * mesh.facesRight).divergence)
 
     # Solve with Picard iterations
     solver = fp.DefaultSolver(tolerance=1e-10, iterations=5000)

[essicolo]

This look promising! However, I couldn't test the results since the patch returned an index error

IndexError: diffusion coefficient tensor is not an appropriate shape for this mesh

Have you changed the mesh to higher dimensions?

[guyer]

I made only the changes I showed. Update your gist and I'll take a look.

[essicolo]

Ok thanks. Gist updated!

[guyer]

You can't use np.where that way. NumPy doesn't know anything about FiPy's lazy mesh variables.

You can do this:

diff --git a/richards.py b/richards.py
index 0ea2b1c..d0beab1 100644
--- a/richards.py
+++ b/richards.py
@@ -7,11 +7,8 @@ from petsc4py import PETSc
 PETSc.Options().setValue("snes_type", "newtonls")
 
 def vg_k(u, k_sat, a, n, m, l=0.5):
-    k = np.where(
-        u>=0,
-        k_sat,
-        k_sat * (1 - (-a * u)**(n-1) * (1 + (-a * u)**n)**(-m))**2 * (1 + (-a * u)**n)**(-m*l)
-    )
+    k = k_sat * ((u >= 0) * 1.
+                 + (u < 0) * (1 - (-a * u)**(n-1) * (1 + (-a * u)**n)**(-m))**2 * (1 + (-a * u)**n)**(-m*l))
     return k
 
 def find_u_for_given_k(target_k, k_sat, a, n, m, l=0.5, u_guess=-1.0):

[essicolo]

Hello Dr Gruyer,

I'm pleased to inform you that the simulation was a success! I've made some modifications to the code, which you can view in this gist. Remarkably, after just a few iterations, the results align closely with the Kisch profile.

In the coming days, I plan to further enhance the code by exploring higher dimensions, incorporating different materials, and introducing more intricate boundary conditions. The foundation you've provided has opened up numerous possibilities, and I'm genuinely enthusiastic about the future potential of this project.

Thank you once again for your invaluable assistance.

--
Essi