-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathcg_dg_advec_diffusion.py
261 lines (220 loc) · 9.05 KB
/
cg_dg_advec_diffusion.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
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
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
# Scheme based on "A finite element method for domain decomposition
# with non-matching grids" by Becker et al. but with a DG scheme to
# solve the advection diffusion equation on half of the domain, and
# a standard CG Poisson solver on the other half.
# Consider a square domain on which we wish to solve the
# advection diffusion equations. The velocity field is given by
# (0.5 - x_1, 0.0) in the bottom half of the domain, and (0.0, 0.0)
# in the top half. We solve the bottom half of the domain with a
# DG advection-diffusion solver, and the top half with a standard
# CG solver. We enforce the Dirichlet boundary condition weakly
# for the DG scheme and strongly for the CG scheme. The assumed
# solution is u = sin(\pi * x_0) * sin(\pi * x_1). In this problem,
# the bottom half can be thought of as a fluid and the top half
# a solid, and the unknown u is the temperature field.
# NOTE: Since the velocity goes to zero at the interface x[1] = 0.5,
# the coupling is due only to the diffusion. No advective interface
# terms have been added
from dolfinx import mesh, fem, io
from mpi4py import MPI
import ufl
from ufl import inner, grad, dot, avg, div, jump
import numpy as np
from petsc4py import PETSc
from utils import (
norm_L2,
convert_facet_tags,
interface_int_entities,
interior_facet_int_entities,
)
from dolfinx.fem.petsc import (
assemble_matrix_block,
create_vector_block,
assemble_vector_block,
)
from meshing import create_divided_square
from poisson_domain_decomp import jump_i, grad_avg_i
def u_e(x, module=np):
"Function to represent the exact solution"
# return module.exp(- ((x[0] - 0.5)**2 + (x[1] - 0.5)**2) / (2 * 0.15**2))
return module.sin(module.pi * x[0]) * module.sin(module.pi * x[1])
# Set some parameters
num_time_steps = 10
k_0 = 3 # Polynomial degree in omega_0
k_1 = 3 # Polynomial degree in omgea_1
delta_t = 1 # TODO Make constant
h = 0.05 # Maximum cell diameter
# Create the mesh
comm = MPI.COMM_WORLD
msh, ct, ft, vol_ids, bound_ids = create_divided_square(comm, h)
# Create sub-meshes of omega_0 and omega_1 so that we can create
# different function spaces over each part of the domain
tdim = msh.topology.dim
submesh_0, sm_0_to_msh = mesh.create_submesh(msh, tdim, ct.find(vol_ids["omega_0"]))[:2]
submesh_1, sm_1_to_msh = mesh.create_submesh(msh, tdim, ct.find(vol_ids["omega_1"]))[:2]
# Define function spaces on each submesh
V_0 = fem.functionspace(submesh_0, ("Discontinuous Lagrange", k_0))
V_1 = fem.functionspace(submesh_1, ("Lagrange", k_1))
W = ufl.MixedFunctionSpace(V_0, V_1)
# Test and trial functions
u = ufl.TrialFunctions(W)
v = ufl.TestFunctions(W)
# We use msh as the integration domain, so we require maps from
# cells in msh to cells in submesh_0 and submesh_1
cell_imap = msh.topology.index_map(tdim)
num_cells = cell_imap.size_local + cell_imap.num_ghosts
msh_to_sm_0 = np.full(num_cells, -1)
msh_to_sm_0[sm_0_to_msh] = np.arange(len(sm_0_to_msh))
msh_to_sm_1 = np.full(num_cells, -1)
msh_to_sm_1[sm_1_to_msh] = np.arange(len(sm_1_to_msh))
# Create integration entities for the interface integral
interface_facets = ft.find(bound_ids["interface"])
domain_0_cells = ct.find(vol_ids["omega_0"])
domain_1_cells = ct.find(vol_ids["omega_1"])
interface_entities, msh_to_sm_0, msh_to_sm_1 = interface_int_entities(
msh, interface_facets, msh_to_sm_0, msh_to_sm_1
)
# Compute integration entities for boundary terms
boundary_entities = [
(
bound_ids["boundary_0"],
fem.compute_integration_domains(
fem.IntegralType.exterior_facet,
msh.topology,
ft.find(bound_ids["boundary_0"]),
),
)
]
# Compute integration entities for the interior facet integrals
# over omega_0. These are needed for the DG scheme
omega_0_int_entities = interior_facet_int_entities(submesh_0, sm_0_to_msh)
# Create measures
dx = ufl.Measure("dx", domain=msh, subdomain_data=ct)
ds = ufl.Measure("ds", domain=msh, subdomain_data=boundary_entities)
dS = ufl.Measure(
"dS",
domain=msh,
subdomain_data=[
(bound_ids["interface"], interface_entities),
(bound_ids["omega_0_int_facets"], omega_0_int_entities),
],
)
x = ufl.SpatialCoordinate(msh)
kappa = [1.0 + 0.1 * ufl.sin(ufl.pi * x[0]) * ufl.sin(ufl.pi * x[1]) for _ in range(2)]
# Define forms
# TODO Add k dependency
gamma_int = 10 * 2 * kappa[0] * kappa[1] / (kappa[0] + kappa[1]) # Penalty param on interface
gamma_dg = 10 * k_0**2 # Penalty parm for DG method
h = ufl.CellDiameter(msh)
n = ufl.FacetNormal(msh)
u_0_n = fem.Function(V_0)
u_1_n = fem.Function(V_1)
w = ufl.as_vector((0.5 - x[1], 0.0))
lmbda = ufl.conditional(ufl.gt(dot(w, n), 0), 1, 0)
alpha = [1.0, 1.0]
# DG scheme in Omega_0
a = (
inner(alpha[0] * u[0] / delta_t, v[0]) * dx(vol_ids["omega_0"])
- alpha[0] * inner(w * u[0], grad(v[0])) * dx(vol_ids["omega_0"])
+ alpha[0]
* inner(
lmbda("+") * dot(w("+"), n("+")) * u[0]("+") - lmbda("-") * dot(w("-"), n("-")) * u[0]("-"),
jump(v[0]),
)
* dS(bound_ids["omega_0_int_facets"])
+ alpha[0] * inner(lmbda * dot(w, n) * u[0], v[0]) * ds(bound_ids["boundary_0"])
+ inner(kappa[0] * grad(u[0]), grad(v[0])) * dx(vol_ids["omega_0"])
- inner(kappa[0] * avg(grad(u[0])), jump(v[0], n)) * dS(bound_ids["omega_0_int_facets"])
- inner(kappa[0] * jump(u[0], n), avg(grad(v[0]))) * dS(bound_ids["omega_0_int_facets"])
+ (gamma_dg / avg(h))
* inner(kappa[0] * jump(u[0], n), jump(v[0], n))
* dS(bound_ids["omega_0_int_facets"])
- inner(kappa[0] * grad(u[0]), v[0] * n) * ds(bound_ids["boundary_0"])
- inner(kappa[0] * grad(v[0]), u[0] * n) * ds(bound_ids["boundary_0"])
+ (gamma_dg / h) * inner(kappa[0] * u[0], v[0]) * ds(bound_ids["boundary_0"])
)
# CG scheme in Omega_1
a += inner(alpha[1] * u[1] / delta_t, v[1]) * dx(vol_ids["omega_1"]) + inner(
kappa[1] * grad(u[1]), grad(v[1])
) * dx(vol_ids["omega_1"])
a += (
-inner(grad_avg_i(u, kappa), jump_i(v, n)) * dS(bound_ids["interface"])
- inner(jump_i(u, n), grad_avg_i(v, kappa)) * dS(bound_ids["interface"])
+ gamma_int / avg(h) * inner(jump_i(u, n), jump_i(v, n)) * dS(bound_ids["interface"])
)
# Compile LHS forms
entity_maps = {submesh_0: msh_to_sm_0, submesh_1: msh_to_sm_1}
a = fem.form(ufl.extract_blocks(a), entity_maps=entity_maps)
# Forms for the right-hand side
# TODO Add time derivative for unsteady problems
f_0 = alpha[0] * dot(w, grad(u_e(ufl.SpatialCoordinate(msh), module=ufl))) - div(
kappa[0] * grad(u_e(ufl.SpatialCoordinate(msh), module=ufl))
)
f_1 = -div(kappa[1] * grad(u_e(ufl.SpatialCoordinate(msh), module=ufl)))
u_D = fem.Function(V_0)
u_D.interpolate(u_e)
L = (
inner(f_0, v[0]) * dx(vol_ids["omega_0"])
- alpha[0] * inner((1 - lmbda) * dot(w, n) * u_D, v[0]) * ds(bound_ids["boundary_0"])
+ inner(alpha[0] * u_0_n / delta_t, v[0]) * dx(vol_ids["omega_0"])
- inner(kappa[0] * u_D * n, grad(v[0])) * ds(bound_ids["boundary_0"])
+ gamma_dg / h * inner(kappa[0] * u_D, v[0]) * ds(bound_ids["boundary_0"])
+ inner(f_1, v[1]) * dx(vol_ids["omega_1"])
+ inner(alpha[1] * u_1_n / delta_t, v[1]) * dx(vol_ids["omega_1"])
)
# Compile RHS forms
L = fem.form(ufl.extract_blocks(L), entity_maps=entity_maps)
# Apply boundary condition. Since the boundary condition is applied on
# V_1, we must convert the facet tags to submesh_1 in order to locate
# the boundary degrees of freedom.
# NOTE: We don't do this for V_0 since the Dirichlet boundary condition
# is enforced weakly by the DG scheme.
ft_sm_1 = convert_facet_tags(msh, submesh_1, sm_1_to_msh, ft)
bound_facets_sm_1 = ft_sm_1.find(bound_ids["boundary_1"])
submesh_1.topology.create_connectivity(tdim - 1, tdim)
bound_dofs = fem.locate_dofs_topological(V_1, tdim - 1, bound_facets_sm_1)
u_bc_1 = fem.Function(V_1)
u_bc_1.interpolate(u_e)
bc_1 = fem.dirichletbc(u_bc_1, bound_dofs)
bcs = [bc_1]
# Assemble the system of equations
A = assemble_matrix_block(a, bcs=bcs)
A.assemble()
b = create_vector_block(L)
# Set up solver
ksp = PETSc.KSP().create(msh.comm)
ksp.setOperators(A)
ksp.setType("preonly")
ksp.getPC().setType("lu")
# Setup files for visualisation
u_0_file = io.VTXWriter(msh.comm, "u_0.bp", [u_0_n._cpp_object], "BP4")
u_1_file = io.VTXWriter(msh.comm, "u_1.bp", [u_1_n._cpp_object], "BP4")
# Time stepping loop
t = 0.0
u_0_file.write(t)
u_1_file.write(t)
x = A.createVecRight()
for n in range(num_time_steps):
t += delta_t
with b.localForm() as b_loc:
b_loc.set(0.0)
assemble_vector_block(b, L, a, bcs=bcs)
# Compute solution
ksp.solve(b, x)
# Recover solution
offset = V_0.dofmap.index_map.size_local * V_0.dofmap.index_map_bs
u_0_n.x.array[:offset] = x.array_r[:offset]
u_1_n.x.array[: (len(x.array_r) - offset)] = x.array_r[offset:]
u_0_n.x.scatter_forward()
u_1_n.x.scatter_forward()
# Write to file
u_0_file.write(t)
u_1_file.write(t)
u_0_file.close()
u_1_file.close()
# Compute errors
e_L2_0 = norm_L2(msh.comm, u_0_n - u_e(ufl.SpatialCoordinate(submesh_0), module=ufl))
e_L2_1 = norm_L2(msh.comm, u_1_n - u_e(ufl.SpatialCoordinate(submesh_1), module=ufl))
e_L2 = np.sqrt(e_L2_0**2 + e_L2_1**2)
if msh.comm.rank == 0:
print(f"e_L2 = {e_L2}")