Parametarization of the constitutive relation for non-linear Poisson Equation and Netwon method does not converge

1

Hi, I am trying to solve time-independent 1D Richards equation (similar to a non-linear Poisson equation) to simulate water flow in porous media with the following constitutive relation.

VGM
VGM1025×339 24.3 KB

I have three issues in my code.

  1. I could not encode the gravitational gradient (\nabla z) in F(u; v). In order to make a vector, I used dl.Expression but it seems not working well. The error says “Cannot determine geometric dimension from expression.”.
  2. The newton method does not converge if I omit the gravitational term (it should converge though). I believe this setting is not severe although the equation is non-linear, so maybe the parameterization of the constitutive relation was not working well in the code.
  3. Related to the previous section, maybe the I should have set the maximum K (conductivity) by inserting conditional statements in the constitutive relation.But, if I do that, I got an error saying that “UFL conditions cannot be evaluated as bool in a Python context.”. Does this mean I need to write the parametrization in C++ syntax?

I would appreciate any help! Note I am using 2019.2.0.dev0.

import dolfin as dl
import ufl

# VGM parameters
VGM_parameters = {
'theta_r': 0.067,
'theta_s': 0.45,
'alpha': 2,
'n': 1.41,
'K_s': 0.108,
'l': 0.5}

# VGM Model for K
def K(psi, theta_r, theta_s, alpha, n, K_s, l):
    if psi >= 0:
        return K_s   # conditional statements
    m = 1-1/n
    S_e = (1 + (-alpha*psi)**n)**(-m)
    theta =  S_e * (theta_s - theta_r) + theta_r
    K = K_s*S_e**l*(1-(1-S_e**(1/m))**m)**2
    return K

n = 100 # 101 points
mesh = dl.UnitIntervalMesh(n) 
degree = 1 # linear basis functions
Vh = dl.FunctionSpace(mesh, "Lagrange", degree)

# boundry condition
def boundary_top(x, on_boundary):
    return on_boundary and abs(x[0] - 1) < dl.DOLFIN_EPS

def boundary_bot(x, on_boundary):
    return on_boundary and abs(x[0]) < dl.DOLFIN_EPS

# u_top = dl.Constant(-0.02)
# u_bot = dl.Constant(-10)

u_D = dl.Expression('-10 + 9.98*x[0]', degree = 1)  # psi = -0.02 for x = 1 (top) and psi = -10 for x = 0 (bottom)

bcs = [dl.DirichletBC(Vh, u_D, boundary_top),
       dl.DirichletBC(Vh, u_D, boundary_bot)]

u = dl.Function(Vh)
v = dl.TestFunction(Vh)
f = dl.Constant(0)
z = dl.Expression('x[0]', degree = 1)  # gravitational

# F = K(u, **VGM_parameters)*dl.inner(dl.grad(u + z), dl.grad(v))*dl.dx - f*v*dl.dx  # Cannot determine geometric dimension from expression.
F = K(u, **VGM_parameters)*dl.inner(dl.grad(u), dl.grad(v))*dl.dx - f*v*dl.dx  # unknown error


dl.solve(F == 0, u, bcs)