Hello FEniCS community,
I am working on a simulation involving the Richards equation for unsaturated flow in soils with fenics 2019.1.0 installed from conda-forge. The general form of Richards’ equation can be written as:
\frac{\partial \theta(u)}{\partial t} = \nabla \cdot [K(u) (\nabla u + \nabla z)]
where
-
\theta(u) is the volumetric water content as a function of pressure head u.
-
K(u) is the hydraulic conductivity as a function of pressure head u.
-
z is the elevation.
For variational formulation,
\int_{\Omega} \frac{\partial \theta(u)}{\partial t} v ~ d\Omega - \int_{\Omega} K(u) \nabla (u + z) \cdot \nabla v ~ d\Omega = \int_{\partial\Omega} v K(u) \nabla (u + z) \cdot \mathbf{n} ~ ds
Rearranging,
\int_{\Omega} \frac{\theta(u) - \theta(u_{old})}{\Delta t} v ~ d\Omega - \int_{\Omega} K(u) \nabla (u + z) \cdot \nabla v ~ d\Omega - \int_{\partial\Omega} v q ~ ds = 0
Here, u_{old} is the pressure head from the previous time step and \Delta t is the time step size.
My expectation is that the pressure head, u, would converge over time to a profile where u decreases linearly with elevation (x[1]), then breaks to a minimal value of u such as K(u) = q. My current results don’t match this expectation. Below is my complete code for the simulation:
import numpy as np
from fenics import *
import matplotlib.pyplot as plt
# Soil parameters
a = 0.2
theta_r = 0.05
theta_s = 0.45
n = 2.0
m = 1-1/n
k_sat = 1E-3
l = 0.5
# Soil water content and hydraulic conductivity
theta = lambda u: conditional(lt(u, 0), theta_r + (theta_s - theta_r) / (1 + (-a * u) ** n) ** m, theta_s)
K = lambda u: conditional(lt(u, 0), k_sat * (1 - (-a * u)**(n-1) * (1 + (-a * u)**n)**(-m)) * (1 + (-a * u)**n)**(-m*l), k_sat)
# Mesh and function space
ymin, ymax = 0, 50
xmin, xmax = 0, 1
nx, ny = 10, 50
mesh = RectangleMesh(Point(xmin, ymin), Point(xmax, ymax), nx, ny)
coordinates = mesh.coordinates()
V = FunctionSpace(mesh, 'CG', 1)
# Initial condition
#u_initial = interpolate(Expression("-x[1]", degree=1), V)
u_initial = Function(V)
u_initial.vector()[:] = 0.0
# Boundary conditions
## Dirichlet
class TopBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], ymax)
class BottomBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], ymin)
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
#boundaries.set_all(0)
TopBoundary().mark(boundaries, 1)
BottomBoundary().mark(boundaries, 2)
bc_bottom = DirichletBC(V, Constant(0.0), boundaries, 2)
## Neumann
q = Constant(1E-4) # as a source in the variational equation
# Variational problem
dt_ = 1000.0
T = 10000.0
## To store the previous solution and solve the next
u_old = Function(V) # function to hold the value from the previous iteration
u_old.assign(u_initial)
## To export the results
u_export = Function(V)
u_export.assign(u_initial)
## To solve the problem
u = Function(V)
v = TestFunction(V)
## Elevation
z_expression = Expression('x[1]',degree=1)
z = project(z_expression, V)
## Variational equation
F = (theta(u) - theta(u_old))/dt_ * v * dx + K(u) * dot(grad(u + z), grad(v)) * dx + q * v * ds(1)
# Solver
## Parameters
solver = NonlinearVariationalSolver(NonlinearVariationalProblem(F, u, bc_bottom))
prm = solver.parameters
prm['newton_solver']['linear_solver'] = 'gmres'
prm['newton_solver']['maximum_iterations'] = 500
prm['newton_solver']['absolute_tolerance'] = 1E-8
prm['newton_solver']['relative_tolerance'] = 1E-7
prm['newton_solver']['relaxation_parameter'] = 0.8
prm['newton_solver']['report'] = True
## Time-stepping and data extraction
u_all_timesteps = []
y_coords_unique = np.unique(mesh.coordinates()[:, 1])
y_indices = np.where(np.isclose(y_coords_unique, 0.5))[0]
## Solving loop
t = 0
while t <= T:
solve(F == 0, u, bc_bottom)
u_old.assign(u)
t += dt_
u_export.assign(u)
u_all_timesteps.append(u_export.copy())
# Plot data
x_center = (coordinates[:, 0].max() + coordinates[:, 0].min()) / 2.0
y_positions = np.unique(coordinates[:, 1])
u_at_center = []
for t in range(int(T/dt_)):
u_values_for_t = []
for y in y_positions:
values = np.array([0.0])
u_all_timesteps[t].eval(values, np.array([x_center, y]))
u_values_for_t.append(values[0])
u_at_center.append(u_values_for_t)
u_at_center = np.array(u_at_center)
for t in range(u_at_center.shape[0]):
plt.plot(u_at_center[t, :], y_positions)
I’ve played around with different solver parameters, and tried a Picard approach, without success. I suspect issues with the variational form. I would deeply appreciate any insights or suggestions you might have.
Thank you for your time and assistance, and best regards,
Essi
