/[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 628 - (show annotations)
Thu Mar 23 02:27:57 2006 UTC (14 years, 7 months ago) by elspeth
File MIME type: text/x-python
File size: 4352 byte(s)
Copyright information added.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26