Solve porous flow through unsaturated soils

Looks like you get some NANs there. Maybe look at the equation one more
time? You might also want to look for terms might might be singular.

Johan

On 03/26/2012 04:00 AM, Serge-Étienne Parent wrote:
> New question #191707 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/191707
>
> Hello,
>
> I'm trying to create an adequate formulation to solve porous flow
 > through unsaturated soils (the problem is described here:
> http://ubuntuone.com/7Mew1Lrfr6ZIqZ551fnO36).
>
> The equation is highly non-linear, so I guess the solution is not straight forward. I tried to modify the file v1_np.py, replacing the default function (1+u)**m by another (van Genuchten, 1980), i.e.:
> ksat*((1-((aVG*u)**(nVG*mVG))*((1+((aVG*u)**nVG))**(-mVG)))**2)/((1+((aVG*u)**nVG))**(mVG*lVG))
>
> *** -------------------------------------------------------------------------
> *** Error: Unable to solve nonlinear system with NewtonSolver.
> *** Reason: Newton solver did not converge. Bummer.
> *** Where: This error was encountered inside NewtonSolver.cpp.
> *** -------------------------------------------------------------------------
>
> Here is the code:
>
> ###################
>
> from dolfin import *
>
> # Create mesh and define function space
> mesh = UnitCube(3, 3, 3)
> V = FunctionSpace(mesh, 'DG', 1)
>
> # Define boundary conditions
> tol = 1E-14
> def left_boundary(x, on_boundary):
> return on_boundary and abs(x[0])< tol
>
> def right_boundary(x, on_boundary):
> return on_boundary and abs(x[0]-1)< tol
>
> Gamma_0 = DirichletBC(V, Constant(0.0), left_boundary)
> Gamma_1 = DirichletBC(V, Constant(1.0), right_boundary)
> bcs = [Gamma_0, Gamma_1]
>
> # van Genuchten (1980) coefficients
> aVG=1
> nVG=0.4
> mVG=0.3
> lVG=0.5
> ksat=5e-7
>
> def q(u):
> kFn = ksat*((1-((aVG*u)**(nVG*mVG))*((1+((aVG*u)**nVG))**(-mVG)))**2)/((1+((aVG*u)**nVG))**(mVG*lVG))
> return kFn
>
> # Define variational problem
> v = TestFunction(V)
> u = TrialFunction(V)
> F = inner(q(u)*grad(u), grad(v))*dx
> u_ = Function(V) # most recently computed solution
> F = action(F, u_)
>
> J = derivative(F, u_, u)
>
> # Compute solution
> problem = NonlinearVariationalProblem(F, u_, bcs, J)
> solver = NonlinearVariationalSolver(problem)
> solver.solve()
>
>

I closely verified the equation, and it is correct. What do you mean by singular term? Undefined variables in the equation? Wrongly formulated variational form? Please forgive my ignorance, I'm quite new with FEM.

If you run the following code, you will notice how non-linear the equation is.

# Load libraries
import numpy as np
import matplotlib.pyplot as plt

# Equation parameters
aVG=3
nVG=2
mVG=1.5
lVG=0.5
ksat=5e-5

# Generated psiRange
psiminRange=0.01
psimaxRange=100
nbPoints=30
psispacer=(np.log10(psimaxRange)-np.log10(psiminRange))/(nbPoints-1)
psiRange=10**np.arange(np.log10(psiminRange),np.log10(psimaxRange+psispacer),psispacer)

# Calculate q from psiRange
q=[]
for i in range(0,len(psiRange)):
    u=psiRange[i]
    q.append(ksat*((1-((aVG*u)**(nVG*mVG))*((1+((aVG*u)**nVG))**(-mVG)))**2)/((1+((aVG*u)**nVG))**(mVG*lVG)))

# Plot with matplotlib
fig = plt.figure()
plt.loglog(psiRange,q,"-")
plt.grid(True)
plt.xlabel('Suction (m)')
plt.ylabel('Hydraulic conductivity (m/s)')
plt.show()

I tried a different function space. Instead of using dofmap=1, I tried dofmap=2, i.e. V = FunctionSpace(mesh, 'DG', 2). The solution is quite long to compute, but dolphin returns no error. I got this message though. Does it seem to be working somehow?

Solving nonlinear variational problem.
  Newton iteration 1: r (abs) = 0.000e+00 (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton solver finished in 1 iterations and 1 linear solver iterations.

Regards,

S.-É. Parent
Université Laval
Canada

On 03/26/2012 02:55 PM, Serge-Étienne Parent wrote:
> Question #191707 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/191707
>
> Status: Answered => Open
>
> Serge-Étienne Parent is still having a problem:
> I closely verified the equation, and it is correct. What do you mean by
> singular term?

If you divide by zero somewhere? If so you might want to move the mesh
slightly, and or use UFL conditionals to change the behavior of the
equation at the singularity.

> Undefined variables in the equation? Wrongly formulated
> variational form? Please forgive my ignorance, I'm quite new with FEM.

Also you can help the newton solver with an intitial guess. If you
happen to know the values at some boundary, you can solve a simple
poisson equation with the DirichletBC at these boundaries, then you feed
that solution as an initial guess into your actual solver.

Johan

> If you run the following code, you will notice how non-linear the
> equation is.
>
>
> # Load libraries
> import numpy as np
> import matplotlib.pyplot as plt
>
> # Equation parameters
> aVG=3
> nVG=2
> mVG=1.5
> lVG=0.5
> ksat=5e-5
>
> # Generated psiRange
> psiminRange=0.01
> psimaxRange=100
> nbPoints=30
> psispacer=(np.log10(psimaxRange)-np.log10(psiminRange))/(nbPoints-1)
> psiRange=10**np.arange(np.log10(psiminRange),np.log10(psimaxRange+psispacer),psispacer)
>
> # Calculate q from psiRange
> q=[]
> for i in range(0,len(psiRange)):
> u=psiRange[i]
> q.append(ksat*((1-((aVG*u)**(nVG*mVG))*((1+((aVG*u)**nVG))**(-mVG)))**2)/((1+((aVG*u)**nVG))**(mVG*lVG)))
>
> # Plot with matplotlib
> fig = plt.figure()
> plt.loglog(psiRange,q,"-")
> plt.grid(True)
> plt.xlabel('Suction (m)')
> plt.ylabel('Hydraulic conductivity (m/s)')
> plt.show()
>

Serge-Étienne, Did you manage to fix this and if so would you be able to share the solution? I too am interested in Richards' Equation.
Thanks
Myles

Here is the last version of the code, which seems to be working. I still have to implement the code to define subdomains (not working so far) and add Neumann BC and time-depency (not tried yet). Please contact me via Launchpad if you are interested to collaborate.

from dolfin import *
import numpy, sys

# Create mesh and define function space
mesh = Mesh("dyke.xml") # download dyke.xml at http://ubuntuone.com/1HvkeZDCgfB2nZRDbsMWPA
V = FunctionSpace(mesh, 'Lagrange', 1) # 1 = degree

# Define boundary conditions
def left_boundary(x, on_boundary):
    return x[0] < 24. and x[1]>=6.-DOLFIN_EPS and on_boundary

def right_boundary(x, on_boundary):
    return x[0] > 37. and x[1]==6.-DOLFIN_EPS and on_boundary

# transform boundary condition from total head to pressure head
Gamma_0 = DirichletBC(V, Expression('8.5-x[1]'), left_boundary)
Gamma_1 = DirichletBC(V, Expression('6.0-x[1]'), right_boundary)
bcs = [Gamma_0, Gamma_1]

# van Genuchten (1980) coefficients: dummy soil
aVG=10
nVG=2
mVG=1-1/nVG
lVG=0.5
ksat=3e-8

def q(u):
    kFn = ksat*((1-((aVG*u)**(nVG*mVG))*((1+((aVG*u)**nVG))**(-mVG)))**2)/((1+((aVG*u)**nVG))**(mVG*lVG))
    return kFn

# Define variational problem
v = TestFunction(V)
u = TrialFunction(V)
F = inner(q(u)*grad(u), grad(v))*dx
u_ = Function(V)
F = action(F, u_)

# J must be a Jacobian (Gateaux derivative in direction of du)
J = derivative(F, u_, u)

# Compute solution
problem = NonlinearVariationalProblem(F, u_, bcs, J)
solver = NonlinearVariationalSolver(problem)
solver.solve()

# Dump solution to file in VTK format
file = File("seep2.pvd")
file << u_
## Hold plot

# Plot solution and mesh
plot(u_, axes=True)
#plot(mesh)
interactive()