Source Code for Module hedge.models.diffusion

  1  # -*- coding: utf8 -*- 
  2  """Operators modeling diffusive 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  import numpy 
 26  import numpy.linalg as la 
 27   
 28  import hedge.data 
 29  from hedge.models import TimeDependentOperator 
 30   
 31   
 32   
 33   
 34 -class StrongHeatOperator(TimeDependentOperator): 
 35 -    def __init__(self, dimensions, coeff=hedge.data.ConstantGivenFunction(1), 
 36              dirichlet_bc=hedge.data.ConstantGivenFunction(), dirichlet_tag="dirichlet", 
 37              neumann_bc=hedge.data.ConstantGivenFunction(), neumann_tag="neumann", 
 38              ldg=True): 
 39          self.dimensions = dimensions 
 40          assert isinstance(dimensions, int) 
 41   
 42          self.coeff = coeff 
 43          self.ldg = ldg 
 44   
 45          self.dirichlet_bc = dirichlet_bc 
 46          self.dirichlet_tag = dirichlet_tag 
 47          self.neumann_bc = neumann_bc 
 48          self.neumann_tag = neumann_tag 
 49   
 50      # fluxes ------------------------------------------------------------------ 
 51 -    def get_weak_flux_set(self, ldg): 
 52          class FluxSet: pass 
 53          fs = FluxSet() 
 54   
 55          from hedge.flux import FluxVectorPlaceholder, make_normal 
 56   
 57          # note here: 
 58   
 59          # local DG is unlike the other kids in that the computation of the flux 
 60          # of u depends *only* on u, whereas the computation of the flux of v 
 61          # (yielding the final right hand side) may also depend on u. That's why 
 62          # we use the layout [u,v], where v is simply omitted for the u flux 
 63          # computation. 
 64   
 65          dim = self.dimensions 
 66          vec = FluxVectorPlaceholder(1+dim) 
 67          fs.u = u = vec[0] 
 68          fs.v = v = vec[1:] 
 69          normal = fs.normal = make_normal(dim) 
 70   
 71          # central 
 72          fs.flux_u = u.avg*normal 
 73          fs.flux_v = numpy.dot(v.avg, normal) 
 74   
 75          # dbdry is "dirichlet boundary" 
 76          # nbdry is "neumann boundary" 
 77          fs.flux_u_dbdry = fs.flux_u 
 78          fs.flux_u_nbdry = fs.flux_u 
 79   
 80          fs.flux_v_dbdry = fs.flux_v 
 81          fs.flux_v_nbdry = fs.flux_v 
 82   
 83          if ldg: 
 84              ldg_beta = numpy.ones((dim,)) 
 85   
 86              fs.flux_u = fs.flux_u - (u.int-u.ext)*0.5*ldg_beta 
 87              fs.flux_v = fs.flux_v + numpy.dot((v.int-v.ext)*0.5, ldg_beta) 
 88   
 89          return fs 
 90   
 91 -    def get_strong_flux_set(self, ldg): 
 92          fs = self.get_weak_flux_set(ldg) 
 93   
 94          u = fs.u 
 95          v = fs.v 
 96          normal = fs.normal 
 97   
 98          fs.flux_u = u.int*normal - fs.flux_u 
 99          fs.flux_v = numpy.dot(v.int, normal) - fs.flux_v 
100          fs.flux_u_dbdry = u.int*normal - fs.flux_u_dbdry 
101          fs.flux_v_dbdry = numpy.dot(v.int, normal) - fs.flux_v_dbdry 
102          fs.flux_u_nbdry = u.int*normal - fs.flux_u_nbdry 
103          fs.flux_v_nbdry = numpy.dot(v.int, normal) - fs.flux_v_nbdry 
104   
105          return fs 
106   
107      # right-hand side --------------------------------------------------------- 
108 -    def grad_op_template(self): 
109          from hedge.optemplate import Field, pair_with_boundary, \ 
110                  InverseMassOperator, make_stiffness, get_flux_operator 
111   
112          stiff = make_stiffness(self.dimensions) 
113   
114          u = Field("u") 
115          sqrt_coeff_u = Field("sqrt_coeff_u") 
116          dir_bc_u = Field("dir_bc_u") 
117          neu_bc_u = Field("neu_bc_u") 
118   
119          fs = self.get_strong_flux_set(self.ldg) 
120          flux_u = get_flux_operator(fs.flux_u) 
121          flux_u_dbdry = get_flux_operator(fs.flux_u_dbdry) 
122          flux_u_nbdry = get_flux_operator(fs.flux_u_nbdry) 
123   
124          return InverseMassOperator() * ( 
125                  stiff * u 
126                  - flux_u*sqrt_coeff_u 
127                  - flux_u_dbdry*pair_with_boundary(sqrt_coeff_u, dir_bc_u, self.dirichlet_tag) 
128                  - flux_u_nbdry*pair_with_boundary(sqrt_coeff_u, neu_bc_u, self.neumann_tag) 
129                  ) 
130   
131 -    def div_op_template(self): 
132          from hedge.optemplate import make_vector_field, pair_with_boundary, \ 
133                  InverseMassOperator, get_flux_operator, make_stiffness 
134   
135          stiff = make_stiffness(self.dimensions) 
136   
137          d = self.dimensions 
138          w = make_vector_field("w", 1+d) 
139          v = w[1:] 
140   
141          dir_bc_w = make_vector_field("dir_bc_w", 1+d) 
142          neu_bc_w = make_vector_field("neu_bc_w", 1+d) 
143   
144          fs = self.get_strong_flux_set(self.ldg) 
145          flux_v = get_flux_operator(fs.flux_v) 
146          flux_v_dbdry = get_flux_operator(fs.flux_v_dbdry) 
147          flux_v_nbdry = get_flux_operator(fs.flux_v_nbdry) 
148   
149          return InverseMassOperator() * ( 
150                  numpy.dot(stiff, v) 
151                  - flux_v * w 
152                  - flux_v_dbdry * pair_with_boundary(w, dir_bc_w, self.dirichlet_tag) 
153                  - flux_v_nbdry * pair_with_boundary(w, neu_bc_w, self.neumann_tag) 
154                  ) 
155   
156      # boundary conditions ----------------------------------------------------- 
157 -    def bind(self, discr): 
158          from hedge.mesh import check_bc_coverage 
159          check_bc_coverage(discr.mesh, [self.dirichlet_tag, self.neumann_tag]) 
160   
161          return self.BoundHeatOperator(self, discr) 
162   
163 -    class BoundHeatOperator: 
164 -        def __init__(self, heat_op, discr): 
165              hop = self.heat_op = heat_op 
166              self.discr = discr 
167   
168              self.sqrt_coeff = numpy.sqrt( 
169                      hop.coeff.volume_interpolant(discr)) 
170              self.dir_sqrt_coeff = numpy.sqrt( 
171                      hop.coeff.boundary_interpolant(discr, hop.dirichlet_tag)) 
172              self.neu_sqrt_coeff = numpy.sqrt( 
173                      hop.coeff.boundary_interpolant(discr, hop.neumann_tag)) 
174   
175              self.neumann_normals = discr.boundary_normals(hop.neumann_tag) 
176   
177              self.grad_c = discr.compile(hop.grad_op_template()) 
178              self.div_c = discr.compile(hop.div_op_template()) 
179   
180 -        def dirichlet_bc_u(self, t, sqrt_coeff_u): 
181              hop = self.heat_op 
182   
183              return ( 
184                      -self.discr.boundarize_volume_field(sqrt_coeff_u, hop.dirichlet_tag) 
185                      +2*self.dir_sqrt_coeff*hop.dirichlet_bc.boundary_interpolant( 
186                          t, self.discr, hop.dirichlet_tag) 
187                      ) 
188   
189 -        def dirichlet_bc_v(self, t, sqrt_coeff_v): 
190              hop = self.heat_op 
191   
192              return self.discr.boundarize_volume_field(sqrt_coeff_v, hop.dirichlet_tag) 
193   
194 -        def neumann_bc_u(self, t, sqrt_coeff_u): 
195              hop = self.heat_op 
196   
197              return self.discr.boundarize_volume_field(sqrt_coeff_u, hop.neumann_tag) 
198   
199 -        def neumann_bc_v(self, t, sqrt_coeff_v): 
200              hop = self.heat_op 
201   
202              from hedge.tools import to_obj_array 
203              return ( 
204                      -self.discr.boundarize_volume_field(sqrt_coeff_v, hop.neumann_tag) 
205                      + 
206                      2*self.neumann_normals* 
207                          hop.neumann_bc.boundary_interpolant( 
208                              t, self.discr, hop.neumann_tag)) 
209   
210 -        def __call__(self, t, u): 
211              from hedge.tools import join_fields 
212   
213              hop = self.heat_op 
214   
215              dtag = hop.dirichlet_tag 
216              ntag = hop.neumann_tag 
217   
218              sqrt_coeff_u = self.sqrt_coeff * u 
219   
220              dir_bc_u = self.dirichlet_bc_u(t, sqrt_coeff_u) 
221              neu_bc_u = self.neumann_bc_u(t, sqrt_coeff_u) 
222   
223              v = self.grad_c( 
224                      u=u, sqrt_coeff_u=sqrt_coeff_u, 
225                      dir_bc_u=dir_bc_u, neu_bc_u=neu_bc_u) 
226   
227              from hedge.tools import ptwise_mul 
228              sqrt_coeff_v = ptwise_mul(self.sqrt_coeff, v) 
229   
230              dir_bc_v = self.dirichlet_bc_v(t, sqrt_coeff_v) 
231              neu_bc_v = self.neumann_bc_v(t, sqrt_coeff_v) 
232   
233              w = join_fields(sqrt_coeff_u, sqrt_coeff_v) 
234              dir_bc_w = join_fields(dir_bc_u, dir_bc_v) 
235              neu_bc_w = join_fields(neu_bc_u, neu_bc_v) 
236   
237              return self.div_c(w=w, dir_bc_w=dir_bc_w, neu_bc_w=neu_bc_w) 
238