/[escript]/trunk/esys2/modellib/py_src/darcy.py
ViewVC logotype

Annotation of /trunk/esys2/modellib/py_src/darcy.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 127 - (hide annotations)
Fri Jul 22 05:11:29 2005 UTC (15 years, 4 months ago) by jgs
File MIME type: text/x-python
File size: 4049 byte(s)
moved modellib code to esys2/modellib/py_src

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     def doInitialization(self,t):
30     self.tn=t
31     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()

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26