Commit f10ff94b authored by Lambert Theisen's avatar Lambert Theisen 🔥

Add GLS stabilization (with new TO, 3 new tests)

- Add Galerking Least Squaes Stabilization (GLS) (easy for linear equation)
- Includes addition of div3d3 and gen3d1 in tensoroperations module
- Includes gls keyword in input.yml with 5 parameters
- 3 new test cases added (heat, stress, r13)
- Method works out nicely, especially for R13 case!
- Further studies needed (use better computer for finer meshes)
parent f87b42b4
Pipeline #320789 passed with stages
in 28 minutes and 31 seconds
......@@ -15,12 +15,19 @@ meshes:
# - fields: List of FEM parameters (shape, degree)
# - shape: Element shape, e.g. Lagrange
# - degree: Element degree, e.g. 2
# - stabilization: Must contain cip
# - stabilization: Must contain cip and gls
# - cip: Collection of Continous Interior Penalty (CIP) parameters
# - enable: Enable CIP stabilization
# - delta_theta: Stabilization of grad(T)*grad(T_test) over edge
# - delta_u: Stabilization of grad(u)*grad(u_test) over edge
# - delta_p: Stabilization of grad(p)*grad(p_test) over edge
# - gls: Collection of Garlerkin Least Squares (GLS) parameters
# - enable: Enable GLS stabilization
# - tau_energy: Stabilization with energy eq. residual
# - tau_heatflux: Stabilization with heatflux eq. residual
# - tau_mass: Stabilization with mass eq. residual
# - tau_momentum: Stabilization with momentum eq. residual
# - tau_stress: Stabilization with stress eq. residual
elements:
theta:
shape: Lagrange
......@@ -43,6 +50,13 @@ stabilization:
delta_theta: 1.0
delta_u: 1.0
delta_p: 0.1
gls:
enable: False
tau_energy: 0.001
tau_heatflux: 0.001
tau_mass: 0.01
tau_momentum: 0.01
tau_stress: 0.01
# Formulation Parameters
# ======================
......
......@@ -15,12 +15,19 @@ meshes:
# - fields: List of FEM parameters (shape, degree)
# - shape: Element shape, e.g. Lagrange
# - degree: Element degree, e.g. 2
# - stabilization: Must contain cip
# - stabilization: Must contain cip and gls
# - cip: Collection of Continous Interior Penalty (CIP) parameters
# - enable: Enable CIP stabilization
# - delta_theta: Stabilization of grad(T)*grad(T_test) over edge
# - delta_u: Stabilization of grad(u)*grad(u_test) over edge
# - delta_p: Stabilization of grad(p)*grad(p_test) over edge
# - gls: Collection of Garlerkin Least Squares (GLS) parameters
# - enable: Enable GLS stabilization
# - tau_energy: Stabilization with energy eq. residual
# - tau_heatflux: Stabilization with heatflux eq. residual
# - tau_mass: Stabilization with mass eq. residual
# - tau_momentum: Stabilization with momentum eq. residual
# - tau_stress: Stabilization with stress eq. residual
elements:
theta:
shape: Lagrange
......@@ -43,6 +50,13 @@ stabilization:
delta_theta: 1.0
delta_u: 1.0
delta_p: 0.1
gls:
enable: False
tau_energy: 0.001
tau_heatflux: 0.001
tau_mass: 0.01
tau_momentum: 0.01
tau_stress: 0.01
# Formulation Parameters
# ======================
......
......@@ -16,12 +16,19 @@ meshes:
# - fields: List of FEM parameters (shape, degree)
# - shape: Element shape, e.g. Lagrange
# - degree: Element degree, e.g. 2
# - stabilization: Must contain cip
# - stabilization: Must contain cip and gls
# - cip: Collection of Continous Interior Penalty (CIP) parameters
# - enable: Enable CIP stabilization
# - delta_theta: Stabilization of grad(T)*grad(T_test) over edge
# - delta_u: Stabilization of grad(u)*grad(u_test) over edge
# - delta_p: Stabilization of grad(p)*grad(p_test) over edge
# - gls: Collection of Garlerkin Least Squares (GLS) parameters
# - enable: Enable GLS stabilization
# - tau_energy: Stabilization with energy eq. residual
# - tau_heatflux: Stabilization with heatflux eq. residual
# - tau_mass: Stabilization with mass eq. residual
# - tau_momentum: Stabilization with momentum eq. residual
# - tau_stress: Stabilization with stress eq. residual
elements:
theta:
shape: Lagrange
......@@ -44,6 +51,13 @@ stabilization:
delta_theta: 1.0
delta_u: 1.0
delta_p: 0.1
gls:
enable: False
tau_energy: 0.001
tau_heatflux: 0.001
tau_mass: 0.01
tau_momentum: 0.01
tau_stress: 0.01
# Formulation Parameters
# ======================
......
......@@ -15,12 +15,19 @@ meshes:
# - fields: List of FEM parameters (shape, degree)
# - shape: Element shape, e.g. Lagrange
# - degree: Element degree, e.g. 2
# - stabilization: Must contain cip
# - stabilization: Must contain cip and gls
# - cip: Collection of Continous Interior Penalty (CIP) parameters
# - enable: Enable CIP stabilization
# - delta_theta: Stabilization of grad(T)*grad(T_test) over edge
# - delta_u: Stabilization of grad(u)*grad(u_test) over edge
# - delta_p: Stabilization of grad(p)*grad(p_test) over edge
# - gls: Collection of Garlerkin Least Squares (GLS) parameters
# - enable: Enable GLS stabilization
# - tau_energy: Stabilization with energy eq. residual
# - tau_heatflux: Stabilization with heatflux eq. residual
# - tau_mass: Stabilization with mass eq. residual
# - tau_momentum: Stabilization with momentum eq. residual
# - tau_stress: Stabilization with stress eq. residual
elements:
theta:
shape: Lagrange
......@@ -43,6 +50,13 @@ stabilization:
delta_theta: 1.0
delta_u: 1.0
delta_p: 0.1
gls:
enable: False
tau_energy: 0.001
tau_heatflux: 0.001
tau_mass: 0.01
tau_momentum: 0.01
tau_stress: 0.01
# Formulation Parameters
# ======================
......
......@@ -15,12 +15,19 @@ meshes:
# - fields: List of FEM parameters (shape, degree)
# - shape: Element shape, e.g. Lagrange
# - degree: Element degree, e.g. 2
# - stabilization: Must contain cip
# - stabilization: Must contain cip and gls
# - cip: Collection of Continous Interior Penalty (CIP) parameters
# - enable: Enable CIP stabilization
# - delta_theta: Stabilization of grad(T)*grad(T_test) over edge
# - delta_u: Stabilization of grad(u)*grad(u_test) over edge
# - delta_p: Stabilization of grad(p)*grad(p_test) over edge
# - gls: Collection of Garlerkin Least Squares (GLS) parameters
# - enable: Enable GLS stabilization
# - tau_energy: Stabilization with energy eq. residual
# - tau_heatflux: Stabilization with heatflux eq. residual
# - tau_mass: Stabilization with mass eq. residual
# - tau_momentum: Stabilization with momentum eq. residual
# - tau_stress: Stabilization with stress eq. residual
elements:
theta:
shape: Lagrange
......@@ -43,6 +50,13 @@ stabilization:
delta_theta: 1.0
delta_u: 1.0
delta_p: 0.1
gls:
enable: False
tau_energy: 0.001
tau_heatflux: 0.001
tau_mass: 0.01
tau_momentum: 0.01
tau_stress: 0.01
# Formulation Parameters
# ======================
......
......@@ -46,12 +46,19 @@ class Input:
# - fields: List of FEM parameters (shape, degree)
# - shape: Element shape, e.g. Lagrange
# - degree: Element degree, e.g. 2
# - stabilization: Must contain cip
# - stabilization: Must contain cip and gls
# - cip: Collection of Continous Interior Penalty (CIP) parameters
# - enable: Enable CIP stabilization
# - delta_theta: Stabilization of grad(T)*grad(T_test) over edge
# - delta_u: Stabilization of grad(u)*grad(u_test) over edge
# - delta_p: Stabilization of grad(p)*grad(p_test) over edge
# - gls: Collection of Garlerkin Least Squares (GLS) parameters
# - enable: Enable GLS stabilization
# - tau_energy: Stabilization with energy eq. residual
# - tau_heatflux: Stabilization with heatflux eq. residual
# - tau_mass: Stabilization with mass eq. residual
# - tau_momentum: Stabilization with momentum eq. residual
# - tau_stress: Stabilization with stress eq. residual
elements:
theta:
shape: Lagrange
......@@ -74,6 +81,13 @@ class Input:
delta_theta: 1.0
delta_u: 1.0
delta_p: 0.01
gls:
enable: False
tau_energy: 0.001
tau_heatflux: 0.001
tau_mass: 0.01
tau_momentum: 0.01
tau_stress: 0.01
# Formulation Parameters
# ======================
......@@ -341,27 +355,59 @@ class Input:
"stabilization": {
"type": "dict",
"required": True,
"keysrules": {"type": "string", "regex": "cip"},
"valuesrules": {
"type": "dict",
"schema": {
"enable": {
"type": "boolean",
"required": True
},
"delta_theta": {
"type": "float",
"required": True
},
"delta_u": {
"type": "float",
"required": True
},
"delta_p": {
"type": "float",
"required": True
},
}
"schema": {
"cip": {
"type": "dict",
"required": True,
"schema": {
"enable": {
"type": "boolean",
"required": True
},
"delta_theta": {
"type": "float",
"required": True
},
"delta_u": {
"type": "float",
"required": True
},
"delta_p": {
"type": "float",
"required": True
},
}
},
"gls": {
"type": "dict",
"required": True,
"schema": {
"enable": {
"type": "boolean",
"required": True
},
"tau_energy": {
"type": "float",
"required": True
},
"tau_heatflux": {
"type": "float",
"required": True
},
"tau_mass": {
"type": "float",
"required": True
},
"tau_momentum": {
"type": "float",
"required": True
},
"tau_stress": {
"type": "float",
"required": True
},
}
},
}
},
"elements": {
......
......@@ -69,11 +69,21 @@ class Solver:
self.cell = self.mesh.ufl_cell()
self.time = time
self.mode = params["mode"]
# CIP
self.use_cip = self.params["stabilization"]["cip"]["enable"]
self.delta_theta = self.params["stabilization"]["cip"]["delta_theta"]
self.delta_u = self.params["stabilization"]["cip"]["delta_u"]
self.delta_p = self.params["stabilization"]["cip"]["delta_p"]
# GLS
self.use_gls = self.params["stabilization"]["gls"]["enable"]
self.tau_energy = self.params["stabilization"]["gls"]["tau_energy"]
self.tau_heatflux = self.params["stabilization"]["gls"]["tau_heatflux"]
self.tau_mass = self.params["stabilization"]["gls"]["tau_mass"]
self.tau_momentum = self.params["stabilization"]["gls"]["tau_momentum"]
self.tau_stress = self.params["stabilization"]["gls"]["tau_stress"]
self.write_pdfs = self.params["postprocessing"]["write_pdfs"]
self.write_vecs = self.params["postprocessing"]["write_vecs"]
self.massflow = self.params["postprocessing"]["massflow"]
......@@ -509,6 +519,11 @@ class Solver:
delta_theta = df.Constant(self.delta_theta)
delta_u = df.Constant(self.delta_u)
delta_p = df.Constant(self.delta_p)
tau_energy = df.Constant(self.tau_energy)
tau_heatflux = df.Constant(self.tau_heatflux)
tau_mass = df.Constant(self.tau_mass)
tau_momentum = df.Constant(self.tau_momentum)
tau_stress = df.Constant(self.tau_stress)
# Define custom measeasures for boundary edges and inner edges
df.dx = df.Measure("dx", domain=mesh, subdomain_data=regions)
......@@ -548,11 +563,15 @@ class Solver:
else:
cpl = 0
# Stabilization switch
# Stabilization
if self.use_cip:
cip = 1
else:
cip = 0
if self.use_gls:
gls = 1
else:
gls = 0
# Setup normal/tangential projections
# => tangential (tx,ty) = (-ny,nx) = perp(n) only for 2D
......@@ -650,7 +669,7 @@ class Solver:
df.inner(v, df.grad(p))
) * df.dx(reg) for reg in regs.keys()])
# 3) CIP Stabilization:
# 3.1) CIP Stabilization:
def j_theta(theta, kappa):
return (
+ delta_theta * h_avg**3 *
......@@ -669,6 +688,56 @@ class Solver:
df.jump(df.grad(p), n_vec) * df.jump(df.grad(q), n_vec)
) * df.dS
# 3.2) GLS Stabilization
def gls_heat(theta, kappa, s, r):
return sum([(
tau_energy * h_msh**1 * (
df.inner(
df.div(s) + cpl * df.div(u) - f_heat,
df.div(r) + cpl * df.div(v)
)
) # energy
+ tau_heatflux * h_msh**1 *
df.inner(
(5 / 2) * df.grad(theta)
- (12 / 5) * regs[reg]["kn"] * df.div(to.stf3d2(df.grad(s)))
- (1 / 6) * regs[reg]["kn"] * 12 * df.grad(df.div(s))
+ 1 / regs[reg]["kn"] * 2 / 3 * s,
(5 / 2) * df.grad(kappa)
- (12 / 5) * regs[reg]["kn"] * df.div(to.stf3d2(df.grad(r)))
- (1 / 6) * regs[reg]["kn"] * 12 * df.grad(df.div(r))
+ 1 / regs[reg]["kn"] * 2 / 3 * r
) # heatflux
) * df.dx(reg) for reg in regs.keys()])
def gls_stress(p, q, u, v, sigma, psi):
return sum([(
tau_mass * h_msh**1.5 *
df.inner(
df.div(v), df.div(u) - f_mass
) # mass
+ tau_momentum * h_msh**1.5 *
df.inner(
df.grad(q) + df.div(psi),
df.grad(p) + df.div(sigma) - f_body
) # momentum
+ tau_stress * h_msh**1.5 *
df.inner(
cpl * (4 / 5) * to.gen3dTF2(df.grad(r))
+ 2 * to.stf3d2(to.gen3d2(df.grad(v)))
- 2 * regs[reg]["kn"] * to.div3d3(
to.stf3d3(to.grad3dOf2(to.gen3dTF2(psi)))
)
+ (1 / regs[reg]["kn"]) * to.gen3dTF2(psi),
cpl * (4 / 5) * to.gen3dTF2(df.grad(s))
+ 2 * to.stf3d2(to.gen3d2(df.grad(u)))
- 2 * regs[reg]["kn"] * to.div3d3(
to.stf3d3(to.grad3dOf2(to.gen3dTF2(sigma)))
)
+ (1 / regs[reg]["kn"]) * to.gen3dTF2(sigma)
) # stress
) * df.dx(reg) for reg in regs.keys()])
# Setup all equations
A = [None] * 5
L = [None] * 5
......@@ -701,20 +770,37 @@ class Solver:
) * q
) * df.ds(bc) for bc in bcs.keys()])
# Combine all equations to compound weak form and add CIP
# Combine all equations to compound weak form and add stabilization
if self.mode == "heat":
self.form_lhs = sum(A[0:2]) + cip * j_theta(theta, kappa)
self.form_rhs = sum(L[0:2])
self.form_lhs = sum(A[0:2]) + (
cip * (j_theta(theta, kappa))
+ gls * df.lhs(gls_heat(theta, kappa, s, r))
)
self.form_rhs = sum(L[0:2]) + (
gls * df.rhs(gls_heat(theta, kappa, s, r))
)
elif self.mode == "stress":
self.form_lhs = sum(A[2:5]) + cip * (
j_u(u, v) + j_p(p, q)
self.form_lhs = sum(A[2:5]) + (
cip * (j_u(u, v) + j_p(p, q))
+ gls * df.lhs(gls_stress(p, q, u, v, sigma, psi))
)
self.form_rhs = sum(L[2:5]) + (
gls * df.rhs(gls_stress(p, q, u, v, sigma, psi))
)
self.form_rhs = sum(L[2:5])
elif self.mode == "r13":
self.form_lhs = sum(A) + cip * (
j_theta(theta, kappa) + j_u(u, v) + j_p(p, q)
self.form_lhs = sum(A) + (
cip * (j_theta(theta, kappa) + j_u(u, v) + j_p(p, q))
+ gls * df.lhs(
gls_heat(theta, kappa, s, r)
+ gls_stress(p, q, u, v, sigma, psi)
)
)
self.form_rhs = sum(L) + (
gls * df.rhs(
gls_heat(theta, kappa, s, r)
+ gls_stress(p, q, u, v, sigma, psi)
)
)
self.form_rhs = sum(L)
def solve(self):
"""
......
......@@ -23,9 +23,14 @@ import ufl
def stf3d2(rank2_2d):
r"""
Return the synthetic 3D symmetric and trace-free part of a 2D 2-tensor.
Return the 3D symmetric and trace-free part of a 2D 2-tensor.
Return the synthetic 3D symmetric and trace-free (dev(sym(.)))
.. warning::
Return a 2-tensor with the same dimensions as the input tensor.
For the :math:`2 \times 2` case, return the 3D symmetric and
trace-free (dev(sym(.)))
:math:`B \in \mathbb{R}^{2 \times 2}`
of the 2D 2-tensor
:math:`A \in \mathbb{R}^{2 \times 2}`.
......@@ -39,8 +44,10 @@ def stf3d2(rank2_2d):
B &= (A)_\mathrm{dev} = \frac{1}{2} (A)_\mathrm{sym}
- \frac{1}{3} \mathrm{tr}(A) I_{2 \times 2}
"""
dim = len(rank2_2d[:, 0])
symm = 1 / 2 * (rank2_2d + ufl.transpose(rank2_2d))
return symm - (1 / 3) * ufl.tr(symm) * ufl.Identity(2)
return symm - (1 / 3) * ufl.tr(symm) * ufl.Identity(dim)
def sym3d3(rank3_3d):
......@@ -123,6 +130,27 @@ def stf3d3(rank3_3d):
return ufl.as_tensor(tracefree_ijk, (i, j, k))
def div3d3(rank3_3d):
r"""
Return the 3D divergence of a 3-tensor.
Return the 3D divergence of a 3-tensor as
.. math::
{(\mathrm{div}(m))}_{ij} = \frac{\partial m_{ijk}}{\partial x_k}
Compare with [SCH2009]_ (p. 92).
.. [SCH2009] Heinz Schade, Klaus Neemann (2009). “Tensoranalysis”.
2. überarbeitete Auflage.
"""
i, j, k = ufl.indices(3)
div_ij = rank3_3d[i, j, 0].dx(0) + rank3_3d[i, j, 1].dx(1)
return ufl.as_tensor(div_ij, (i, j))
def gen3dTF2(rank2_2d):
r"""
Generate a 3D tracefree 2-tensor from a 2D 2-tensor.
......@@ -170,6 +198,19 @@ def gen3dTF2(rank2_2d):
])
def gen3d1(rank1_2d):
r"""
Return synthetic 3d version of 2d vector with zero last component.
Return synthetic 3d version of 2d vector
:math:`s_{\mathrm{in}} = (s_x, s_y)`
as
:math:`s_{\mathrm{out}} = (s_x, s_y, 0)`
.
"""
return df.as_vector([rank1_2d[0], rank1_2d[1], 0])
def gen3d2(rank2_2d):
r"""
Generate a 3D 2-tensor from a 2D 2-tensor (add zeros to last dimensions).
......
# General
# =======
# - output_folder: Used as output folder
output_folder: heat_01_coeffs_p1p1_gls
# Meshes
# ======
# - meshes: List of input meshes in h5 format to run simulations on
meshes:
- ../mesh/ring0.h5
- ../mesh/ring1.h5
- ../mesh/ring2.h5
- ../mesh/ring3.h5
- ../mesh/ring4.h5
# - ../mesh/ring5.h5
# - ../mesh/ring6.h5
# - ../mesh/ring7.h5
# Numerical Parameters
# ====================
# - elements: Must contain the fields: theta, s, p, u, sigma
# - fields: List of FEM parameters (shape, degree)
# - shape: Element shape, e.g. Lagrange
# - degree: Element degree, e.g. 2
# - stabilization: Must contain cip and gls
# - cip: Collection of Continous Interior Penalty (CIP) parameters
# - enable: Enable CIP stabilization
# - delta_theta: Stabilization of grad(T)*grad(T_test) over edge
# - delta_u: Stabilization of grad(u)*grad(u_test) over edge
# - delta_p: Stabilization of grad(p)*grad(p_test) over edge
# - gls: Collection of Garlerkin Least Squares (GLS) parameters
# - enable: Enable GLS stabilization
# - tau_energy: Stabilization with energy eq. residual
# - tau_heatflux: Stabilization with heatflux eq. residual
# - tau_mass: Stabilization with mass eq. residual
# - tau_momentum: Stabilization with momentum eq. residual
# - tau_stress: Stabilization with stress eq. residual
elements:
theta:
shape: Lagrange
degree: 1
s:
shape: Lagrange
degree: 1
p:
shape: Lagrange
degree: 1
u:
shape: Lagrange
degree: 1
sigma:
shape: Lagrange
degree: 1
stabilization:
cip:
enable: False
delta_theta: 1.0
delta_u: 1.0
delta_p: 0.1
gls:
enable: True
tau_energy: 0.001
tau_heatflux: 0.001
tau_mass: 0.01
tau_momentum: 0.01
tau_stress: 0.01
# Formulation Parameters
# ======================
# - nsd: Number of spatial dimensions == 2
# - mode: Formulation mode, one of heat, stress, r13
# - heat_source: Heat source function for mode==heat||r13
# - mass_source: Mass source function for mode==stress||r13