27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182 | def viscosity_term(eq: Equation, trial: Argument | Indexed | Function) -> Form:
r"""Viscosity term $-nabla * (mu nabla u)$ in the momentum equation.
Using the symmetric interior penalty method (Epshteyn & Rivière, 2007), the weak
form becomes
$$
{:( -int_Omega nabla * (mu grad u) phi dx , = , int_Omega mu (grad phi) * (grad u) dx ),
( , - , int_(cc"I" uu cc"I"_v) "jump"(phi bb n) * "avg"(mu grad u) dS
- int_(cc"I" uu cc"I"_v) "jump"(u bb n) * "avg"(mu grad phi) dS ),
( , + , int_(cc"I" uu cc"I"_v) sigma "avg"(mu) "jump"(u bb n) * "jump"(phi bb n) dS )
:}
$$
where σ is a penalty parameter.
Epshteyn, Y., & Rivière, B. (2007).
Estimation of penalty parameters for symmetric interior penalty Galerkin methods.
Journal of Computational and Applied Mathematics, 206(2), 843-872.
"""
mu = eq.approximation.mu
stress = eq.stress
F = inner(nabla_grad(eq.test), stress) * eq.dx
# A solution-dependent `mu` adds the penalty derivative below. Test each
# coefficient behind `trial` to include every mixed field, such as pressure.
# Spatial variation alone does not make `mu` solution-dependent.
mu_nonlinear = any(depends_on(mu, c) for c in extract_coefficients(trial))
sigma = interior_penalty_factor(eq)
sigma *= FacetArea(eq.mesh) / avg(CellVolume(eq.mesh))
if not is_continuous(eq.trial_space):
if mu_nonlinear:
raise NotImplementedError(
"Symmetric SIPG interior-facet (dS) terms for a solution-dependent "
"viscosity are not implemented for discontinuous velocity elements."
)
trial_tensor_jump = eq.approximation.deviatoric_tensor_from_grad(
tensor_jump(eq.n, trial)
)
F += (
sigma
* inner(tensor_jump(eq.n, eq.test), avg(mu) * trial_tensor_jump)
* eq.dS
)
F -= inner(avg(mu * nabla_grad(eq.test)), trial_tensor_jump) * eq.dS
F -= inner(tensor_jump(eq.n, eq.test), avg(stress)) * eq.dS
# The symmetrising term of a weak velocity boundary condition is the
# transpose of the flux term, which makes it the derivative of the stress in
# the direction of the test function, $D\sigma(u)[\phi]$. It is taken from
# the stress expression this equation carries rather than rebuilt from the
# approximation, because only that expression holds the state a particular
# formulation puts into the stress: the internal variables of the
# viscoelastic solvers, or the stress carried over from the previous step.
# Differentiating it picks up the right shear coefficient in each case, and
# for a solution-dependent viscosity it also picks up the $D\mu[\phi]$
# contribution that makes the boundary Jacobian symmetric.
weak_velocity_bcs = any(bc.keys() & {"u", "un"} for bc in eq.bcs.values())
if weak_velocity_bcs:
if not any(depends_on(stress, c) for c in extract_coefficients(trial)):
raise ValueError(
"The stress supplied to viscosity_term does not depend on the "
"trial function, so the symmetrising term of the weak velocity "
"boundary conditions would be identically zero and their "
"Jacobian would not be symmetric. Build the stress from the "
"same function the residual is evaluated at."
)
tangent_stress = expand_derivatives(derivative(stress, trial, eq.test))
# NOTE: Unspecified boundaries result in free stress (i.e. free in all directions).
# NOTE: "un" can be combined with "stress" provided the stress component is
# tangential (e.g. no normal flow with wind)
for bc_id, bc in eq.bcs.items():
if "u" in bc and any(bc_type in bc for bc_type in ["stress", "un"]):
raise ValueError(
'"stress" or "un" cannot be specified if "u" is already given.'
)
if "normal_stress" in bc and any(bc_type in bc for bc_type in ["u", "un"]):
raise ValueError(
'"u" or "un" cannot be specified if "normal_stress" is already given.'
)
if "u" in bc:
w = trial - bc["u"]
jump_gradient = outer(eq.n, w)
# Penalty term, similar to the above term for the DG dS integrals.
# The approximation converts the boundary jump tensor into a stress
# the same way it converts a velocity gradient, so the penalty
# carries every part of the stress the boundary condition
# constrains, volumetric part included.
F += (
2
* sigma
* inner(
outer(eq.n, eq.test),
eq.approximation.stress_from_grad(jump_gradient),
)
* eq.ds(bc_id)
)
# Symmetrising term, the transpose of the flux integration by parts.
F -= dot(w, dot(tangent_stress, eq.n)) * eq.ds(bc_id)
F -= inner(outer(eq.n, eq.test), stress) * eq.ds(bc_id)
# Derivative of the penalty term through mu: the penalty functional
# is $\sigma_{pen}\,\langle G, S(G)\rangle$ with
# $G = n \otimes w$ (jump_gradient) and $S$ the stress from a
# gradient-like tensor, so its exact first variation (the residual)
# picks up this extra term whenever $\mu$ itself depends on the
# trial. Only the deviatoric part of $S$ scales with $\mu$, so only
# that part appears here. This makes the resulting Newton Jacobian
# (the second variation) symmetric by construction.
if mu_nonlinear:
dmu = expand_derivatives(derivative(mu, trial, eq.test))
jump_tensor = eq.approximation.deviatoric_tensor_from_grad(jump_gradient)
F += sigma * dmu * inner(jump_gradient, jump_tensor) * eq.ds(bc_id)
if "un" in bc:
un_jump = dot(eq.n, trial) - bc["un"]
w = un_jump * eq.n
jump_gradient = outer(eq.n, w)
# Penalty term, as in the "u" branch but with the jump restricted to
# its normal component. The trace of the jump tensor is the normal
# jump itself, so an approximation with a bulk modulus penalises the
# normal jump volumetrically as well as deviatorically.
F += (
2
* sigma
* inner(
outer(eq.n, eq.test),
eq.approximation.stress_from_grad(jump_gradient),
)
* eq.ds(bc_id)
)
# Symmetrising term, as in the "u" branch but with the jump restricted
# to its normal component (free-slip/free-stress tangential direction).
F -= dot(w, dot(tangent_stress, eq.n)) * eq.ds(bc_id)
# We only keep the normal part of stress; the tangential part is assumed to
# be zero stress (i.e. free slip) or prescribed via "stress".
F -= dot(eq.n, eq.test) * dot(eq.n, dot(stress, eq.n)) * eq.ds(bc_id)
# Derivative of the penalty term through mu, as in the "u" branch
# above, restricted to the normal component of the jump.
if mu_nonlinear:
dmu = expand_derivatives(derivative(mu, trial, eq.test))
jump_tensor = eq.approximation.deviatoric_tensor_from_grad(jump_gradient)
F += sigma * dmu * inner(jump_gradient, jump_tensor) * eq.ds(bc_id)
if "stress" in bc: # a momentum flux, a.k.a. "force"
# Here we need only the third term because we assume jump_u = 0
# (u_ext = trial) and stress = n . (mu . stress_tensor).
F -= dot(eq.test, bc["stress"]) * eq.ds(bc_id)
if "normal_stress" in bc:
F += dot(eq.test, bc["normal_stress"] * eq.n) * eq.ds(bc_id)
return F
|