hedge.models.advection

1 # -*- coding: utf8 -*- 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

35 -class AdvectionOperatorBase(TimeDependentOperator):

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

55 - def weak_flux(self):

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, # outflow 70 u.ext, # inflow 71 )) 72 else: 73 raise ValueError, "invalid flux type"

74

75 - def max_eigenvalue(self):

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

116 -class StrongAdvectionOperator(AdvectionOperatorBase):

117 - def flux(self):

118 from hedge.flux import make_normal, FluxScalarPlaceholder 119 120 u = FluxScalarPlaceholder(0) 121 normal = make_normal(self.dimensions) 122 123 return u.int * numpy.dot(normal, self.v) - self.weak_flux()

124

125 - def op_template(self):

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

148 -class WeakAdvectionOperator(AdvectionOperatorBase):

149 - def flux(self):

150 return self.weak_flux()

151

152 - def op_template(self):

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 # boundary conditions ------------------------------------------------- 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

183 -class VariableCoefficientAdvectionOperator(TimeDependentOperator):

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

211 - def flux(self, ):

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, # outflow 241 numpy.dot(normal, v.ext) * u.ext, # inflow 242 )) 243 return f 244 else: 245 raise ValueError, "invalid flux type"

246

247 - def op_template(self):

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 # boundary conditions ------------------------------------------------- 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):

289 compiled_op_template = discr.compile(self.op_template()) 290 291 from hedge.mesh import check_bc_coverage, TAG_ALL 292 check_bc_coverage(discr.mesh, [TAG_ALL]) 293 294 def rhs(t, u): 295 v = self.advec_v.volume_interpolant(t, discr) 296 297 if self.bc_u_f is not "None": 298 bc_u = self.bc_u_f.boundary_interpolant(t, discr, tag=TAG_ALL) 299 return compiled_op_template(u=u, v=v, bc_u=bc_u) 300 else: 301 return compiled_op_template(u=u, v=v)

302 303 return rhs

304

305 - def max_eigenvalue(self, t, discr):

306 # Gives the max eigenvalue of a vector of eigenvalues. 307 # As the velocities of each node is stored in the velocity-vector-field 308 # a pointwise dot product of this vector has to be taken to get the 309 # magnitude of the velocity at each node. From this vector the maximum 310 # values limits the timestep. 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

Source Code for Module hedge.models.advection