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

revision 3500 by gross, Mon Apr 11 03:09:06 2011 UTC revision 3501 by gross, Wed Apr 13 04:07:53 2011 UTC
# Line 38  from linearPDEs import LinearPDE, Linear Line 38  from linearPDEs import LinearPDE, Linear
38  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
39
40
41  class DarcyFlowVeryNew(object):
42    class DarcyFlowNew(object):
43     """     """
44     solves the problem     solves the problem
45
# Line 138  class DarcyFlowVeryNew(object): Line 139  class DarcyFlowVeryNew(object):
139               is along the *x_i* axis.               is along the *x_i* axis.
140
141        """        """
142        if location_of_fixed_pressure!=None: self.location_of_fixed_pressure=location_of_fixed_pressure        if location_of_fixed_pressure!=None: self.location_of_fixed_pressure=util.wherePositive(location_of_fixed_pressure)
143        if location_of_fixed_flux!=None: self.location_of_fixed_flux=location_of_fixed_flux        if location_of_fixed_flux!=None: self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
144
145        if self.useVPIteration:        if self.useVPIteration:
146           if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=self.location_of_fixed_pressure)           if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=self.location_of_fixed_pressure)
# Line 219  class DarcyFlowVeryNew(object): Line 220  class DarcyFlowVeryNew(object):
220        if self.useVPIteration:        if self.useVPIteration:
221          # get u:          # get u:
222              rtol=self.getTolerance()              rtol=self.getTolerance()
223          self.__pde_k.setValue(Y=0.5*util.tensor_mult(self.__permeability_inv,g),X=escript.Data())              p_init=0*p0_b
225                                      X=p_init*util.kronecker(self.__pde_k.getDomain()))
226          du=self.__pde_k.getSolution()          du=self.__pde_k.getSolution()
227          self.__pde_p.setValue(Y=f-util.div(du), X=0.5*(g-du))          self.__pde_p.setValue(Y=f-util.div(du), X=0.5*(g-du))
228          p=GMRES(self.__pde_p.getSolution(),          p=GMRES(self.__pde_p.getSolution(),
229                  self.__Aprod_GMRES,                  self.__Aprod_GMRES,
230              p0_b*0,              p_init,
231              self.__inner_GMRES,              self.__inner_GMRES,
232              atol=0,              atol=0,
233              rtol=rtol,              rtol=rtol,
# Line 278  class DarcyFlowVeryNew(object): Line 281  class DarcyFlowVeryNew(object):
281        :rtype: ``float``        :rtype: ``float``
282        """        """
283        return self.__rtol        return self.__rtol
284
285  class DarcyFlow(object):  class DarcyFlow(object):
286     """     """
287     solves the problem     solves the problem
# Line 290  class DarcyFlow(object): Line 294  class DarcyFlow(object):
294     :note: The problem is solved in a least squares formulation.     :note: The problem is solved in a least squares formulation.
295     """     """
296
297     def __init__(self, domain, useReduced=False, adaptSubTolerance=True, solveForFlux=False, useVPIteration=True, weighting_scale=0.1):     def __init__(self, domain, useReduced=False, adaptSubTolerance=True, solveForFlux=False, useVPIteration=True, weighting_scale=1.):
298        """        """
299        initializes the Darcy flux problem        initializes the Darcy flux problem
300        :param domain: domain of the problem        :param domain: domain of the problem
# Line 405  class DarcyFlow(object): Line 409  class DarcyFlow(object):
409
410        """        """
411        if self.useVPIteration:        if self.useVPIteration:
412           if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)           if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=util.wherePositive(location_of_fixed_pressure))
413           if location_of_fixed_flux!=None: self.__pde_k.setValue(q=location_of_fixed_flux)           if location_of_fixed_flux!=None: self.__pde_k.setValue(q=util.wherePositive(location_of_fixed_flux))
414        else:        else:
415           q=self.__pde_k.getCoefficient("q")           q=self.__pde_k.getCoefficient("q")
416           if q.isEmpty(): q=self.__pde_k.createCoefficient("q")           if q.isEmpty(): q=self.__pde_k.createCoefficient("q")
417           if location_of_fixed_pressure!=None: q[self.domain.getDim()]=location_of_fixed_pressure           if location_of_fixed_pressure!=None: q[self.domain.getDim()]=util.wherePositive(location_of_fixed_pressure)
418           if location_of_fixed_flux!=None: q[:self.domain.getDim()]=location_of_fixed_flux           if location_of_fixed_flux!=None: q[:self.domain.getDim()]=util.wherePositive(location_of_fixed_flux)
419           self.__pde_k.setValue(q=q)           self.__pde_k.setValue(q=q)
420
421        # flux is rescaled by the factor mean value(perm_inv)*length where length**self.domain.getDim()=vol(self.domain)        # flux is rescaled by the factor mean value(perm_inv)*length where length**self.domain.getDim()=vol(self.domain)
422          V=util.vol(self.domain)
423        if permeability!=None:        if permeability!=None:
424       perm=util.interpolate(permeability,self.__pde_k.getFunctionSpaceForCoefficient("A"))       perm=util.interpolate(permeability,self.__pde_k.getFunctionSpaceForCoefficient("A"))
425           V=util.vol(self.domain)
426       if perm.getRank()==0:       if perm.getRank()==0:
427          perm_inv=(1./perm)          perm_inv=(1./perm)
428              if self.useVPIteration:              if self.useVPIteration:
# Line 441  class DarcyFlow(object): Line 446  class DarcyFlow(object):
446       self.__permeability=perm       self.__permeability=perm
447       self.__permeability_inv=perm_inv       self.__permeability_inv=perm_inv
448
449       self.__l2 =(util.longestEdge(self.domain)**2*util.length(self.__permeability_inv))*self.weighting_scale       self.__l2 = V**(1./self.domain.getDim())*util.length(self.__permeability_inv)*self.domain.getSize()*self.weighting_scale
450           if self.useVPIteration:           if self.useVPIteration:
451          if  self.solveForFlux:          if  self.solveForFlux:
452             self.__pde_k.setValue(D=self.__permeability_inv)             self.__pde_k.setValue(D=self.__permeability_inv)

Legend:
 Removed from v.3500 changed lines Added in v.3501