Commit 59873450 authored by Samuel Burbulla's avatar Samuel Burbulla
Browse files

Add mesh files.

parent 481819a8
name = 'grid.msh'
import sys
if len(sys.argv) > 1:
name = sys.argv[1]
else:
name = 'grid.msh'
tikz = []
readNodes = False
......@@ -37,8 +42,8 @@ for line in file.readlines():
# interface edges
if len(data) == 7:
_, _, _, _, p, a, b = line.split(' ')
if int(p) >= 10:
_, _, _, o, p, a, b = line.split(' ')
if int(o) >= 10 or int(p) >= 10:
tikz += ['\draw[very thick] ('+a+') -- ('+b+');']
out = open(name+'.tikz', 'w')
......
......@@ -27,7 +27,7 @@ two indices :math:`i` and :math:`j` that indicate the two vertices of the edge.
\tikzset{edge/.style={midway, sloped, circle, fill=white, inner sep=1pt}}
\draw[thick] (0,0) node[vertex] {\tiny 0}
-- (2,0) node[vertex] {\tiny 1} node[edge] {\tiny 2(1)}
-- (2,0) node[vertex] {\tiny 1} node[edge] {\tiny 2(0)}
-- (0,2) node[vertex] {\tiny 2} node[edge] {\tiny 0(2)}
-- (0,0) node[edge] {\tiny 1(1)};
......
......@@ -57,6 +57,7 @@ breathe_default_project = 'dune-mmesh'
tikz_proc_suite = 'GhostScript'
tikz_tikzlibraries = 'calc'
tikz_resolution = 2000
nbsphinx_execute = 'never'
......@@ -70,7 +71,7 @@ html_theme = "sphinx_rtd_theme"
# Add any paths that contain custom static files (such as style sheets) here,
# relative to this directory. They are copied after the builtin static files,
# so a file named "default.css" will overwrite the builtin "default.css".
html_static_path = ['_static']
# html_static_path = ['_static']
# -- Options for Latex output -------------------------------------------------
latex_documents = [
......
......@@ -16,7 +16,13 @@ We implemented a few examples to display how Dune-MMesh can be used in different
examples/poroelasticity
examples/twophase
.. toctree::
:maxdepth: 1
:hidden:
examples/grids/grids
Some examples for the creation of grid files can be found in :ref:`grids`.
.. [CMR+18] C. Chalons, J. Magiera, C. Rohde, M. Wiebe. A Finite-Volume Tracking Scheme for Two-Phase Compressible Flow. Theory, Numerics and Applications of Hyperbolic Problems I, pp. 309--322, 2018.
......
......@@ -17,7 +17,7 @@
"source": [
"## Grid creation\n",
"\n",
"We use a .msh file that contains an interface $\\Gamma = [0.25, 0.75] \\times {0.5}$ embedded in a domain $\\Omega = [0,1]^2$."
"We use the [horizontal](grids/horizontal.rst) grid file that contains an interface $\\Gamma = [0.25, 0.75] \\times {0.5}$ embedded in a domain $\\Omega = [0,1]^2$."
]
},
{
......
%% Cell type:markdown id:organic-pierre tags:
# Coupling to the interface
This is an example of how to solve coupled problems on the bulk and interface grid.
%% Cell type:markdown id:ambient-charge tags:
## Grid creation
We use a .msh file that contains an interface $\Gamma = [0.25, 0.75] \times {0.5}$ embedded in a domain $\Omega = [0,1]^2$.
We use the [horizontal](grids/horizontal.rst) grid file that contains an interface $\Gamma = [0.25, 0.75] \times {0.5}$ embedded in a domain $\Omega = [0,1]^2$.
%% Cell type:code id:round-motel tags:
``` python
from dune.grid import reader
from dune.mmesh import mmesh
dim = 2
file = "grids/horizontal.msh"
gridView = mmesh((reader.gmsh, file), dim)
igridView = gridView.hierarchicalGrid.interfaceGrid
```
%% Cell type:markdown id:equivalent-spanking tags:
## Solve a problem on the bulk grid
%% Cell type:markdown id:general-persian tags:
Let us solve the Poisson equation
\begin{align}
-\Delta u = f & \qquad \text{in}\ \Omega
\end{align}
on the bulk grid and on the interface grid. We use the manufactured solution
\begin{align}
u &= \sin(4 \pi x y), \\
f &= -\operatorname{div}( \nabla u ).
\end{align}
%% Cell type:code id:bearing-lighting tags:
``` python
from ufl import *
from dune.ufl import DirichletBC
from dune.fem.space import lagrange
from dune.fem.scheme import galerkin
from dune.fem.function import integrate
space = lagrange(gridView, order=3)
u = TrialFunction(space)
v = TestFunction(space)
x = SpatialCoordinate(space)
exact = sin(x[0]*x[1]*4*pi)
f = -div(grad(exact))
a = inner(grad(u), grad(v)) * dx
b = f * v * dx
scheme = galerkin([a == b, DirichletBC(space, exact)], solver=("suitesparse", "umfpack"))
uh = space.interpolate(0, name="solution")
scheme.solve(target=uh)
def L2(u1, u2):
return sqrt(integrate(u1.grid, (u1-u2)**2, order=5))
L2(uh, exact)
```
%%%% Output: execute_result
5.628259764335839e-07
%% Cell type:markdown id:combined-hurricane tags:
## Solve a problem on the interface
We can solve similar problem on the interface $\Gamma$.
%% Cell type:code id:through-munich tags:
``` python
ispace = lagrange(igridView, order=3)
iuh = ispace.interpolate(0, name="isolution")
iu = TrialFunction(ispace)
iv = TestFunction(ispace)
ix = SpatialCoordinate(ispace)
iexact = sin(0.5*ix[dim-2]*4*pi)
iF = -div(grad(iexact))
ia = inner(grad(iu), grad(iv)) * dx
ib = iF * iv * dx
ischeme = galerkin([ia == ib, DirichletBC(ispace, iexact)])
ischeme.solve(target=iuh)
L2(iuh, iexact)
```
%%%% Output: execute_result
5.807742030532398e-08
%% Cell type:markdown id:joint-workstation tags:
We can use the `plotPointData` function to visualize the solution of both grids.
%% Cell type:code id:experimental-hopkins tags:
``` python
import matplotlib.pyplot as plt
from dune.fem.plotting import plotPointData as plot
figure = plt.figure(figsize=(3,3))
plot(uh, figure=figure, gridLines=None)
plot(iuh, figure=figure, linewidth=0.04, colorbar=None)
```
%%%% Output: display_data
[Hidden Image Output]
%% Cell type:markdown id:electric-taiwan tags:
## Couple bulk to surface
We are able to compute traces of discrete functions on $\Omega$ along $\Gamma$.
%% Cell type:code id:union-independence tags:
``` python
from dune.mmesh import trace
tr = avg(trace(uh))
ib = inner(grad(tr), grad(iv)) * dx
iuh.interpolate(0)
ischeme = galerkin([ia == ib, DirichletBC(ispace, avg(trace(uh)))])
ischeme.solve(target=iuh)
L2(iuh, iexact)
```
%%%% Output: execute_result
4.266679509138358e-08
%% Cell type:markdown id:invisible-certificate tags:
## Couple surface to bulk
Similarly, we can evaluate a discrete function on $\Gamma$ at the skeleton of the triangulation of $\Omega$.
%% Cell type:code id:inside-workshop tags:
``` python
from dune.mmesh import skeleton
sk = skeleton(iuh)
b = avg(sk) * avg(v) * dS
uh.interpolate(0)
scheme = galerkin([a == b, DirichletBC(space, 0)])
scheme.solve(target=uh)
figure = plt.figure(figsize=(3,3))
plot(uh, figure=figure, gridLines=None)
```
%%%% Output: display_data
[Hidden Image Output]
%% Cell type:markdown id:exposed-shadow tags:
## Compute jump of gradient of traces
We can use trace and skeleton within UFL expressions. This is useful when using jump and average operators, but also to compute gradients and more.
%% Cell type:code id:analyzed-change tags:
``` python
from dune.mmesh import normals
n = normals(igridView)
jmp = jump(grad(trace(uh)), n)
```
......
This diff is collapsed.
**********
circle.geo
**********
.. include:: circle.geo
:literal:
This file can be meshed using :code:`gmsh -2 -format msh2 circle.geo`.
.. include:: circle.msh.rst
.. _grids:
**********
Grid files
**********
Some examples for the creation of grid files:
.. toctree::
circle
horizontal
tjunction
vertical
This diff is collapsed.
*************
horizontal.py
*************
.. include:: horizontal.py
:literal:
.. include:: horizontal.msh.rst
This diff is collapsed.
# A rectangle grid with a horizontal centered interface
# A rectangle grid with a T-shaped junction
name = "tjunction.msh"
h = 0.1
......
************
tjunction.py
************
.. include:: tjunction.py
:literal:
.. include:: tjunction.msh.rst
This diff is collapsed.
# A rectangle grid with a vertical interface
def create(name, h, hf=None):
if hf is None:
hf = h
import gmsh
gmsh.initialize()
gmsh.option.setNumber("General.Verbosity", 0)
gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)
gmsh.model.add(name)
p1 = gmsh.model.geo.addPoint( 0, 0, 0, h, 1)
p2 = gmsh.model.geo.addPoint( 0.5, 0, 0, hf, 2)
p3 = gmsh.model.geo.addPoint( 1, 0, 0, h, 3)
p4 = gmsh.model.geo.addPoint( 1, 1, 0, h, 4)
p5 = gmsh.model.geo.addPoint( 0.5, 1, 0, hf, 5)
p6 = gmsh.model.geo.addPoint( 0, 1, 0, h, 6)
l1 = gmsh.model.geo.addLine(p1, p2, 1)
l2 = gmsh.model.geo.addLine(p2, p3, 2)
l3 = gmsh.model.geo.addLine(p3, p4, 3)
l4 = gmsh.model.geo.addLine(p4, p5, 4)
l5 = gmsh.model.geo.addLine(p5, p6, 5)
l6 = gmsh.model.geo.addLine(p6, p1, 6)
lf = gmsh.model.geo.addLine(p2, p5, 10)
gmsh.model.geo.addCurveLoop([1, 10, 5, 6], 1)
gmsh.model.geo.addPlaneSurface([1], 0)
gmsh.model.geo.addCurveLoop([2, 3, 4, -10], 2)
gmsh.model.geo.addPlaneSurface([2], 1)
gmsh.model.geo.synchronize()
gmsh.model.mesh.generate(dim=2)
gmsh.write(name)
gmsh.finalize()
name = "vertical.msh"
h = 0.1
hf = 0.1
import gmsh
gmsh.initialize()
gmsh.option.setNumber("General.Verbosity", 0)
gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)
gmsh.model.add(name)
p1 = gmsh.model.geo.addPoint( 0, 0, 0, h, 1)
p2 = gmsh.model.geo.addPoint( 0.5, 0, 0, hf, 2)
p3 = gmsh.model.geo.addPoint( 1, 0, 0, h, 3)
p4 = gmsh.model.geo.addPoint( 1, 1, 0, h, 4)
p5 = gmsh.model.geo.addPoint( 0.5, 1, 0, hf, 5)
p6 = gmsh.model.geo.addPoint( 0, 1, 0, h, 6)
l1 = gmsh.model.geo.addLine(p1, p2, 1)
l2 = gmsh.model.geo.addLine(p2, p3, 2)
l3 = gmsh.model.geo.addLine(p3, p4, 3)
l4 = gmsh.model.geo.addLine(p4, p5, 4)
l5 = gmsh.model.geo.addLine(p5, p6, 5)
l6 = gmsh.model.geo.addLine(p6, p1, 6)
lf = gmsh.model.geo.addLine(p2, p5, 10)
gmsh.model.geo.addCurveLoop([1, 10, 5, 6], 1)
gmsh.model.geo.addPlaneSurface([1], 0)
gmsh.model.geo.addCurveLoop([2, 3, 4, -10], 2)
gmsh.model.geo.addPlaneSurface([2], 1)
gmsh.model.geo.synchronize()
gmsh.model.mesh.generate(dim=2)
gmsh.write(name)
gmsh.finalize()
***********
vertical.py
***********
.. include:: vertical.py
:literal:
.. include:: vertical.msh.rst
......@@ -17,7 +17,7 @@
"source": [
"## Grid creation\n",
"\n",
"We use a .msh file that contains an interface $\\Gamma = {0.5} \\times [0, 1]$ embedded in a domain $\\Omega = [0,1]^2$ and construct an _adaptive_ leaf grid view."
"We use the [vertical](grids/vertical.rst) grid file that contains an interface $\\Gamma = {0.5} \\times [0, 1]$ embedded in a domain $\\Omega = [0,1]^2$ and construct an _adaptive_ leaf grid view."
]
},
{
......
%% Cell type:markdown id:assured-parts tags:
# Moving and adapting
This is an example of how to move the interface and adapt the mesh.
%% Cell type:markdown id:racial-mandate tags:
## Grid creation
We use a .msh file that contains an interface $\Gamma = {0.5} \times [0, 1]$ embedded in a domain $\Omega = [0,1]^2$ and construct an _adaptive_ leaf grid view.
We use the [vertical](grids/vertical.rst) grid file that contains an interface $\Gamma = {0.5} \times [0, 1]$ embedded in a domain $\Omega = [0,1]^2$ and construct an _adaptive_ leaf grid view.
%% Cell type:code id:spread-indian tags:
``` python
from dune.grid import reader
from dune.mmesh import mmesh
from dune.fem.view import adaptiveLeafGridView as adaptive
dim = 2
file = "grids/vertical.msh"
gridView = adaptive( mmesh((reader.gmsh, file), dim) )
hgrid = gridView.hierarchicalGrid
igridView = hgrid.interfaceGrid
```
%% Cell type:markdown id:permanent-lesson tags:
## Problem
Let us consider the following transport problem.
\begin{align}
u_t + \operatorname{div} f(u) = 0, & \qquad \text{in } \Omega \times [0,T], \\
u(\cdot, 0) = u_0, & \qquad \text{in } \Omega
\end{align}
where
\begin{align}
f(u) &= [1,0]^T u, \\
u_0(x,y) &= (0.5+x) \chi_{x<0.5} \\
\text{ and } T &= 0.2.
\end{align}
The interface is supposed to move with the transport speed, i.e. $s = [1,0]^T$.
\begin{align*}
\renewcommand{\jump}[1]{[\mskip-5mu[ #1 ]\mskip-5mu]}
\end{align*}
%% Cell type:code id:younger-howard tags:
``` python
from ufl import *
from dune.ufl import Constant
t = 0
tEnd = 0.4
dt = 0.04
def speed():
return as_vector([1.0, 0.0])
def movement(x):
return as_vector([1.0, 0.0])
def f(u):
return speed() * u
def u0(x):
return conditional(x[0] < 0.5, 0.5+x[0], 0.0)
def uexact(x, t):
return u0( x - t * speed() )
```
%% Cell type:markdown id:homeless-coordinate tags:
## Finite Volume Moving Mesh Method
We use a Finite Volume Moving Mesh method to keep the discontinuity sharp. It can be formulated by
\begin{align}
\int_\Omega (u^{n+1} |det(\Psi)| - u^n) v\ dx + \Delta t \int_\mathcal{F} \big( g(u^{n+1}, n) - h(u^{n+1}, n) \big) \jump{v}\ dS = 0
\end{align}
where $\Psi = x + \Delta t s$.
The numerical fluxes $g(u, n)$ and $h(u, n)$ are assumed to be consistent with the flux functions $f(u) \cdot n$ and $u (s \cdot n)$, respectively.
%% Cell type:code id:encouraging-harvard tags:
``` python
from dune.fem.space import finiteVolume
space = finiteVolume(gridView)
u = TrialFunction(space)
v = TestFunction(space)
x = SpatialCoordinate(space)
n = FacetNormal(space)
uh = space.interpolate(u0(x), name="uh")
uh_old = uh.copy()
```
%% Cell type:code id:fourth-mortality tags:
``` python
import numpy as np
from dune.geometry import vertex
from dune.mmesh import edgeMovement
def getShifts():
mapper = igridView.mapper({vertex: 1})
shifts = np.zeros((mapper.size, dim))
for v in igridView.vertices:
shifts[ mapper.index(v) ] = as_vector(movement( v.geometry.center ))
return shifts
em = edgeMovement(gridView, getShifts())
time = Constant(t, name="time")
def g(u, n):
sgn = inner(speed(), n('+'))
return inner( conditional( sgn > 0, f( u('+') ), f( u('-') ) ), n('+') )
def gBnd(u, n):
sgn = inner(speed(), n)
return inner( conditional( sgn > 0, f(u), f(uexact(x, time)) ), n )
def h(u, n):
sgn = inner(em('+'), n('+'))
return conditional( sgn > 0, sgn * u('+'), sgn * u('-') )
```
%% Cell type:code id:superb-passage tags:
``` python
from dune.fem.scheme import galerkin
tau = Constant(dt, name="tau")
detPsi = abs(det(nabla_grad(x + tau * em)))
a = (u * detPsi - uh_old) * v * dx
a += tau * (g(u, n) - h(u, n)) * jump(v) * dS
a += tau * gBnd(u, n) * v * ds
scheme = galerkin([a == 0], solver=("suitesparse","umfpack"))
```
%% Cell type:markdown id:ideal-bradley tags:
## Timeloop
Now, we can start with the time loop. We need to set the `fem.adaptation.method` parameter to `callback` in order to use the adaptation startegy of Dune-MMesh.
%% Cell type:code id:alive-responsibility tags:
``` python
from dune.fem import parameter, adapt
parameter.append( { "fem.adaptation.method": "callback" } )
from dune.fem.plotting import plotPointData as plot
import matplotlib.pyplot as plt
fig, axs = plt.subplots(1, 5, figsize=(12,3))
i = 0
while t < tEnd:
hgrid.markElements()
hgrid.ensureInterfaceMovement(getShifts()*dt)
adapt([uh])
hgrid.moveInterface(getShifts()*dt)
em.assign(edgeMovement(gridView, getShifts()))
t += dt
time.assign(t)
uh_old.assign(uh)
scheme.solve(target=uh)
i += 1
if i % 2 == 0:
plot(uh, figure=(fig, axs[i//2-1]), clim=[0,1], colorbar=None)
```
%%%% Output: display_data
[Hidden Image Output]
......
......@@ -31,7 +31,7 @@
"metadata": {},
"source": [
"## Test problem\n",
"We consider the domains $\\Omega_1 = [0,0.5] \\times [0,1]$, $\\Omega_2 = [0.5,1] \\times [0,1]$ and $\\Gamma = \\{ 0.5 \\} \\times [0,1]$ and the physical parameters $\\mathbf{K} = \\mathbf{I}$, $\\mathbf{K}_\\parallel^\\Gamma = \\mathbf{I}$, $K_\\perp^\\Gamma = 2$ and $d = 1$."
"We consider the [vertical](grids/vertical.rst) grid file defining the domains $\\Omega_1 = [0,0.5] \\times [0,1]$, $\\Omega_2 = [0.5,1] \\times [0,1]$ and $\\Gamma = \\{ 0.5 \\} \\times [0,1]$ and the physical parameters $\\mathbf{K} = \\mathbf{I}$, $\\mathbf{K}_\\parallel^\\Gamma = \\mathbf{I}$, $K_\\perp^\\Gamma = 2$ and $d = 1$."
]
},
{
......
%% Cell type:markdown id: tags:
# Single-phase flow in fractured porous media
%% Cell type:markdown id: tags:
Consider the following mixed-dimensional model for single-phase flow in fractured porous media.
\begin{align}
\renewcommand{\jump}[1]{[\mskip-5mu[ #1 ]\mskip-5mu]}
\renewcommand{\avg}[1]{\{\mskip-5mu\{ #1 \}\mskip-5mu\}}
\end{align}
\begin{align}
-\nabla\cdot\left(\mathbf{K}_i\nabla p_i\right) &= q_i &\qquad &\text{in } \Omega_i, \quad &i&=1,2 ,\\
-\nabla_\parallel \cdot \left( d\mathbf{K}^\Gamma_\parallel \nabla_\parallel \left( d p^\Gamma\right) \right) &= q^\Gamma - \jump{\mathbf{K}\nabla p} &\qquad &\text{in }\Gamma ,\\
-\avg{\mathbf{K} \nabla p} \cdot {n} &= K^\Gamma_\perp \left( p_1 - p_2 \right) &\qquad &\text{on }\Gamma ,\\
-\jump{\mathbf{K} \nabla p} &= 12 K_\perp^\Gamma \left( \avg{p} - p^\Gamma \right) &\qquad &\text{on } \Gamma ,\\
p_i &= g_i &\qquad &\text{on } \rho_i , &i&=1,2 , \\
p^\Gamma &= g^\Gamma &\qquad &\text{on } \partial\Gamma. \\
\end{align}
%% Cell type:markdown id: tags:
## Test problem
We consider the domains $\Omega_1 = [0,0.5] \times [0,1]$, $\Omega_2 = [0.5,1] \times [0,1]$ and $\Gamma = \{ 0.5 \} \times [0,1]$ and the physical parameters $\mathbf{K} = \mathbf{I}$, $\mathbf{K}_\parallel^\Gamma = \mathbf{I}$, $K_\perp^\Gamma = 2$ and $d = 1$.
We consider the [vertical](grids/vertical.rst) grid file defining the domains $\Omega_1 = [0,0.5] \times [0,1]$, $\Omega_2 = [0.5,1] \times [0,1]$ and $\Gamma = \{ 0.5 \} \times [0,1]$ and the physical parameters $\mathbf{K} = \mathbf{I}$, $\mathbf{K}_\parallel^\Gamma = \mathbf{I}$, $K_\perp^\Gamma = 2$ and $d = 1$.
%% Cell type:code id: tags:
``` python
from dune.grid import reader
from dune.mmesh import mmesh, domainMarker
from ufl import *
dim = 2
file = "grids/vertical.msh"
gridView = mmesh((reader.gmsh, file), dim)
igridView = gridView.hierarchicalGrid.interfaceGrid
K = Identity(dim)
K_g = Identity(dim)
K_perp = 2
d = 1
```
%% Cell type:markdown id: tags:
We use a manufactured solution given by
\begin{align*}
p_1 \left( {x} \right) &= \sin\left( 4x_1 \right) \cos\left( \pi x_2 \right) ,\\
p_2 \left( {x} \right) &= \cos\left( 4x_1 \right) \cos\left( \pi x_2 \right) ,\\
p^\Gamma \left( {s}\right) &= \tfrac{3}{4} \left( \cos\left( 2\right) + \sin\left( 2\right) \right) \cos\left( \pi x_2 \right)
\end{align*}
and the source terms
\begin{align*}
q_1 \left( {x} \right) &= \left( 16 + \pi^2 \right) \sin\left( 4x_1 \right) \cos\left( \pi x_2 \right) ,\\
q_2 \left( {x} \right) &= \left( 16 + \pi^2 \right) \cos\left( 4x_1 \right) \cos\left( \pi x_2 \right),\\
q^\Gamma \left( {s} \right) &= \left( \cos\left( 2\right) + \sin\left( 2 \right) \right) \left( 4 + \tfrac{3}{4} d^2 \pi^2 \right) \cos\left( \pi s_2 \right).
\end{align*}
%% Cell type:code id: tags:
``` python
dm = domainMarker(gridView)
x = SpatialCoordinate( Cell("triangle", dim) )
x_g = SpatialCoordinate( Cell("interval", dim) )
p1 = sin(4*x[0])*cos(pi*x[1])
p2 = cos(4*x[0])*cos(pi*x[1])
pexact = dm*p1 + (1-dm)*p2
p_gexact = 3./4.*(cos(2.)+sin(2.))*cos(pi*x_g[1])
q = (16.+pi**2) * pexact
q_g = 4./3.*p_gexact*(4. + 3./4.* d**2 * pi**2)
g = pexact
g_g = p_gexact
```
%% Cell type:markdown id: tags:
## Mixed-dimensional interior-penalty DG formulation
An IPDG formulation for the problem above can be written as:
Find $( p_h , p_h^\Gamma ) \in \mathcal{S}_h^{PM} \times \mathcal{S}_h^\Gamma$, s.t.
\begin{align}
\mathcal{B}_{DG} \big( ( p_h , p_h^\Gamma ) , ( \phi_h , \phi_h^\Gamma )\big) = \mathcal{L}_{DG} ( \phi_h , \phi_h^\Gamma )
\end{align}
for all $( \phi_h , \phi_h^\Gamma ) \in \mathcal{S}_h^{PM} \times \mathcal{S}_h^\Gamma$, where
\begin{align}
\begin{split}
\mathcal{B}_{DG} \big( ( p_h , p_h^\Gamma ) , ( \phi_h , \phi_h^\Gamma )\big) &:= \mathcal{B}_{DG}^{PM} ( p_h , \phi_h ) + \mathcal{B}^\Gamma_{DG} ( p_h^\Gamma , \phi_h^\Gamma ) + \mathcal{K}_{DG} \big( ( p_h, p_h^\Gamma ) , ( \phi_h , \phi_h^\Gamma )\big), \\
\mathcal{L}_{DG} ( \phi_h , \phi_h^\Gamma ) &:= \mathcal{L}^{PM}_{DG} ( \phi_h ) + \mathcal{L}^\Gamma_{DG} ( \phi_h^\Gamma ),
\end{split}
\end{align}
where the forms are defined below.
%% Cell type:code id: tags:
``` python
from dune.fem.space import dglagrange
order = 2
space = dglagrange( gridView, order=order)
space_g = dglagrange(igridView, order=order)
p = TrialFunction(space)
phi = TestFunction(space)
p_g = TrialFunction(space_g)
phi_g = TestFunction(space_g)
n = FacetNormal(space)
n_g = FacetNormal(space_g)
from dune.mmesh import interfaceIndicator
I = interfaceIndicator(igridView)
```
%% Cell type:markdown id: tags:
## Bulk forms
\begin{align}
\mathcal{B}^{PM}_{DG} \left( p_h , \phi_h \right) &:= \int_\Omega \mathbf{K} \nabla_h p_h \cdot \nabla_h \phi_h \,{d} V \\
&\quad + \int_{\mathcal{F}^{\circ}\setminus\Gamma} \left( \mu \jump{p_h}\cdot\jump{\phi_h} - \avg{\mathbf{K} \nabla_h p_h} \cdot \jump{\phi_h} - \jump{p_h} \cdot \avg{\mathbf{K} \nabla_h \phi_h }\right) \,{d}\sigma \\
&\quad + \int_{\partial\Omega} \mu p_h \phi_h \,{d}\sigma - \int_{\partial\Omega} \left( p_h \mathbf{K} \nabla_h \phi_h + \phi_h \mathbf{K} \nabla_h p_h \right) \cdot {d}{\sigma}, \\
\mathcal{L}^{PM}_{DG} \left( \phi_h\right) &:= \int_\Omega q \phi_h \,{d} V +\int_{\partial\Omega} \mu g \phi_h \,{d}\sigma - \int_{\partial\Omega} g \mathbf{K} \nabla_h \phi_h \cdot {d}{\sigma}, \\
\end{align}
%% Cell type:code id: tags:
``` python
mu = 1e3
B = dot(dot(grad(p), K), grad(phi)) * dx
B += mu * jump(p) * jump(phi) * (1-I)*dS
B -= jump(phi) * dot(avg(dot(grad( p ), K)), n('+')) * (1-I)*dS
B -= jump( p ) * dot(avg(dot(grad(phi), K)), n('+')) * (1-I)*dS
B += mu * p * phi * ds
B -= p * dot(dot(grad(phi), K), n) * ds
B -= phi * dot(dot(grad(p), K), n) * ds
L = q * phi * dx
L += mu * g * phi * ds
L -= g * dot(dot(grad(phi), K), n) * ds
```
%% Cell type:markdown id: tags:
## Interface forms
\begin{align}
\mathcal{B}^\Gamma_{DG} \left( p_h^\Gamma , \phi_h^\Gamma \right) &:= \int_\Gamma \mathbf{K}_\parallel^\Gamma \nabla_{\parallel h}\left( d p_h^\Gamma \right) \cdot \nabla_{\parallel h} \left( d \phi_h^\Gamma \right) \,{d}\sigma
+ \int_{\mathcal{E}_\Gamma^\circ} d \Big[ \mu^\Gamma d \jump{ p_h^\Gamma} \cdot \jump{\phi_h^\Gamma} \\
&\quad - \avg{\mathbf{K}_\parallel^\Gamma \nabla_{\parallel h} \left( d p_h^\Gamma \right) } \cdot \jump{\phi_h^\Gamma} - \jump{p_h^\Gamma} \cdot \avg{\mathbf{K}_\parallel^\Gamma \nabla_{\parallel h}\left( d \phi_h^\Gamma\right) } \Big] \,{d}r \\
&\quad + \int_{\partial\Gamma} \mu^\Gamma d^2 p_h^\Gamma \phi_h^\Gamma \,{d} r - \int_{\partial \Gamma} d \left[ p_h^\Gamma \mathbf{K}_\parallel^\Gamma \nabla_{\parallel h} \left( d \phi_h^\Gamma \right) + \phi_h^\Gamma \mathbf{K}_\parallel^\Gamma \nabla_{\parallel h} \left( d p_h^\Gamma\right)\right] \cdot {d}{r},\\
\mathcal{L}_{DG}^\Gamma \left( \phi_h^\Gamma \right) &:= \int_\Gamma q^\Gamma \phi_h^\Gamma \,{d}\sigma + \int_{\partial\Gamma} \mu^\Gamma d^2 g^\Gamma \phi_h^\Gamma \,{d} r - \int_{\partial\Gamma} d g^\Gamma \mathbf{K}_\parallel^\Gamma \nabla_{\parallel h} \left( d \phi_h^\Gamma \right) \cdot {d} {r}, \\
\end{align}
%% Cell type:code id: tags:
``` python
B_g = d * dot(dot(grad(d*p_g), K_g), grad(d*phi_g)) * dx
B_g += d * mu * d * jump(p_g) * jump(phi_g) * dS
B_g -= d * jump(phi_g) * dot(avg(dot(grad(d * p_g ), K_g)), n_g('+')) * dS
B_g -= d * jump( p_g ) * dot(avg(dot(grad(d * phi_g), K_g)), n_g('+')) * dS
B_g += mu * d**2 * p_g * phi_g * ds
B_g -= d**2 * p_g * dot(dot(grad(d * phi_g), K_g), n_g) * ds
B_g -= d**2 * phi_g * dot(dot(grad(d * p_g), K_g), n_g) * ds
L_g = d * q_g * phi_g * dx
L_g += mu * d**2 * g_g * phi_g * ds
L_g -= d**2 * g_g * dot(dot(grad(d * phi_g), K_g), n_g) * ds
```
%% Cell type:markdown id: tags:
## Coupling form
\begin{align}
\mathcal{K}_{DG} \left( \left( p_h , p_h^\Gamma \right) , \left( \phi_h , \phi_h^\Gamma \right) \right) & := \int_\Gamma K_\perp^\Gamma \jump{p_h} \cdot \jump{\phi_h} \,{d} \sigma + \int_\Gamma \beta^\Gamma \left( \avg{p_h} - p_h^\Gamma \right) \left( \avg{\phi_h} - \phi_h^\Gamma \right) \,{d}\sigma .
\end{align}
%% Cell type:code id: tags:
``` python
from dune.mmesh import trace, skeleton
ph = space.interpolate(0, name="pressure")
ph_g = space_g.interpolate(0, name="pressure_g")
xi = 2./3.
beta = 4. * K_perp / (2. * xi - 1)
K_DG = K_perp * jump(p) * jump(phi) * I*dS
K_DG += beta * ( avg(p) * avg(phi) - avg(skeleton(ph_g)) * avg(phi) ) * I*dS
K_DG_g = beta * ( -avg(trace(ph)) * phi_g + p_g * phi_g ) * dx
```
%% Cell type:markdown id: tags:
## Solve
We use the iterative solution strategy provided by Dune-MMesh.
%% Cell type:code id: tags:
``` python
from dune.fem.scheme import galerkin
from dune.mmesh import iterativeSolve
scheme = galerkin([B + K_DG == L ], solver=("suitesparse", "umfpack"))
scheme_g = galerkin([B_g + K_DG_g == L_g], solver=("suitesparse", "umfpack"))
res = iterativeSolve(schemes=(scheme, scheme_g), targets=(ph, ph_g), verbose=True)
```
%%%% Output: stream
0: [ 2.36e+02 1.62e+00 ]
1: [ 4.08e+00 6.62e-02 ]
2: [ 2.34e-01 4.01e-03 ]
3: [ 1.43e-02 2.46e-04 ]
4: [ 8.77e-04 1.51e-05 ]
5: [ 5.38e-05 9.27e-07 ]
6: [ 3.30e-06 5.69e-08 ]
7: [ 2.03e-07 3.49e-09 ]
8: [ 1.25e-08 2.14e-10 ]
9: [ 7.64e-10 1.32e-11 ]
%% Cell type:markdown id: tags:
## Plot
%% Cell type:code id: tags:
``` python
import matplotlib.pyplot as plt
from dune.fem.plotting import plotPointData
fig, axs = plt.subplots(1, 2, figsize=(8, 3))
axs[0].set_title("Solution")
plotPointData(ph, figure=(fig, axs[0]))
plotPointData(ph_g, figure=(fig, axs[0]), linewidth=0.04, colorbar=None)
phexact = space.interpolate(pexact, name="pexact")
ph_gexact = space_g.interpolate(p_gexact, name="p_gexact")
axs[1].set_title("Exact")
plotPointData(phexact, figure=(fig, axs[1]))
plotPointData(ph_gexact, figure=(fig, axs[1]), linewidth=0.04, colorbar=None)
```
%%%% Output: display_data
[Hidden Image Output]
......
......@@ -29,7 +29,8 @@
"\\begin{align}\n",
"p &= p_\\Gamma, \\\\\n",
"\\left( \\sigma(\\mathbf{u}) - \\alpha p I \\right) \\cdot \\mathbf{n} &= -p_\\Gamma \\mathbf{n}.\n",
"\\end{align}"
"\\end{align}\n",
"We consider a [T-junction](grids/tjunction.rst)."
]
},
{
......
%% Cell type:markdown id: tags:
# Poroelasticity
%% Cell type:markdown id: tags:
We consider quasi-static linear Biot-equations.
Find $\mathbf{u}, p$ s.t
\begin{align}
- \operatorname{div}( \mathbf{K} \nabla p ) &= 0, \\
- \operatorname{div} \left( \sigma( \mathbf{u} ) - \alpha p I \right) &= 0,
\end{align}
where
\begin{align}
\sigma( \mathbf{u} ) &= \lambda \operatorname{tr} \epsilon ( \mathbf{u} ) I + 2 \mu \epsilon( \mathbf{u} ), \\
\epsilon( \mathbf{u} ) &= \frac{1}{2} \left( \nabla \mathbf{u} + (\nabla \mathbf{u})^T \right) \\
\end{align}
and $\lambda, \mu > 0$ are the Lamé constants, $\alpha$ the Biot-Willis constant and $\mathbf{K}$ the permeability.
We include an interface $\Gamma$ with a fixed pressure $p_\Gamma$ and impose interior boundary conditions
\begin{align}
p &= p_\Gamma, \\
\left( \sigma(\mathbf{u}) - \alpha p I \right) \cdot \mathbf{n} &= -p_\Gamma \mathbf{n}.
\end{align}
We consider a [T-junction](grids/tjunction.rst).
%% Cell type:code id: tags:
``` python
from dune.grid import reader
from dune.mmesh import mmesh
file = "grids/tjunction.msh"
gridView = mmesh((reader.gmsh, file), 2)
igridView = gridView.hierarchicalGrid.interfaceGrid
```
%% Cell type:code id: tags:
``` python
from ufl import *
from dune.ufl import Constant
lamb = Constant( 0.5, name="lambda")
mu = Constant( 1, name="mu")
alpha = Constant( 1, name="alpha")
p_gamma = Constant( 1e2, name="p_gamma")
K = as_matrix([[1e-3, 0], [0, 1e-6]])
epsilon = lambda u: 0.5 * (nabla_grad(u) + nabla_grad(u).T)
sigma = lambda u: lamb * div(u) * Identity(2) + 2 * mu * epsilon(u)
```
%% Cell type:code id: tags:
``` python
from dune.fem.space import dglagrange
space = dglagrange(gridView, dimRange=3, order=1)
trial = TrialFunction(space)
test = TestFunction(space)
p, ux, uy = split(trial)
pp, uux, uuy = split(test)
u = as_vector([ux, uy ])
uu = as_vector([uux, uuy])
x = SpatialCoordinate(space)
n = FacetNormal(space)
h = FacetArea(space)
```
%% Cell type:markdown id: tags:
## Interior Penalty Discontinuous Galerkin (IPDG) Scheme
%% Cell type:code id: tags:
``` python
from dune.mmesh import skeleton, trace, interfaceIndicator
I = interfaceIndicator(igridView)
beta = Constant(1e3, name="beta")
# Pressure
a = inner(K * grad(p), grad(pp)) * dx
a += beta / h * inner(jump(p), jump(pp)) * (1-I)*dS
a -= inner(avg(K*grad(p )), n('+')) * jump(pp) * (1-I)*dS
a -= inner(avg(K*grad(pp)), n('+')) * jump(p ) * (1-I)*dS
a += beta / h * (p - 0) * pp * ds
a -= inner(K*grad(p ), n) * pp * ds
a -= inner(K*grad(pp), n) * (p - 0) * ds
# Pressure is continuous at the interface
a += beta / h * (p('+') - p_gamma) * pp('+') * I*dS
a -= inner(K*grad( p('+')), n('+')) * pp('+') * I*dS
a -= inner(K*grad(pp('+')), n('+')) * (p('+') - p_gamma) * I*dS
a += beta / h * (p('-') - p_gamma) * pp('-') * I*dS
a -= inner(K*grad( p('-')), n('-')) * pp('-') * I*dS
a -= inner(K*grad(pp('-')), n('-')) * (p('-') - p_gamma) * I*dS
# Displacement
a += inner(sigma(u), epsilon(uu)) * dx
a += alpha * inner(grad(p), uu) * dx
a += beta / h * inner(jump(u), jump(uu)) * (1-I)*dS
a -= dot(dot(avg(sigma(u )), n('+')), jump(uu)) * (1-I)*dS
a -= dot(dot(avg(sigma(uu)), n('+')), jump(u )) * (1-I)*dS
a += beta / h * inner(u - as_vector([0,0]), uu) * ds
a -= dot(dot(sigma(u ), n), uu) * ds
a -= dot(dot(sigma(uu), n), u - as_vector([0,0])) * ds
# Normal stress is -p at the interface
a -= -p_gamma * inner(uu('+'), n('+')) * I*dS
a -= -p_gamma * inner(uu('-'), n('-')) * I*dS
```
%% Cell type:code id: tags:
``` python
from dune.fem.scheme import galerkin
scheme = galerkin([a == 0], solver=("suitesparse", "umfpack"))
solution = space.interpolate([0,0,0], name="solution")
res = scheme.solve(target=solution)
```
%% Cell type:code id: tags:
``` python
import matplotlib.pyplot as plt
from dune.fem.plotting import plotPointData as plot
fig, axs = plt.subplots(1, 4, figsize=(20,5))
axs[0].set_title('Pressure')
plot(solution[0], figure=(fig, axs[0]), gridLines=None)
axs[1].set_title('Displacement (x)')
plot(solution[1], figure=(fig, axs[1]), gridLines=None)
axs[2].set_title('Displacement (y)')
plot(solution[2], figure=(fig, axs[2]), gridLines=None)
axs[2].set_title('Displacement')
plot(solution, vectors=[1,2], figure=(fig, axs[3]), gridLines=None)
```
%%%% Output: display_data
[Hidden Image Output]
......
......@@ -5,7 +5,7 @@
"id": "special-yorkshire",
"metadata": {},
"source": [
"# Two-Phase Navier-Stokes with Surface Tension"
"# Two-Phase Navier-Stokes with surface tension"
]
},
{
......@@ -33,7 +33,7 @@
"\n",
"Here, $\\nu_i, i=1,2,$ are the kinematic viscositities and $\\rho_i, i=1,2,$ the densities of the two phases, $g$ is the gravity, $\\sigma$ is the surface tension and $\\kappa$ is the signed mean curvature of the interface.\n",
"\n",
"Now, we use Dune-MMesh to compute the result of a rising bubble."
"Now, we use Dune-MMesh to compute the result of a rising bubble. Let us consider the [circle](grids/circle.rst) grid."
]
},
{
......
%% Cell type:markdown id:special-yorkshire tags:
# Two-Phase Navier-Stokes with Surface Tension
# Two-Phase Navier-Stokes with surface tension
%% Cell type:markdown id:announced-attention tags:
Let's consider the incompressible Navier-Stokes equations with two immiscible phases and surface tension.
A domain $\Omega \subset \mathbb{R}^2$ is assumed to be separeted into two phases $\Omega_i(t), i=1,2,$ by a sharp interface $\Gamma(t)$.
Find $(u, p)$ and $\Gamma(t)$ s.t.
\begin{align}
\renewcommand{\jump}[1]{[\mskip-5mu[ #1 ]\mskip-5mu]}
\end{align}
\begin{align}
u_t + (u \cdot \nabla u) u - \nu_i \Delta u + \nabla p &= g, &\qquad &\text{in } \Omega_i(t), \quad &i&=1,2,\\
\operatorname{div} u &= 0, &\qquad &\text{in } \Omega_i(t), \quad &i&=1,2,\\
\jump{\rho (\nu n (\nabla u + (\nabla u)^T) n - p ) } &= \sigma \kappa, &\qquad &\text{on } \Gamma(t),\\
\jump{u} &= 0, &\qquad &\text{on } \Gamma(t),\\
\dot x &= u, &\qquad &x \in \Gamma(t),\\
u(0) &= u_0, &\qquad &\text{in } \Omega_i(0), \quad &i&=1,2,\\
\Gamma(0) &= \Gamma_0.
\end{align}
Here, $\nu_i, i=1,2,$ are the kinematic viscositities and $\rho_i, i=1,2,$ the densities of the two phases, $g$ is the gravity, $\sigma$ is the surface tension and $\kappa$ is the signed mean curvature of the interface.
Now, we use Dune-MMesh to compute the result of a rising bubble.
Now, we use Dune-MMesh to compute the result of a rising bubble. Let us consider the [circle](grids/circle.rst) grid.
%% Cell type:code id:offshore-thomson tags:
``` python
from dune.grid import reader
from dune.mmesh import mmesh, trace, skeleton, domainMarker
from dune.fem.view import adaptiveLeafGridView as adaptive
dim = 2
file = "grids/circle.msh"
gridView = adaptive( mmesh((reader.gmsh, file), dim) )
hgrid = gridView.hierarchicalGrid
igridView = adaptive( hgrid.interfaceGrid )
```
%% Cell type:code id:brown-marsh tags:
``` python
from dune.ufl import Constant
mu0 = Constant( 1e-2, name="mu0")
mu1 = Constant( 1e-1, name="mu1")
rho0 = Constant( 3, name="rho0")
rho1 = Constant( 1, name="rho1")
sigma = Constant( 3e-1, name="sigma")
g = Constant( 10, name="g")
dt = 0.001
T = 1.5
```
%% Cell type:code id:constitutional-professor tags:
``` python
from dune.mmesh import domainMarker
dm = domainMarker(gridView)
mu = (1-dm) * mu0 + dm * mu1
rho = (1-dm) * rho0 + dm * rho1
nu = mu / rho
```
%% Cell type:code id:pleasant-possession tags:
``` python
from ufl import *
from dune.fem.space import dglagrange, product
uspace = dglagrange(gridView, dimRange=2, order=2)
pspace = dglagrange(gridView, order=1)
u = TrialFunction(uspace)
uu = TestFunction(uspace)
p = TrialFunction(pspace)
pp = TestFunction(pspace)
uh = uspace.interpolate([0,0], name="velocity")
uh1 = uspace.interpolate([0,0], name="velocity1")
ph = pspace.interpolate(0, name="pressure")
```
%% Cell type:markdown id:damaged-perth tags:
## Curvature
The mean curvature of the interface times its normal can be computed by
\begin{align*}
\int_{\Gamma(t)} \kappa \cdot \phi + \nabla x \cdot \nabla \phi \ dS = 0, &\qquad &\text{in } \Gamma(t).\\
\end{align*}
%% Cell type:code id:returning-buying tags:
``` python
from dune.fem.space import lagrange
from dune.fem.scheme import galerkin
kspace = lagrange(igridView, dimRange=dim, order=1)
k = TrialFunction(kspace)
kk = TestFunction(kspace)
curvature = kspace.interpolate([0]*dim, name="curvature")
ix = SpatialCoordinate(kspace)
C = inner(k, kk) * dx
C -= inner(grad(ix), grad(kk)) * dx
kscheme = galerkin([C == 0])
res = kscheme.solve(curvature)
```
%% Cell type:markdown id:intelligent-lightweight tags:
## Navier-Stokes equations
As solution strategy, we use an operator split known as Chorin's method.
%% Cell type:code id:crude-scale tags:
``` python
from dune.mmesh import interfaceIndicator
from dune.fem.scheme import galerkin
from dune.ufl import Constant
x = SpatialCoordinate(pspace)
n = FacetNormal(pspace)
I = interfaceIndicator(igridView)
uBC = as_vector([0,0])
tau = Constant(dt, name="tau")
penu = Constant(1e3, name="penaltyu")
penp = Constant(1e3, name="penaltyp")
a1 = inner(u - uh, uu) / tau * dx
a1 += inner(grad(uh) * uh, uu) * dx
a1 += nu * inner(grad(u), grad(uu)) * dx
a1 += g * uu[1] * dx
a1 += penu * inner(jump(u), jump(uu)) * dS
a1 += penu * inner(u - uBC, uu) * ds
a2 = inner(grad(p), grad(pp)) * dx
a2 += div(uh1) * pp / tau * dx
a2 += penp * jump(p * rho) * jump(pp) * dS
a2 += 1e-6 * p * pp * dx
S = sigma * inner(avg(skeleton(curvature)), n('+'))
S += dot(jump(mu * (nabla_grad(uh1) + nabla_grad(uh1).T), n('+')), n('+'))
a2 -= penp * S * jump(pp) * I * dS
a3 = inner(u - uh1, uu) * dx
a3 += tau * inner(grad(ph), uu) * dx
a3 += penu * inner(jump(u), jump(uu)) * dS
a3 += penu * inner(u - uBC, uu) * ds
A1 = galerkin([a1 == 0], solver=('suitesparse', 'umfpack'))
A2 = galerkin([a2 == 0], solver=('suitesparse', 'umfpack'))
A3 = galerkin([a3 == 0], solver=('suitesparse', 'umfpack'))
```
%% Cell type:markdown id:removed-modification tags:
## Moving
We can obtain the interface movement from the trace of the bulk velocity $u$.
%% Cell type:code id:turned-definition tags:
``` python
import numpy as np
def getShifts(uh):
mapper = igridView.hierarchicalGrid.leafView.indexSet
shifts = np.zeros((igridView.size(dim-1), dim))
for e in igridView.elements:
for v in e.subEntities(igridView.dimension):
x = e.geometry.toLocal(v.geometry.center)
shifts[ mapper.index(v) ] += 0.5 * trace(uh)(e, x)
return shifts
```
%% Cell type:markdown id:facial-temple tags:
## Timeloop
%% Cell type:code id:former-allen tags:
``` python
from dune.fem import parameter, adapt
parameter.append( { "fem.adaptation.method": "callback" } )
from dune.fem.plotting import plotPointData as plot
import matplotlib.pyplot as plt
N = 4
i = 0
fig, axs = plt.subplots(1, N, figsize=(16,6))
uh.interpolate([0,0])
ph.interpolate(0)
step = 0
t = 0
while t < T+dt:
kscheme.solve(curvature)
A1.solve(target=uh1)
A2.solve(target=ph)
A3.solve(target=uh)
hgrid.markElements()
hgrid.ensureInterfaceMovement( 2*dt*getShifts(uh) )
adapt([uh, uh1, ph, dm])
adapt([curvature])
hgrid.moveInterface( dt*getShifts(uh) )
t += dt
gridView.writeVTK("twophase-"+str(step), pointdata=[ph, uh], celldata=[dm],
nonconforming=True, subsampling=1)
igridView.writeVTK("interface-"+str(step), pointdata=[curvature])
step += 1
if int(N * t/T) > i:
plot(dm, figure=(fig, axs[i]), gridLines=None)
plot(uh, figure=(fig, axs[i]), gridLines=None, vectors=[0,1])
i += 1
```
%%%% Output: display_data
[Hidden Image Output]
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment