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

Contents of /trunk/modellib/py_src/darcy.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 155 - (show annotations)
Wed Nov 9 02:02:19 2005 UTC (13 years, 10 months ago) by jgs
File MIME type: text/x-python
File size: 4063 byte(s)
move all directories from trunk/esys2 into trunk and remove esys2

1 # $Id$
2
3
4 from esys.escript import *
5 from esys.escript.modelframe import Model,IterationDivergenceError
6 from esys.escript.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):
30 self.tn=0
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