1
2 """Operators modeling advective phenomena."""
3
4 from __future__ import division
5
6 __copyright__ = "Copyright (C) 2009 Andreas Kloeckner"
7
8 __license__ = """
9 This program is free software: you can redistribute it and/or modify
10 it under the terms of the GNU General Public License as published by
11 the Free Software Foundation, either version 3 of the License, or
12 (at your option) any later version.
13
14 This program is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 GNU General Public License for more details.
18
19 You should have received a copy of the GNU General Public License
20 along with this program. If not, see U{http://www.gnu.org/licenses/}.
21 """
22
23
24
25
26 import numpy
27 import numpy.linalg as la
28
29 from hedge.models import TimeDependentOperator
30 import hedge.data
31
32
33
34
36 flux_types = [
37 "central",
38 "upwind",
39 "lf"
40 ]
41
42 - def __init__(self, v,
43 inflow_tag="inflow",
44 inflow_u=hedge.data.make_tdep_constant(0),
45 outflow_tag="outflow",
46 flux_type="central"
47 ):
48 self.dimensions = len(v)
49 self.v = v
50 self.inflow_tag = inflow_tag
51 self.inflow_u = inflow_u
52 self.outflow_tag = outflow_tag
53 self.flux_type = flux_type
54
56 from hedge.flux import make_normal, FluxScalarPlaceholder, IfPositive
57
58 u = FluxScalarPlaceholder(0)
59 normal = make_normal(self.dimensions)
60
61 if self.flux_type == "central":
62 return u.avg*numpy.dot(normal, self.v)
63 elif self.flux_type == "lf":
64 return u.avg*numpy.dot(normal, self.v) \
65 + 0.5*la.norm(self.v)*(u.int - u.ext)
66 elif self.flux_type == "upwind":
67 return (numpy.dot(normal, self.v)*
68 IfPositive(numpy.dot(normal, self.v),
69 u.int,
70 u.ext,
71 ))
72 else:
73 raise ValueError, "invalid flux type"
74
76 return la.norm(self.v)
77
78 - def bind(self, discr):
79 compiled_op_template = discr.compile(self.op_template())
80
81 from hedge.mesh import check_bc_coverage
82 check_bc_coverage(discr.mesh, [self.inflow_tag, self.outflow_tag])
83
84 def rhs(t, u):
85 bc_in = self.inflow_u.boundary_interpolant(t, discr, self.inflow_tag)
86 return compiled_op_template(u=u, bc_in=bc_in)
87
88 return rhs
89
90 - def bind_interdomain(self,
91 my_discr, my_part_data,
92 nb_discr, nb_part_data):
93 from hedge.partition import compile_interdomain_flux
94 compiled_op_template, from_nb_indices = compile_interdomain_flux(
95 self.op_template(), "u", "nb_bdry_u",
96 my_discr, my_part_data, nb_discr, nb_part_data,
97 use_stupid_substitution=True)
98
99 from hedge.tools import with_object_array_or_scalar, is_zero
100
101 def nb_bdry_permute(fld):
102 if is_zero(fld):
103 return 0
104 else:
105 return fld[from_nb_indices]
106
107 def rhs(t, u, u_neighbor):
108 return compiled_op_template(u=u,
109 nb_bdry_u=with_object_array_or_scalar(nb_bdry_permute, u_neighbor))
110
111 return rhs
112
113
114
115
124
126 from hedge.optemplate import Field, pair_with_boundary, \
127 get_flux_operator, make_nabla, InverseMassOperator
128
129 u = Field("u")
130 bc_in = Field("bc_in")
131
132 nabla = make_nabla(self.dimensions)
133 m_inv = InverseMassOperator()
134
135 flux_op = get_flux_operator(self.flux())
136
137 return (
138 -numpy.dot(self.v, nabla*u)
139 + m_inv*(
140 flux_op * u
141 + flux_op * pair_with_boundary(u, bc_in, self.inflow_tag)
142 )
143 )
144
145
146
147
151
153 from hedge.optemplate import \
154 Field, \
155 pair_with_boundary, \
156 get_flux_operator, \
157 make_minv_stiffness_t, \
158 InverseMassOperator, \
159 BoundarizeOperator
160
161 u = Field("u")
162
163
164 from hedge.optemplate import BoundarizeOperator
165 bc_in = Field("bc_in")
166 bc_out = BoundarizeOperator(self.outflow_tag)*u
167
168 minv_st = make_minv_stiffness_t(self.dimensions)
169 m_inv = InverseMassOperator()
170
171 flux_op = get_flux_operator(self.flux())
172
173 return numpy.dot(self.v, minv_st*u) - m_inv*(
174 flux_op*u
175 + flux_op * pair_with_boundary(u, bc_in, self.inflow_tag)
176 + flux_op * pair_with_boundary(u, bc_out, self.outflow_tag)
177 )
178
179
180
181
182
184 """A class for space- and time-dependent DG-advection operators.
185
186 `advec_v` is a callable expecting two arguments `(x, t)` representing space and time,
187 and returning an n-dimensional vector representing the velocity at x.
188 `bc_u_f` is a callable expecting `(x, t)` representing space and time,
189 and returning an 1-dimensional vector representing the state on the boundary.
190 Both `advec_v` and `bc_u_f` conform to the
191 `hedge.data.ITimeDependentGivenFunction` interface.
192 """
193
194 flux_types = [
195 "central",
196 "upwind",
197 "lf"
198 ]
199
200 - def __init__(self,
201 dimensions,
202 advec_v,
203 bc_u_f="None",
204 flux_type="central"
205 ):
206 self.dimensions = dimensions
207 self.advec_v = advec_v
208 self.bc_u_f = bc_u_f
209 self.flux_type = flux_type
210
212 from hedge.flux import \
213 make_normal, \
214 FluxScalarPlaceholder, \
215 FluxVectorPlaceholder, \
216 IfPositive, flux_max, norm
217
218 d = self.dimensions
219
220 w = FluxVectorPlaceholder((1+d)+1)
221 u = w[0]
222 v = w[1:d+1]
223 c = w[1+d]
224
225 normal = make_normal(self.dimensions)
226
227 if self.flux_type == "central":
228 return (u.int*numpy.dot(v.int, normal )
229 + u.ext*numpy.dot(v.ext, normal)) * 0.5
230 elif self.flux_type == "lf":
231 n_vint = numpy.dot(normal, v.int)
232 n_vext = numpy.dot(normal, v.ext)
233 return 0.5 * (n_vint * u.int + n_vext * u.ext) \
234 - 0.5 * (u.ext - u.int) \
235 * flux_max(c.int, c.ext)
236
237 elif self.flux_type == "upwind":
238 return (
239 IfPositive(numpy.dot(normal, v.avg),
240 numpy.dot(normal, v.int) * u.int,
241 numpy.dot(normal, v.ext) * u.ext,
242 ))
243 return f
244 else:
245 raise ValueError, "invalid flux type"
246
248 from hedge.optemplate import \
249 Field, \
250 pair_with_boundary, \
251 get_flux_operator, \
252 make_minv_stiffness_t, \
253 InverseMassOperator,\
254 make_vector_field, \
255 ElementwiseMaxOperator, \
256 BoundarizeOperator
257
258
259 from hedge.tools import join_fields, \
260 ptwise_dot
261
262 u = Field("u")
263 v = make_vector_field("v", self.dimensions)
264 c = ElementwiseMaxOperator()*ptwise_dot(1, 1, v, v)
265 w = join_fields(u, v, c)
266
267
268 from hedge.mesh import TAG_ALL
269 bc_c = BoundarizeOperator(TAG_ALL) * c
270 bc_u = Field("bc_u")
271 bc_v = BoundarizeOperator(TAG_ALL) * v
272 if self.bc_u_f is "None":
273 bc_w = join_fields(0, bc_v, bc_c)
274 else:
275 bc_w = join_fields(bc_u, bc_v, bc_c)
276
277 minv_st = make_minv_stiffness_t(self.dimensions)
278 m_inv = InverseMassOperator()
279
280 flux_op = get_flux_operator(self.flux())
281
282 result = numpy.dot(minv_st, v*u) - m_inv*(
283 flux_op * w
284 + flux_op * pair_with_boundary(w, bc_w, TAG_ALL)
285 )
286 return result
287
288 - def bind(self, discr):
302
303 return rhs
304
306
307
308
309
310
311
312 from hedge.tools import ptwise_dot
313 v = self.advec_v.volume_interpolant(t, discr)
314 return (ptwise_dot(1, 1, v, v)**0.5).max()
315