Hi there, I’m a complete novice at Fenics, so this may be very simple to address. The problem I am eventually trying to solve is a 1D time-dependent Shallow Ice Approximation on a tilted slope. However, right now I am trying to solve a 1D steady-state SIA on flat ground. (Gotta start somewhere, right?)
The issue I am having is in creating the variational formula. The PDE contains both the Thickness of the ice as well as the derivative of the height with respect to x.
My PDE is 0=\nabla\cdot(-BH^5(H_x)^3)
Where B is a bunch of constants and H_x is the derivative of H wrt x.
I multiply this by v and integrate by parts getting: 0=BH^5(H_x)^3\cdot v_x
My computational domain \Omega is 0 to 1000 split up into a variable number of nodes.
So I think my problem is then \int_{\Omega}\nabla \cdot (-BH^5H_x^3)(v)=\int_{\Omega} (BH^5H_x^3)(v_x)
My question is how exactly to represent these derivatives in the variational formula. Since this is in one dimension, the divergence I am removing by integrating by parts functions as a simple derivative. I have gone through the tutorial but been unable to figure this out. Any help is appreciated.
Could you carefully write the mathematics using LaTeX and the enclosed environment $ $, e.g.
$ \nabla (a + b) = \nabla a + \nabla b $
will yield
\nabla (a + b) = \nabla a + \nabla b
Also could you state the computational domain? Is it simply a 1D interval?
For analogy consider the following quick formulations: Let \Omega := (0, 1), find u \in V such that
\int_\Omega -\nabla^2 u \, v \; \mathrm{d} x = \int_\Omega f v \; \mathrm{d} x \quad \forall v \in V \\
\int_\Omega \nabla u \cdot \nabla v \; \mathrm{d} x - \int_{\partial\Omega} \nabla u \cdot \mathbf{n} \, v \; \mathrm{d} s = \int_\Omega f v \; \mathrm{d} x \quad \forall v \in V
is equivalent to
\int_\Omega \frac{\mathrm{d} u}{\mathrm{d} x} \frac{\mathrm{d} v}{\mathrm{d} x} \; \mathrm{d} x - \int_{\partial\Omega} \frac{\mathrm{d} u}{\mathrm{d} x} n_x \, v \; \mathrm{d} s = \int_\Omega f v \; \mathrm{d} x \quad \forall v \in V
Based on a glance at your weak form, other than the missing boundary terms it looks correct. Also this problem looks “very nonlinear”. If this is the very first problem you’re solving with FEniCS I would recommend starting much simpler with the Poisson problem, then the nonlinear-Poisson problem, then you’re probably in good shape to solve this one.
I’m also assuming that (H_x)^3 = H_x^3 = (\frac{\partial H}{\partial x})^3 \neq \frac{\partial^3 H}{\partial x^3}.
Yes that’s correct, it’s the first derivative raised to the third power. I have previously worked through the Poisson and non linear poisson, using the latter as my base to set up this problem. My uncertainty is in how to express the derivatives of my Function in my weak form. I wasn’t able to figure out this question from the other tutorial examples unfortunately.
Sure thing, here’s what I have so far. I have been able to progress a bit, and now the code compiles but I run into issues with newtonsolver, which I believe is what you were warning me of with the fact that my problem is very nonlinear.
from fenics import *
from mshr import *
import numpy as np
import matplotlib.pyplot as plt
#Constants
A = 10e-16 # temperature constant
rho = 911 # density
g = 9.81 # gravity
# Horizontal range
Lx = 1000 #
# nodes
Nx = 51 # Nodes along x
interval = IntervalMesh(Nx-1, 0, Lx)
# function space
V = FunctionSpace(interval, "Lagrange", 1)
#Here we set the boundary condition
H_D = Expression('100', degree=0)
#Here we set up our boundary definition
def boundary(x, on_boundary):
return on_boundary
#Here we apply our boundary definition, Dirichilet because we're specifying value
bc = DirichletBC(V, H_D, boundary)
B=(2/5)*A*(rho*g)**3
H = Function(V)
v = TestFunction(V)
H_x= grad(H)
F = (B*H**5)*dot((H_x), (H_x))*dot((H_x), grad(v))*dx
solve(F == 0, H, bc)
plot(H)
The exponent H**5 is horrific… think about it in the context of your boundary condition that H = 100 on the left and right of the domain, that value is going to be 10^{10}. I couldn’t get this to work so added the parameter H**n where 1 < n < 2 seems to work.
Further to the previous point, the gradient changes (\nabla H)^2 could easily blow out numerical precision.
Perhaps you could consider some kind of non-dimensionalisation?
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
#Constants
A = 10e-16 # temperature constant
rho = 911 # density
g = 9.81 # gravity
# Horizontal range
Lx = 1000 #
# nodes
Nx = 51 # Nodes along x
interval = IntervalMesh(Nx-1, 0, Lx)
# function space
V = FunctionSpace(interval, "Lagrange", 1)
#Here we set the boundary condition
H_D = Expression('100', degree=0)
#Here we set up our boundary definition
def boundary(x, on_boundary):
return on_boundary
#Here we apply our boundary definition, Dirichilet because we're specifying value
bc = DirichletBC(V, H_D, boundary)
B=(2/5)*A*(rho*g)**3
H = Function(V)
H_x= grad(H)
# Some non-trivial intial guess
H.interpolate(Expression("100.0 + x[0]", degree=1))
v = TestFunction(V)
n = Constant(1.0)
F = (B*H**n)*dot(H_x, H_x)*dot(H_x, grad(v))*dx
# Setup the problem and solve
solve(F == 0, H, bc)
plot(H)
plt.show()
Okay, I have managed to get it to converge. The trick was setting the bound that no negative values for H are permitted, as well as noting that since my length scale is km my initial guess for H should have a max value on the order of 1, not 100.
I will post my code below so anyone in the future with this issue can see how I solved it. Note that I changed the solver from newton to snes which requires PETSc. I also changed from the high level solve() to explicitly writing out the solver and problem to make it more clear which parameters I was setting. Changing the initial guess has small effects on the solution, but only minor.
If you have any further comments on this code, I’d be grateful, Thanks Nate for all the help.
from fenics import *
from mshr import *
import numpy as np
import matplotlib.pyplot as plt
#Constants
A = 10e-16 # temperature constant
rho = 911 # density
g = 9.81 # gravity
# Horizontal range
Lx = 1000 #
# nodes
Nx = 51
interval = IntervalMesh(Nx-1, 0, Lx)
# Define 1D function space
V = FunctionSpace(interval, "Lagrange", 1)
#Here we set the boundary condition
H_D = Expression('0', degree=0)
#Here we set up our boundary definition
def boundary(x, on_boundary):
return on_boundary
#Here we apply our boundary definition, Dirichilet because we're specifyng value
bc = DirichletBC(V, H_D, boundary)
B=(2/5)*A*(rho*g)**3
H = Function(V)
v = TestFunction(V)
H_x= grad(H)
#set bounds
lower = Function(V)
upper = Function(V)
lower = interpolate(Constant(0.0),V)
upper = interpolate(Constant(4.0),V)
#Intial guess for H
H.interpolate(Expression("1+x[0]",degree=2))
#H.interpolate(Expression("pow(x[0],.5)", degree=2))
#Create weak form
F = (B*H**5)*dot((H_x), (H_x))*dot((H_x), grad(v))*dx
#compute Jacobian
J=derivative(F, H)
problem=NonlinearVariationalProblem(F, H, bc, J=J)
solver=NonlinearVariationalSolver(problem)
problem.set_bounds(lower.vector(),upper.vector())
solver.parameters['nonlinear_solver'] = 'snes'
solver.parameters['snes_solver']['maximum_iterations']=1000
solver.parameters['snes_solver']['linear_solver']='mumps'
solver.parameters['snes_solver']['preconditioner'] = 'ilu'
(iter, converged)=solver.solve()
plot(H)
plt.show()