hedge.models.poisson

38 """Implements the Local Discontinuous Galerkin (LDG) Method for elliptic 39 operators. 40 41 See P. Castillo et al., 42 Local discontinuous Galerkin methods for elliptic problems", 43 Communications in Numerical Methods in Engineering 18, no. 1 (2002): 69-75. 44 """ 45

46 - def __init__(self, dimensions, diffusion_tensor=None, 47 dirichlet_bc=hedge.data.ConstantGivenFunction(), dirichlet_tag="dirichlet", 48 neumann_bc=hedge.data.ConstantGivenFunction(), neumann_tag="neumann", 49 flux="ip"):

50 """Initialize the weak Poisson operator. 51 52 @arg flux: Either C{"ip"} or C{"ldg"} to indicate which type of flux is 53 to be used. IP tends to be faster, and is therefore the default. 54 """ 55 self.dimensions = dimensions 56 assert isinstance(dimensions, int) 57 58 self.flux_type = flux 59 60 # treat diffusion tensor 61 if diffusion_tensor is None: 62 diffusion_tensor = hedge.data.ConstantGivenFunction( 63 numpy.eye(dimensions)) 64 65 self.diffusion_tensor = diffusion_tensor 66 67 self.dirichlet_bc = dirichlet_bc 68 self.dirichlet_tag = dirichlet_tag 69 self.neumann_bc = neumann_bc 70 self.neumann_tag = neumann_tag

71 72 # fluxes ------------------------------------------------------------------

73 - def get_weak_flux_set(self, flux):

74 class FluxSet: pass 75 fs = FluxSet() 76 77 if flux == "ldg": 78 ldg_terms = True 79 elif flux == "ip": 80 ldg_terms = False 81 else: 82 raise "Invalid flux type '%s'" % flux 83 84 from hedge.flux import FluxVectorPlaceholder, \ 85 make_normal, PenaltyTerm 86 from numpy import dot 87 88 dim = self.dimensions 89 vec = FluxVectorPlaceholder(1+dim) 90 fs.u = u = vec[0] 91 fs.v = v = vec[1:] 92 normal = make_normal(dim) 93 94 # central flux 95 fs.flux_u = u.avg*normal 96 fs.flux_v = dot(v.avg, normal) 97 98 if ldg_terms: 99 # ldg terms 100 ldg_beta = numpy.array([1]*dim) 101 102 fs.flux_u = fs.flux_u - (u.int-u.ext)*0.5*ldg_beta 103 fs.flux_v = fs.flux_v + dot((v.int-v.ext)*0.5, ldg_beta) 104 105 # penalty term 106 stab_term = 10 * PenaltyTerm() * (u.int - u.ext) 107 fs.flux_v -= stab_term 108 109 # boundary fluxes 110 fs.flux_u_dbdry = normal * u.ext 111 fs.flux_v_dbdry = dot(v.int, normal) - stab_term 112 113 fs.flux_u_nbdry = normal * u.int 114 fs.flux_v_nbdry = dot(normal, v.ext) 115 116 return fs

117 118 # operator application, rhs prep ------------------------------------------

119 - def grad_op_template(self):

120 from hedge.optemplate import Field, pair_with_boundary, get_flux_operator, \ 121 make_stiffness_t, InverseMassOperator 122 123 stiff_t = make_stiffness_t(self.dimensions) 124 m_inv = InverseMassOperator() 125 126 u = Field("u") 127 128 fs = self.get_weak_flux_set(self.flux_type) 129 130 flux_u = get_flux_operator(fs.flux_u) 131 flux_u_dbdry = get_flux_operator(fs.flux_u_dbdry) 132 flux_u_nbdry = get_flux_operator(fs.flux_u_nbdry) 133 134 return m_inv * ( 135 - (stiff_t * u) 136 + flux_u*u 137 + flux_u_dbdry*pair_with_boundary(u, 0, self.dirichlet_tag) 138 + flux_u_nbdry*pair_with_boundary(u, 0, self.neumann_tag) 139 )

140

141 - def div_op_template(self, apply_minv):

142 from hedge.optemplate import make_vector_field, pair_with_boundary, \ 143 make_stiffness_t, InverseMassOperator, get_flux_operator 144 145 d = self.dimensions 146 w = make_vector_field("w", 1+d) 147 v = w[1:] 148 dir_bc_w = make_vector_field("dir_bc_w", 1+d) 149 neu_bc_w = make_vector_field("neu_bc_w", 1+d) 150 151 stiff_t = make_stiffness_t(d) 152 m_inv = InverseMassOperator() 153 154 fs = self.get_weak_flux_set(self.flux_type) 155 156 flux_v = get_flux_operator(fs.flux_v) 157 flux_v_dbdry = get_flux_operator(fs.flux_v_dbdry) 158 flux_v_nbdry = get_flux_operator(fs.flux_v_nbdry) 159 160 result = ( 161 -numpy.dot(stiff_t, v) 162 + flux_v * w 163 + flux_v_dbdry * pair_with_boundary(w, dir_bc_w, self.dirichlet_tag) 164 + flux_v_nbdry * pair_with_boundary(w, neu_bc_w, self.neumann_tag) 165 ) 166 167 if apply_minv: 168 return InverseMassOperator() * result 169 else: 170 return result

171 172 @memoize_method

173 - def grad_bc_op_template(self):

174 from hedge.optemplate import Field, pair_with_boundary, \ 175 InverseMassOperator, get_flux_operator 176 177 flux_u_dbdry = get_flux_operator( 178 self.get_weak_flux_set(self.flux_type).flux_u_dbdry) 179 180 return InverseMassOperator() * ( 181 flux_u_dbdry*pair_with_boundary(0, Field("dir_bc_u"), 182 self.dirichlet_tag))

183 184 # bound operator ----------------------------------------------------------

185 - class BoundPoissonOperator(hedge.iterative.OperatorBase):

186 - def __init__(self, poisson_op, discr):

187 hedge.iterative.OperatorBase.__init__(self) 188 self.discr = discr 189 190 pop = self.poisson_op = poisson_op 191 192 self.grad_c = discr.compile(pop.grad_op_template()) 193 self.div_c = discr.compile(pop.div_op_template(False)) 194 self.minv_div_c = discr.compile(pop.div_op_template(True)) 195 self.grad_bc_c = discr.compile(pop.grad_bc_op_template()) 196 197 self.neumann_normals = discr.boundary_normals(poisson_op.neumann_tag) 198 199 if isinstance(pop.diffusion_tensor, hedge.data.ConstantGivenFunction): 200 self.diffusion = self.neu_diff = pop.diffusion_tensor.value 201 else: 202 self.diffusion = pop.diffusion_tensor.volume_interpolant(discr) 203 self.neu_diff = pop.diffusion_tensor.boundary_interpolant(discr, 204 poisson_op.neumann_tag) 205 206 # Check whether use of Poincaré mean-value method is required. 207 # This only is requested for periodic BC's over the entire domain. 208 # Partial periodic BC mixed with other BC's does not need the 209 # special treatment. 210 211 from hedge.mesh import TAG_ALL 212 self.poincare_mean_value_hack = ( 213 len(self.discr.get_boundary(TAG_ALL).nodes) 214 == len(self.discr.get_boundary(poisson_op.neumann_tag).nodes))

215 216 @property

217 - def dtype(self):

218 return self.discr.default_scalar_type

219 220 @property

221 - def shape(self):

222 nodes = len(self.discr) 223 return nodes, nodes

224 225 # actual functionality

226 - def grad(self, u):

227 return self.grad_c(u=u)

228

229 - def div(self, v, u=None, apply_minv=True):

230 """Compute the divergence of v using an LDG operator. 231 232 The divergence computation is unaffected by the scaling 233 effected by the diffusion tensor. 234 235 @param apply_minv: Bool specifying whether to compute a complete 236 divergence operator. If False, the final application of the inverse 237 mass operator is skipped. This is used in L{op}() in order to reduce 238 the scheme M{M^{-1} S u = f} to M{S u = M f}, so that the mass operator 239 only needs to be applied once, when preparing the right hand side 240 in @L{prepare_rhs}. 241 """ 242 from hedge.tools import join_fields 243 244 dim = self.discr.dimensions 245 246 if u is None: 247 u = self.discr.volume_zeros() 248 w = join_fields(u, v) 249 250 dir_bc_w = join_fields(0, [0]*dim) 251 neu_bc_w = join_fields(0, [0]*dim) 252 253 if apply_minv: 254 div_tpl = self.minv_div_c 255 else: 256 div_tpl = self.div_c 257 258 return div_tpl(w=w, dir_bc_w=dir_bc_w, neu_bc_w=neu_bc_w)

259

260 - def op(self, u, apply_minv=False):

261 from hedge.tools import ptwise_dot 262 263 # Check if poincare mean value method has to be applied. 264 if self.poincare_mean_value_hack: 265 # ∫(Ω) u dΩ 266 state_int = self.discr.integral(u) 267 # calculate mean value: (1/|Ω|) * ∫(Ω) u dΩ 268 mean_state = state_int / self.discr.mesh_volume() 269 m_mean_state = mean_state * self.discr._mass_ones() 270 #m_mean_state = mean_state * m 271 else: 272 m_mean_state = 0 273 274 return self.div( 275 ptwise_dot(2, 1, self.diffusion, self.grad(u)), 276 u, apply_minv=apply_minv) \ 277 - m_mean_state

278 279 __call__ = op 280

281 - def prepare_rhs(self, rhs):

282 """Prepare the right-hand side for the linear system op(u)=rhs(f). 283 284 In matrix form, LDG looks like this: 285 286 Mv = Cu + g 287 Mf = Av + Bu + h 288 289 where v is the auxiliary vector, u is the argument of the operator, f 290 is the result of the grad operator, g and h are inhom boundary data, and 291 A,B,C are some operator+lifting matrices. 292 293 M f = A Minv(Cu + g) + Bu + h 294 295 so the linear system looks like 296 297 M f = A Minv Cu + A Minv g + Bu + h 298 M f - A Minv g - h = (A Minv C + B)u (*) 299 --------rhs------- 300 301 So the right hand side we're putting together here is really 302 303 M f - A Minv g - h 304 305 Finally, note that the operator application above implements 306 the equation (*) left-multiplied by Minv, so that the 307 right-hand-side becomes 308 309 f - Minv( A Minv g + h) 310 """ 311 dim = self.discr.dimensions 312 313 pop = self.poisson_op 314 315 dtag = pop.dirichlet_tag 316 ntag = pop.neumann_tag 317 318 dir_bc_u = pop.dirichlet_bc.boundary_interpolant(self.discr, dtag) 319 vpart = self.grad_bc_c(dir_bc_u=dir_bc_u) 320 321 from hedge.tools import ptwise_dot 322 diff_v = ptwise_dot(2, 1, self.diffusion, vpart) 323 324 def neu_bc_v(): 325 return ptwise_dot(2, 1, self.neu_diff, 326 self.neumann_normals* 327 pop.neumann_bc.boundary_interpolant(self.discr, ntag))

328 329 from hedge.tools import join_fields 330 w = join_fields(0, diff_v) 331 dir_bc_w = join_fields(dir_bc_u, [0]*dim) 332 neu_bc_w = join_fields(0, neu_bc_v()) 333 334 from hedge.optemplate import MassOperator 335 336 return (MassOperator().apply(self.discr, 337 rhs.volume_interpolant(self.discr)) 338 - self.div_c(w=w, dir_bc_w=dir_bc_w, neu_bc_w=neu_bc_w))

Source Code for Module hedge.models.poisson