modellib/py_src/darcy.py

# $Id$


from esys.escript import *
from esys.modelframe import Model,IterationDivergenceError
from esys.linearPDEs import AdvectivePDE,LinearPDE

class DarcyFlow(Model):
       """ """

       def __init__(self,debug=False):
           Model.__init__(self,debug=debug)
           self.declareParameter(domain=None, \
                                 tend=0., \
                                 dt=0.1, \
                                 temperature=1., \
                                 density=1., \
                                 c_p=1., \
                                 thermal_permabilty=1., \
                                 reference_temperature=0., \
                                 radiation_coefficient=0., \
                                 thermal_source=0., \
                                 location_fixed_temperature=Data(), \
                                 iterate=True, \
                                 tol=1.e-8, \
                                 implicit=True)
           self.iter=0

       def doInitialization(self):
           self.tn=0
           self.pde=LinearPDE(self.domain)
           self.pde.setSymmetryOn()

       def doIterationInitialization(self,dt):
            self.iter=0
            self.T_last=self.temperature
            self.diff=1.e400

       def doIterationStep(self,dt):
          T=self.temperature
          diff=self.diff
          dim=self.pde.getDim()
          self.iter+=1
          rhocp=self.density*self.c_p
          self.pde.setValue(A=self.thermal_permabilty*kronecker(dim), \
                              D=rhocp/dt, \
                              Y=self.thermal_source+rhocp/dt*self.T_last, \
                              d=self.radiation_coefficient, \
                              y=self.radiation_coefficient*self.reference_temperature, \
                              q=self.location_fixed_temperature, \
                              r=self.T_last)
          if isinstance(self,TemperatureAdvection): self.pde.setValue(C=self.velocity[:dim]*rhocp)
          self.pde.setTolerance(self.tol*1.e-2)
          self.temperature=self.pde.getSolution()
          self.diff=Lsup(T-self.temperature)
          if diff<=self.diff:
              raise IterationDivergenceError,"no improvement in the temperature iteration"

       def terminate(self):
          if not self.implicit:
              return True
          elif self.iter<1:
              return False
          else:
             return self.diff<self.tol*Lsup(self.temperature)

       def doIterationFinalization(self,dt):
          self.tn+=dt

       def getSafeTimeStepSize(self,dt):
           return self.dt

       def finalize(self):
            return self.tn>=self.tend

if __name__=="__main__":
   from esys.modelframe import Link,Simulation,ExplicitSimulation
   from esys.visualization import WriteVTK
   from esys.materials import SimpleMaterialTable
   from esys.geometry import RectangularDomain,ScalarConstrainer
   from esys.input import InterpolatedTimeProfile,GausseanProfile

   dom=RectangularDomain()
   constraints=ScalarConstrainer()
   constraints.domain=Link(dom)
   constraints.top=1
   constraints.bottom=1
   
   mt=MaterialTable()
   
   pf=InterpolatedTimeProfile()
   pf.t=[0.,0.25,0.5,0.75]
   pf.values=[0.,1.,1.,0.]
   
   q=GausseanProfile()
   q.domain=Link(dom)
   q.width=0.05
   q.x_c=numarray.array([0.5,0.5,0.5])
   q.r=0.01
   q.A=Link(pf,"out")
   
   tt=DarcyFlow()
   tt.domain=Link(dom)
   tt.tend=1.
   tt.dt=0.1
   tt.temperature=0.
   tt.density=Link(mt)
   tt.c_p=Link(mt)
   tt.thermal_permabilty=Link(mt)
   tt.reference_temperature=0.
   tt.radiation_coefficient=Link(mt)
   tt.thermal_source=Link(q,"out")
   tt.location_fixed_temperature=Link(constraints,"location_of_constraint")
   tt.implicit=True
   
   vis=WriteVTK()
   vis.scalar=Link(tt,"temperature")
   
   s=ExplicitSimulation([dom,constraints,pf,q,Simulation([mt,tt],debug=True),vis],debug=True)
   # s=Simulation([dom,constraints,pf,q,Simulation([mt,tt]),vis])
   s.writeXML()
   s.run()
1	jgs	127	# $Id$
2
3
4			from esys.escript import *
5			from esys.modelframe import Model,IterationDivergenceError
6			from esys.linearPDEs import AdvectivePDE,LinearPDE
7
8			class DarcyFlow(Model):
9			""" """
10
11			def __init__(self,debug=False):
12			Model.__init__(self,debug=debug)
13			self.declareParameter(domain=None, \
14			tend=0., \
15			dt=0.1, \
16			temperature=1., \
17			density=1., \
18			c_p=1., \
19			thermal_permabilty=1., \
20			reference_temperature=0., \
21			radiation_coefficient=0., \
22			thermal_source=0., \
23			location_fixed_temperature=Data(), \
24			iterate=True, \
25			tol=1.e-8, \
26			implicit=True)
27			self.iter=0
28
29	jgs	147	def doInitialization(self):
30			self.tn=0
31	jgs	127	self.pde=LinearPDE(self.domain)
32			self.pde.setSymmetryOn()
33
34			def doIterationInitialization(self,dt):
35			self.iter=0
36			self.T_last=self.temperature
37			self.diff=1.e400
38
39			def doIterationStep(self,dt):
40			T=self.temperature
41			diff=self.diff
42			dim=self.pde.getDim()
43			self.iter+=1
44			rhocp=self.density*self.c_p
45			self.pde.setValue(A=self.thermal_permabilty*kronecker(dim), \
46			D=rhocp/dt, \
47			Y=self.thermal_source+rhocp/dt*self.T_last, \
48			d=self.radiation_coefficient, \
49			y=self.radiation_coefficient*self.reference_temperature, \
50			q=self.location_fixed_temperature, \
51			r=self.T_last)
52			if isinstance(self,TemperatureAdvection): self.pde.setValue(C=self.velocity[:dim]*rhocp)
53			self.pde.setTolerance(self.tol*1.e-2)
54			self.temperature=self.pde.getSolution()
55			self.diff=Lsup(T-self.temperature)
56			if diff<=self.diff:
57			raise IterationDivergenceError,"no improvement in the temperature iteration"
58
59			def terminate(self):
60			if not self.implicit:
61			return True
62			elif self.iter<1:
63			return False
64			else:
65			return self.diff<self.tol*Lsup(self.temperature)
66
67			def doIterationFinalization(self,dt):
68			self.tn+=dt
69
70			def getSafeTimeStepSize(self,dt):
71			return self.dt
72
73			def finalize(self):
74			return self.tn>=self.tend
75
76			if __name__=="__main__":
77			from esys.modelframe import Link,Simulation,ExplicitSimulation
78			from esys.visualization import WriteVTK
79			from esys.materials import SimpleMaterialTable
80			from esys.geometry import RectangularDomain,ScalarConstrainer
81			from esys.input import InterpolatedTimeProfile,GausseanProfile
82
83			dom=RectangularDomain()
84			constraints=ScalarConstrainer()
85			constraints.domain=Link(dom)
86			constraints.top=1
87			constraints.bottom=1
88
89			mt=MaterialTable()
90
91			pf=InterpolatedTimeProfile()
92			pf.t=[0.,0.25,0.5,0.75]
93			pf.values=[0.,1.,1.,0.]
94
95			q=GausseanProfile()
96			q.domain=Link(dom)
97			q.width=0.05
98			q.x_c=numarray.array([0.5,0.5,0.5])
99			q.r=0.01
100			q.A=Link(pf,"out")
101
102			tt=DarcyFlow()
103			tt.domain=Link(dom)
104			tt.tend=1.
105			tt.dt=0.1
106			tt.temperature=0.
107			tt.density=Link(mt)
108			tt.c_p=Link(mt)
109			tt.thermal_permabilty=Link(mt)
110			tt.reference_temperature=0.
111			tt.radiation_coefficient=Link(mt)
112			tt.thermal_source=Link(q,"out")
113			tt.location_fixed_temperature=Link(constraints,"location_of_constraint")
114			tt.implicit=True
115
116			vis=WriteVTK()
117			vis.scalar=Link(tt,"temperature")
118
119			s=ExplicitSimulation([dom,constraints,pf,q,Simulation([mt,tt],debug=True),vis],debug=True)
120			# s=Simulation([dom,constraints,pf,q,Simulation([mt,tt]),vis])
121			s.writeXML()
122			s.run()
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision