/[escript]/trunk/escriptcore/py_src/flows.py
ViewVC logotype

Diff of /trunk/escriptcore/py_src/flows.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3905 by gross, Tue Jun 5 08:33:41 2012 UTC revision 3909 by gross, Wed Jun 13 07:03:06 2012 UTC
# Line 109  class DarcyFlow(object): Line 109  class DarcyFlow(object):
109    
110        self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))        self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
111        self.__g=escript.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))        self.__g=escript.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
112          self.__permeability_invXg=escript.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
113          self.__permeability_invXg_ref=util.numpy.zeros((self.domain.getDim()),util.numpy.float64)
114          self.ref_point_id=None
115          self.ref_point=util.numpy.zeros((self.domain.getDim()),util.numpy.float64)
116        self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))        self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
117        self.location_of_fixed_flux = escript.Vector(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))        self.location_of_fixed_flux = escript.Vector(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
118        self.perm_scale=1.        self.perm_scale=1.
# Line 137  class DarcyFlow(object): Line 141  class DarcyFlow(object):
141    
142        """        """
143        if location_of_fixed_pressure!=None:        if location_of_fixed_pressure!=None:
144             self.location_of_fixed_pressure=util.wherePositive(location_of_fixed_pressure)             self.location_of_fixed_pressure=util.wherePositive(util.interpolate(location_of_fixed_pressure, self.__pde_p.getFunctionSpaceForCoefficient("q")))
145               self.ref_point_id=self.location_of_fixed_pressure.maxGlobalDataPoint()
146               if not self.location_of_fixed_pressure.getTupleForGlobalDataPoint(*self.ref_point_id)[0] > 0: raise ValueError("pressure needs to be fixed at least one point.")
147               self.ref_point=self.__pde_p.getFunctionSpaceForCoefficient("q").getX().getTupleForGlobalDataPoint(*self.ref_point_id)
148               if self.verbose: print(("DarcyFlow: reference point at %s."%(self.ref_point,)))
149             self.__pde_p.setValue(q=self.location_of_fixed_pressure)             self.__pde_p.setValue(q=self.location_of_fixed_pressure)
150        if location_of_fixed_flux!=None:        if location_of_fixed_flux!=None:
151            self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)            self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
# Line 184  class DarcyFlow(object): Line 192  class DarcyFlow(object):
192          else:          else:
193               if not g.getShape()==(self.domain.getDim(),): raise ValueError("illegal shape of g")               if not g.getShape()==(self.domain.getDim(),): raise ValueError("illegal shape of g")
194          self.__g=g          self.__g=g
195            self.__permeability_invXg=util.tensor_mult(self.__permeability_inv,self.__g * (1./self.perm_scale ))
196            self.__permeability_invXg_ref=util.integrate(self.__permeability_invXg)/util.vol(self.domain)
197        if f !=None:        if f !=None:
198           f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))           f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
199           if f.isEmpty():                 if f.isEmpty():      
# Line 240  class DarcyFlow(object): Line 250  class DarcyFlow(object):
250        :rtype: ``tuple`` of `escript.Data`.        :rtype: ``tuple`` of `escript.Data`.
251    
252        """        """
253        self.__pde_p.setValue(X=self.__g * 1./self.perm_scale,        p0=util.interpolate(p0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
254          if self.ref_point_id == None:
255              p_ref=0
256          else:
257              p_ref=p0.getTupleForGlobalDataPoint(*self.ref_point_id)[0]
258          p0_hydrostatic=p_ref+util.inner(self.__permeability_invXg_ref, self.__pde_p.getFunctionSpaceForCoefficient("q").getX() - self.ref_point)
259          g_2=self.__g - util.tensor_mult(self.__permeability, self.__permeability_invXg_ref * self.perm_scale)
260          self.__pde_p.setValue(X=g_2 * 1./self.perm_scale,
261                              Y=self.__f * 1./self.perm_scale,                              Y=self.__f * 1./self.perm_scale,
262                              y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux * 1./self.perm_scale ),                              y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux * 1./self.perm_scale ),
263                              r=p0)                              r=p0 - p0_hydrostatic)
264        p=self.__pde_p.getSolution()        pp=self.__pde_p.getSolution()
265        u = self.getFlux(p, u0)        u = self._getFlux(pp, u0)
266        return u,p        return u,pp + p0_hydrostatic
267                
268     def getFlux(self,p, u0=None):     def getFlux(self,p, u0=None):
269          """          """
# Line 261  class DarcyFlow(object): Line 278  class DarcyFlow(object):
278          :return: flux          :return: flux
279          :rtype: `escript.Data`          :rtype: `escript.Data`
280          """          """
281            p=util.interpolate(p, self.__pde_p.getFunctionSpaceForCoefficient("q"))
282            if self.ref_point_id == None:
283                p_ref=0
284            else:
285                p_ref=p.getTupleForGlobalDataPoint(*self.ref_point_id)[0]
286            p_hydrostatic=p_ref+util.inner(self.__permeability_invXg_ref, self.__pde_p.getFunctionSpaceForCoefficient("q").getX() - self.ref_point)
287            return self._getFlux(p-p_hydrostatic, u0)
288    
289       def _getFlux(self,pp, u0=None):
290            """
291            returns the flux for a given pressure ``p`` where the flux is equal to ``u0``
292            on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
293            Notice that ``g`` is used, see `setValue`.
294    
295            :param p: pressure.
296            :type p: scalar value on the domain (e.g. `escript.Data`).
297            :param u0: flux on the locations of the domain marked be ``location_of_fixed_flux``.
298            :type u0: vector values on the domain (e.g. `escript.Data`) or ``None``
299            :return: flux
300            :rtype: `escript.Data`
301            """
302          if self.solver  == self.EVAL:          if self.solver  == self.EVAL:
303             u = self.__g-self.perm_scale * util.tensor_mult(self.__permeability,util.grad(p))             u = self.__g - util.tensor_mult(self.__permeability, self.perm_scale * (util.grad(pp) + self.__permeability_invXg_ref))
304          elif self.solver  == self.POST or self.solver  == self.SMOOTH:          elif self.solver  == self.POST or self.solver  == self.SMOOTH:
305              self.__pde_v.setValue(Y=util.tensor_mult(self.__permeability_inv,self.__g * 1./self.perm_scale)-util.grad(p))              self.__pde_v.setValue(Y= self.__permeability_invXg - (util.grad(pp) + self.__permeability_invXg_ref))
306              if u0 == None:              if u0 == None:
307                 self.__pde_v.setValue(r=escript.Data())                 self.__pde_v.setValue(r=escript.Data())
308              else:              else:

Legend:
Removed from v.3905  
changed lines
  Added in v.3909

  ViewVC Help
Powered by ViewVC 1.1.26