| 1 |
# $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() |
Parent Directory
| 
Revision Log