BC for a hyperbolic 1st order system of PDE


142
views
0
8 weeks ago by
Ta S  
Hi everyone,

I'm trying to solve a system of two 1st order hyperbolic PDEs with boundary conditions on parts of the boundary. I'm not sure how to implement the boundary conditions properly.

The equations I want to solve:
\[ \partial_x A + \partial_z B = 0 \]
\[ \partial_x B - \partial_z A = 0 \]
on a box with 0<x<L and 0<z<H. Boundary conditions are
\[ A(x,0) = B(x,0) = 0 \]
\[ A(L,z) = B(L,z) = 0 \]
\[ A(0,z) = z \]
Introducing the vectors
\[ a = (1,0),b = (0,1) \]
the equations above can be written as:
\[\nabla \cdot \vec{a}A + \nabla \cdot \vec{b}B = 0 \]
\[ \nabla \cdot \vec{a}B - \nabla \cdot \vec{b}A = 0 \]
I introduce two test functions v_1 and v_2 and multiply each equation with one test function, then I integrate over the domain and add the equations:
\[\int_\Omega \left(\nabla \cdot \vec{a}A + \nabla \cdot \vec{b}B \right) v_1 \mathrm{d}x + \int_\Omega \left(\nabla \cdot \vec{a}B - \nabla \cdot \vec{b}A \right) v_2 \mathrm{d}x = 0 \]

I tried solving this with a mixed function space for my two scalar functions A,B:
from fenics import *
from mshr import *

# Define domain and create mesh
H = Constant(0.1)   # in km
L = Constant(1.0) 
domain = Rectangle(Point(0, 0), Point(1,0.1))
mesh = generate_mesh(domain, 64)

# Define the function space
#P1 = FiniteElement('P', triangle, 2)
P1 = FiniteElement('DG', triangle, 1)
element = MixedElement([P1, P1])
V = FunctionSpace(mesh, element)

# Define test and trial functions
v1, v2 = TestFunctions(V)
u = Function(V)
A, B = split(u)

# Define boundaries
surface = 'near(x[1], 0)'
glacier = 'near(x[0], 1.0)'
front = 'near(x[0], 0)'

# Define front stress
front_normal_stress = Expression('-x[1]',degree=1)

# Define boundary conditions
bc_glacier = DirichletBC(V, Constant((0, 0)), glacier)
bc_surface = DirichletBC(V, Constant((0, 0)), surface)
bc_front = DirichletBC(V.sub(0), front_normal_stress, front)
bc = [bc_glacier,bc_surface,bc_front]

# Define the variational problem
a = Constant((1,0))
b = Constant((0,1))
F = (div(A*a)+div(B*b))*v1*dx + (div(B*a)-div(A*b))*v2*dx

# Compute solution
u = Function(V)
solve(F==0, u)​


The solution that I get out of this is zero everywhere. That is a solution of the PDE, but it doesn't fulfill the boundary condition at the front. What is happening? Why is FeniCS ignoring this bc?

I treid to implement the bc as a boundary integral in F after integration by parts:

F = dot(a*A + b*B,grad(v1))*dx + dot(a*B - b*A,grad(v2))*dx + front_normal_stress*v1*ds

Then I get the following error message:
*** Error:   Unable to successfully call PETSc function 'MatSetValuesLocal'.
*** Reason:  PETSc error code is: 63 (Argument out of range).
*** Where:   This error was encountered inside /tmp/dolfin/dolfin/la/PETScMatrix.cpp.​
One obvios problem that I see with this approach is that I have no bc at all at the bottom boundary, but in this case FeniCs will try to calculate the boundary integral there as well because it doesn't know this.

Any ideas?
Community: FEniCS Project

1 Answer


2
8 weeks ago by
pf4d  
First, you have not imposed your boundary conditions in the "solve" method, so you were solving the zero-Neumann problem, which is zero everywhere for your BVP.

Next, you don't need to specify the unit-vector stuff, you can specify the weak form directly from the PDE.

Next, you redefined "u" prior to calling solve, so you were solving for some variable that was not even in the weak form.

Finally, you have a typo in your math above.  The code specifies the boundary conditions
\[
\begin{align}
A(x, z=0) & = B(x, z=0) = 0 \\
A(x=1, z) & = B(x=1, z) = 0 \\
B(x=0, z) &= -z
\end{align}
\]

Here is a solution that works :

from fenics import *
from mshr   import generate_mesh, Rectangle

# Define domain and create mesh
H       = 0.1
L       = 1.0
domain  = Rectangle(Point(0, 0), Point(L, H))
mesh    = generate_mesh(domain, 64)

# Define the function space
P1      = FiniteElement('CG', triangle, 1)
element = MixedElement([P1, P1])
V       = FunctionSpace(mesh, element)

# Define test and trial functions
v1, v2 = TestFunctions(V)
phi    = TrialFunction(V)
u      = Function(V)
A, B   = u

# Define boundaries
def surface(x, on_boundary):  return on_boundary and near(x[1], 0)
def glacier(x, on_boundary):  return on_boundary and near(x[0], L)
def front(x, on_boundary):    return on_boundary and near(x[0], 0)

# Define front stress
front_normal_stress = Expression('-x[1]', degree=1)

# Define boundary conditions
bc_surface = DirichletBC(V, Constant((0, 0)), surface)
bc_glacier = DirichletBC(V, Constant((0, 0)), glacier)
bc_front   = DirichletBC(V.sub(0), front_normal_stress, front)
bc         = [bc_glacier, bc_surface,bc_front]

# Define the variational problem
F = + (A.dx(0) + B.dx(1))*v1*dx \
    + (B.dx(0) - A.dx(1))*v2*dx

# define the Jacobian :
J = derivative(F, u, phi)

# compute solution using Newton's method :
solve(F==0, u, bc, J=J)

print u.vector().min(), u.vector().max()

File('u.pvd') << u
​

I hope I didn't just do your glaciological-modeling homework for you.

Thanks a lot!

The missing bc in the solve call is a typo, I was trying a lot of different things and in my original attempts, I had it included. 
I saw u being redefined before solving in many examples I looked at, but that probably only makes sense if you've initialised it as a TrialFunction first, right?
I didn't know I could use partial derivatives, because all the examples in the manual use nabla operators, so thanks for that.

Anyway, I'm really glad it works now. I'm new to FeniCS and I spent another week triyng to get another FEM package to work that apparently couldn't solve hyperbolic equations at all... So I'm really happy I can get back to doing actual physics now :)

Thanks a lot and greetings from PIK
written 8 weeks ago by Ta S  
You are welcome and correct regarding the TrialFunction.

What is PIK?
written 8 weeks ago by pf4d  
I saw that you are at AWI in Bremerhaven, I'm at the Potsdam Institude of Climate Impact Research ;-)
written 8 weeks ago by Ta S  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »