Transcript Baek
Motivation – Why deal.II?
Adaptive Mesh Refinement (AMR)
Advantages of AMR
Start with solving on coarse grid
Compute error
Refine mesh until error < tolerance
Same accuracy with less cost (uniformly fine meshes
are expensive!)
Achieve much higher resolution at areas of interest
Great for solving geophysical problems since the
action is localized
Flexibility
C++ program library, not software
Can be written from scratch for complex work
Well documented tutorials
www.dealii.org - step 22
Subduction Model
Subduction zone is a region
where two tectonic plates collide
and one sinks beneath the other
Compute mantle flow field and
pressure in the interior of Earth
at subduction zones by solving
Stokes flow equations
This simplified model includes
two trapezoidal geometries, each
representing the oceanic corner
and arc corner
Not coupled – two systems are
treated separately
Dip angle of 45°
Oceanic corner
Turcotte and Schubert
(2002)
Arc
corner
Stokes Flow Equations
u p f
u 0
u b
n 0
with °
in
in
on 1
on 2
Non-dimensional model
μ = 1 (uniform viscosity)
b x2 + b y2 = 1
(constant plate velocity)
Incompressible flow
No external forces, thus flow is
entirely driven by subducting
plate velocity (f = 0)
Fluid velocity (u) is set equal
to plate velocity (b) on Γ1
Natural boundary conditions of
σn = 0 everywhere else
2
u fluid velocity
p pressure
dynamic viscosity
f external forces
b plate velocity
stress
Model Geometry
u p f
u 0
u b
n 0
with °
in
in
on 1
on 2
Oceanic Corner
ub
ub
n 0
n 0
2
Arc Corner
u fluid velocity
u0
p pressure
dynamic viscosity
f external forces
b plate velocity
stress
n 0
ub
n 0
Weak Formulation
Write equations in a vector form
u p f
u 0
u b
n 0
in
in
on 1
on 2
with °
and u pI
where I is an identity matrix
2
u p f
u 0 ,
and construct a dot product with a
vector-valued test function
v
.
q
Integrating over the domain and integrating
by parts yields
(v,u) (divv, p) (q, divu) (v, f ) .
Strongly imposed Dirichlet boundary
conditions do not appear in the weak form.
Solving with Schur Complement
The weak form leads to a linear system
Ax b, which can be written as
A BT U F
P 0 ,
B 0
where A is the mass matrix on the velocity space,
B is the divergence operator, BT corressponds to the
gradient, and U, P are the values of velocity and pressure.
Solve this system by forming the Schur complement:
BA 1 BT P BA 1F,
AU F BT P.
This can be written as
A
U I
0 F
BT
0 ,
1
1 T
P
I
BA
0 BA B
such that
1
I
U A
0 F
BT
0 .
P
1
1 T
BA
I
0
BA
B
When solving for the Schur complement
it is important to have a good conditioner
for the matrix A, since it is related to the
Laplace operator and is badly conditioned.
A direct sparse LU decomposition with the
UMFPACK direct solver is used as a
preconditioner.
Flow Field in Oceanic Corner
Flow Field in Arc Corner
Pressure in Oceanic Corner
Pressure in Oceanic Corner
Pressure in Arc Corner
“Coupled” Subduction Model
u p 0
u 0
u 0
n 0
n p ext n
with ° 2
in
in
on 1
on 2
on 3
p ext pressure at boundary computed
from the oceanic corner
Coupling two sides
Solve stokes flow
equations for the “oceanic
corner” as previously done
Instead of prescribing a
moving wall condition at
the ramp (Γ3) of the “arc
corner” geometry, set the
normal stress at the
boundary equal to the
external pressure in
normal direction, which is
obtained from the
“oceanic corner” solution
“Coupled” Subduction Model
Oceanic Corner
u p 0
u 0
u 0
n 0
n p ext n
with ° 2
in
in
on 1
ub
ub
on 2
on 3
n 0
p ext pressure at boundary computed
from the oceanic corner
With this natural boundary condition
the weak form becomes
(v,u) (divv, p) (q, divu)
(v, f ) (v, p ext n) 3 .
Arc Corner
u0
n 0
n p ext n
n 0
Flow Field in Arc Corner