/[escript]/trunk/finley/test/python/localization.py
ViewVC logotype

Contents of /trunk/finley/test/python/localization.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3283 - (show annotations)
Mon Oct 18 22:39:28 2010 UTC (10 years, 4 months ago) by gross
File MIME type: text/x-python
File size: 13985 byte(s)
AMG reengineered, BUG is Smoother fixed.


1 ########################################################
2 #
3 # Copyright (c) 2003-2010 by University of Queensland
4 # Earth Systems Science Computational Center (ESSCC)
5 # http://www.uq.edu.au/esscc
6 #
7 # Primary Business: Queensland, Australia
8 # Licensed under the Open Software License version 3.0
9 # http://www.opensource.org/licenses/osl-3.0.php
10 #
11 ########################################################
12 """
13 this is a localization a simulation over domain [0,L] X [0,L] x [0,H]
14 with a plastic layer above a viscous layer of thickness H_VISC.
15 The yield condition is perturbed along a line at the boundary between
16 viscous and plastic layer to trigger localization.
17 """
18 __copyright__="""Copyright (c) 2003-2010 by University of Queensland
19 Earth Systems Science Computational Center (ESSCC)
20 http://www.uq.edu.au/esscc
21 Primary Business: Queensland, Australia"""
22 __license__="""Licensed under the Open Software License version 3.0
23 http://www.opensource.org/licenses/osl-3.0.php"""
24 __url__="https://launchpad.net/escript-finley"
25
26 from esys.escript import *
27 from esys.escript.unitsSI import DEG
28 from esys.escript.models import StokesProblemCartesian
29 from esys.finley import Rectangle, Brick, LoadMesh
30 from math import pi, ceil
31 import sys
32 import time
33
34 # ======================= Default Values ==================================================
35 DIM=2 # spatial dimension
36 H=1. # height
37 L=4*H # length
38 NE=10 # number of elements in H-direction.
39 H_VISC=0.2*H # height of viscous zone (must aline with element boundary)
40
41 ETA_TOP=10.
42 ETA_BOTTOM=ETA_TOP/10.
43 FRICTION_ANGLE=70*DEG
44 C=10000. # =None interesting value is calculated.
45 V0=0.5 # =None interesting value is calculated.
46 COMPRESSION=False
47 RHO=1.
48 G=1.*0
49 W_WEAK=0.03*L # width of weak zone (>H/NE)
50 ALPHA_WEAK=65*DEG # angle of week zone against x-axis
51 OFFSET_X_WEAK=L*0.3 # offset of weak zone
52 OFFSET_Z_WEAK=H_VISC # offset of weak zone
53 WEAK_FRACTION=1.0 # reduction factor for weak zone
54 STOP_FRAC=0.2 # iteratiun stops when height ~ (1-STOP_FRAC) * H or height ~ (1.+STOP_FRAC)*H
55 RESTART = False
56 VIS_DIR="./data"
57 TOL=1.e-4
58 MAX_ITER=100
59 VERBOSE=True
60
61 #
62 # derived values:
63 #
64 if C==None:
65 C=RHO*G*H*(1-sin(FRICTION_ANGLE))/cos(FRICTION_ANGLE)
66 if V0==None:
67 V0=RHO*G*H**2/ETA_TOP
68 if COMPRESSION:
69 DIRECTION=-1.
70 T_END=(STOP_FRAC+2)/(STOP_FRAC+1)*L/V0
71 else:
72 DIRECTION=1.
73 T_END=STOP_FRAC/(1-STOP_FRAC)*L/V0
74 FRICTION_ANGLE_WEAK=FRICTION_ANGLE*WEAK_FRACTION
75 C_WEAK=C*WEAK_FRACTION
76 NE_L=int(ceil(L/H*NE))
77 H_VISC=int((H_VISC*NE)/H+0.5)*(H/NE)
78 OFFSET_X_WEAK=int((OFFSET_X_WEAK*NE_L)/L+0.5)*(L/NE_L)
79 OFFSET_Z_WEAK=int((OFFSET_Z_WEAK*NE)/H+0.5)*(H/NE)
80 #
81 # print input
82 #
83 print "spatial dimenstion DIM = ",DIM
84 print "height H =",H
85 print "length L = ",L
86 print "vertical number of element NE = ", NE
87 print "horizontal number of element NE_L = ", NE_L
88 print "total number of element = ", NE_L**(DIM-1)*NE
89 print "height of viscosity layer H_VISC = ",H_VISC
90 print "viscosity in top layer ETA_TOP = ",ETA_TOP
91 print "viscosity in bottom/viscous layer ETA_BOTTOM = ",ETA_BOTTOM
92 print "friction angle FRICTION_ANGLE =",FRICTION_ANGLE/DEG
93 print "C =",C
94 print "width of weak zone W_WEAK = ",W_WEAK
95 print "elements in weak zone = ",W_WEAK/H*NE
96 print "strike of weak zone ALPHA_WEAK = ",ALPHA_WEAK/DEG
97 print "x-offset of weak zone OFFSET_X_WEAK = ",OFFSET_X_WEAK
98 print "z-offset of weak zone OFFSET_Z_WEAK = ",OFFSET_Z_WEAK
99 print "friction angle in weak zone FRICTION_ANGLE_WEAK =",FRICTION_ANGLE_WEAK/DEG
100 print "C in weak zone C_WEAK =",C_WEAK
101 print "density RHO = ",RHO
102 print "Gravity G = ",G
103 if COMPRESSION:
104 print "Direction: compression"
105 else:
106 print "Direction: expansion"
107 print "surface velocity = ",V0
108 print "end time T_END = ",T_END
109 mkDir(VIS_DIR)
110 #
111 # set up geomtry:
112 #
113 #
114 if RESTART:
115 dom=LoadMesh("mesh.nc")
116 if not dom.getDim() == DIM:
117 raise ValueError,"Expected DIM and dimension of mesh in restart file don't match."
118
119 FF=open("stamp.csv","r").read().split(",")
120 t=float(FF[0]) # time stamp
121 t_vis=float(FF[1]) #
122 n_vis=int(FF[2])
123 n=int(FF[3]) # time step counter
124 counter_vis=int(FF[4])
125 else:
126 if DIM==2:
127 dom=Rectangle(NE_L,NE,l0=L,l1=H,order=-1,optimize=True)
128 else:
129 dom=Brick(NE_L,NE_L,NE,l0=L,l1=L,l2=H,order=-1,optimize=True)
130 t=0
131 n=0
132 t_vis=0
133 n_vis=0
134 counter_vis=0
135 #
136 # material properties:
137 #
138 x=Function(dom).getX()
139 mask_visc=whereNegative(x[DIM-1]-H_VISC)
140 if DIM == 3:
141 offset=[OFFSET_X_WEAK, 0, OFFSET_Z_WEAK]
142 strike=numpy.array([cos(ALPHA_WEAK),sin(ALPHA_WEAK),0.])
143 d2=length(x-offset)**2-inner(x-offset,strike)**2
144 else:
145 offset=[OFFSET_X_WEAK, OFFSET_Z_WEAK]
146 d2=length(x-offset)**2
147 factor_weak=exp(-2*d2/W_WEAK**2)
148 c=(C_WEAK-C)*factor_weak+C
149 friction_angle=(FRICTION_ANGLE_WEAK-FRICTION_ANGLE)*factor_weak+FRICTION_ANGLE
150 #
151 # velocity constraints:
152 #
153 x=dom.getX()
154 x0=x[0]
155 v=DIRECTION*V0/2*(1.-2*(sup(x0)-x0)/(sup(x0)-inf(x0)))*unitVector(0,DIM)
156 fixed_v_mask=(whereZero(x0-inf(x0))+whereZero(x0-sup(x0)))*unitVector(0,DIM) + \
157 whereZero(x[DIM-1])*unitVector(DIM-1,DIM)
158 if DIM==3: fixed_v_mask+=(whereZero(x[1])+whereZero(x[1]-L))*unitVector(1,DIM)
159
160 p=Scalar(0.,ReducedSolution(dom))
161 gamma_dot=sqrt(2.)*length(deviatoric(symmetric(grad(v))))
162 tau=0.
163 tau_lsup=0
164
165 flow=StokesProblemCartesian(dom)
166 flow.setTolerance(TOL/10.)
167
168 while n<1:
169
170 print "========= Time step %d ======="%( n+1,)
171 m=0
172 error=1.
173 dtau_lsup =1.
174 while dtau_lsup > TOL * tau_lsup:
175 print "--- iteration step %d ---- "%m
176 if n==1 and m==1:
177 eta_top=ETA_TOP
178 else:
179 tau_Y=clip(c*cos(friction_angle)+p*sin(friction_angle), minval=0.)
180 eta_top=clip(safeDiv(tau_Y,gamma_dot),maxval=ETA_TOP)
181 eta_eff=eta_top*(1-mask_visc) + ETA_BOTTOM*mask_visc
182 saveVTK("test.%s.vtu"%m,p=p,eta_eff=eta_eff,m=mask_visc)
183 flow.initialize(fixed_u_mask=fixed_v_mask,eta=eta_eff,f=-RHO*G*unitVector(DIM-1,DIM))
184 v,p=flow.solve(v,-3.*p,max_iter=MAX_ITER,verbose=VERBOSE,usePCG=True)
185 p*=-(1./3.)
186 gamma_dot=sqrt(2.)*length(deviatoric(symmetric(grad(v))))
187 tau, tau_old = eta_eff*gamma_dot, tau
188 dtau_lsup=Lsup(tau-tau_old)
189 tau_lsup=Lsup(tau)
190 print "increment tau = ",dtau_lsup,tau_lsup
191 m+=1
192 if m>MAX_ITER:
193 raise ValueError,"no convergence."
194 print "iteration complted after %d steps"%(m,)
195 saveVTK("test.vtu",v=v, p=p, eta=eta_eff, tau=tau);
196
197 print "Max velocity =", Lsup(v)
198 #Courant condition
199 #dt=0.4*h/(Lsup(velocity))
200 n+=1
201
202 1/0
203
204
205
206 #
207 # set up heat problem:
208 #
209 heat=TemperatureCartesian(dom,useBackwardEuler=False)
210 print "<%s> Temperature transport has been set up."%time.asctime()
211 heat.getSolverOptions().setTolerance(T_TOL)
212 heat.getSolverOptions().setVerbosity(VERBOSE)
213 fixed_T_at=whereZero(x[DIM-1])+whereZero(H-x[DIM-1])
214 heat.setInitialTemperature(clip(T,T_0))
215 heat.setValue(rhocp=1,k=1,given_T_mask=fixed_T_at)
216 #
217 # velocity constraints:
218 #
219 fixed_v_mask=Vector(0,Solution(dom))
220 faces=Scalar(0.,Solution(dom))
221 for d in range(DIM):
222 if d == DIM-1:
223 ll = H
224 else:
225 ll = L
226 if CREATE_TOPOGRAPHY and d==DIM-1:
227 fixed_v_mask+=whereZero(x[d])*unitVector(d,DIM)
228 else:
229 s=whereZero(x[d])+whereZero(x[d]-ll)
230 faces+=s
231 fixed_v_mask+=s*unitVector(d,DIM)
232 #
233 # set up velovity problem
234 #
235 flow=IncompressibleIsotropicFlowCartesian(dom, stress=stress, v=v, p=p, t=t, numMaterials=2, verbose=VERBOSE)
236 flow.setDruckerPragerLaw(tau_Y=TAU_Y/P_REF+BETA*(1.-Function(dom).getX()[DIM-1]))
237
238 flow.setElasticShearModulus(MUE)
239 flow.setTolerance(FLOW_TOL)
240 flow.setEtaTolerance(FLOW_TOL)
241 flow.setExternals(fixed_v_mask=fixed_v_mask)
242 print "<%s> Flow solver has been set up."%time.asctime()
243 #
244 # topography set-up
245 #
246 boundary=FunctionOnBoundary(dom).getX()[DIM-1]
247 top_boundary_mask=whereZero(boundary-sup(boundary))
248 surface_area=integrate(top_boundary_mask)
249 if CREATE_TOPOGRAPHY:
250 mts=Mountains(dom,eps=TOPO_SMOOTH)
251 mts.setTopography(topography)
252 print "<%s> topography has been set up."%time.asctime()
253 #
254 # let the show begin:
255 #
256 t1 = time.time()
257 print "<%s> Start time step %s (t=%s)."%(time.asctime(),n,t)
258 while t<T_END:
259 if CREATE_TOPOGRAPHY: topography_old=topography
260 v_old, p_old, stress_old=v, p, stress
261 T_old=T
262 #======= solve for velovity ====================================================================
263 eta_N=exp(Ar*((1.+V_REF*(1-Function(dom).getX()[DIM-1]))/(T_OFFSET_REF+interpolate(T,Function(dom)))-1./T_OFFSET_REF))
264 print "viscosity range :", inf(eta_N), sup(eta_N)
265 flow.setPowerLaws([eta_N, eta_N ], [ 1., TAU_0], [1,N])
266 flow.setExternals(F=Ra*T*unitVector(DIM-1,DIM))
267 # if dt<=0 or not CREATE_TOPOGRAPHY:
268 if not CREATE_TOPOGRAPHY:
269 flow.update(dt, iter_max=100, verbose=False)
270 else:
271 topography_last=topography
272 Topo_norm, error_Topo=1,1
273 i=0
274 print "DDDDD : ====",dt
275 while error_Topo > TOPO_TOL * Topo_norm:
276 flow.setStatus(t, v_old, p_old, stress_old)
277 flow.setExternals(f=-SURFACE_LOAD*(topography-dt*v)*unitVector(DIM-1,DIM)*top_boundary_mask, restoration_factor=SURFACE_LOAD*dt*top_boundary_mask)
278 flow.update(dt, iter_max=100, verbose=False)
279 v=flow.getVelocity()
280 mts.setTopography(topography_old)
281 mts.setVelocity(v)
282 topography_last, topography=topography, mts.update(dt, allow_substeps=True)
283 error_Topo=sqrt(integrate(((topography-topography_last)*top_boundary_mask)**2))
284 Topo_norm=sqrt(integrate((topography*top_boundary_mask)**2))
285 print "DDDDD :", "input=",integrate(v*top_boundary_mask)[DIM-1],"output=", integrate(topography*top_boundary_mask)/Lsup(topography), error_Topo, Topo_norm
286 print "topography update step %s error = %e, norm = %e."%(i, error_Topo, Topo_norm), Lsup(v)
287 i+=1
288 if i > TOPO_ITER_MAX:
289 raise RuntimeError,"topography did not converge after %s steps."%TOPO_ITER_MAX
290 v=flow.getVelocity()
291 for d in range(DIM):
292 print "range %d-velocity"%d,inf(v[d]),sup(v[d])
293 print "Courant = ",inf(dom.getSize()/(length(v)+1e-19)), inf(dom.getSize()**2)
294 print "<%s> flow solver completed."%time.asctime()
295 n+=1
296 t+=dt
297 # print "influx= ",integrate(inner(v,dom.getNormal())), sqrt(integrate(length(v)**2,FunctionOnBoundary(dom))), integrate(1., FunctionOnBoundary(dom))
298 print "<%s> Time step %s (t=%s) completed."%(time.asctime(),n,t)
299 #======= setup Temperature problem ====================================================================
300 #
301 heat.setValue(v=v,Q=CHI_REF*flow.getTau()**2/flow.getCurrentEtaEff())
302 dt=heat.getSafeTimeStepSize()
303 print "<%s> New time step size is %e"%(time.asctime(),dt)
304 if n == 10: 1/0
305 #======= set-up topography ==================================================================================
306 if CREATE_TOPOGRAPHY:
307 dt=min(mts.getSafeTimeStepSize()*0.5,dt)
308 print "<%s> New time step size is %e"%(time.asctime(),dt)
309 print "<%s> Start time step %s (t=%s)."%(time.asctime(),n+1,t+dt)
310 #
311 # solve temperature:
312 #
313 T=heat.getSolution(dt)
314 print "Temperature range ",inf(T),sup(T)
315 print "<%s> temperature update completed."%time.asctime()
316 #======= anaysis ==================================================================================
317 #
318 # .... Nusselt number
319 #
320 dTdz=grad(T)[DIM-1]
321 Nu=1.-integrate(v[DIM-1]*T)/integrate(dTdz)
322 eta_bar=integrate(flow.getTau())/integrate(flow.getTau()/flow.getCurrentEtaEff())
323 Ra_eff= (t_REF*RHO_0*G*H*(T_1-T_0)*ALPHA_0)/eta_bar
324 print "nusselt number = %e"%Nu
325 print "av. eta = %e"%eta_bar
326 print "effective Rayleigh number = %e"%Ra_eff
327 if CREATE_TOPOGRAPHY:
328 topo_level=integrate(topography*top_boundary_mask)/surface_area
329 valleys_deep=inf(topography)
330 mountains_heigh=sup(topography)
331 print "topography level = ",topo_level
332 print "valleys deep = ", valleys_deep
333 print "mountains_heigh = ", mountains_heigh
334 diagnostics_file.write("%e %e %e %e %e %e %e\n"%(t,Nu, topo_level, valleys_deep, mountains_heigh, eta_bar, Ra_eff))
335 else:
336 diagnostics_file.write("%e %e %e %e\n"%(t,Nu, eta_bar, Ra_eff))
337 #
338 # .... visualization
339 #
340 if t>=t_vis or n>n_vis:
341 if CREATE_TOPOGRAPHY:
342 saveVTK(os.path.join(VIS_DIR,"state.%d.vtu"%counter_vis),T=T,v=v,eta=flow.getCurrentEtaEff(), topography=topography_old, vex=mts.getVelocity())
343 else:
344 saveVTK(os.path.join(VIS_DIR,"state.%d.vtu"%counter_vis),T=T,v=v,eta=flow.getCurrentEtaEff())
345 print "<%s> Visualization file %d for time step %e generated."%(time.asctime(),counter_vis,t)
346 counter_vis+=1
347 t_vis+=DT_VIS
348 n_vis+=DN_VIS
349 # =========================
350 #
351 # create restart files:
352 #
353 if DN_RESTART>0:
354 if (n-1)%DN_RESTART == 0:
355 c=(n-1)/DN_RESTART
356 old_restart_dir="%s_%s_"%(PREFIX_RESTART,c-1)
357 new_restart_dir="%s_%s_"%(PREFIX_RESTART,c)
358 mkDir(new_restart_dir)
359 dom.MPIBarrier()
360 file(os.path.join(new_restart_dir,"stamp.%d"%dom.getMPIRank()),"w").write("%e; %e; %s; %s; %s; %e;\n"%(t, t_vis, n_vis, n, counter_vis, dt))
361 v.dump(os.path.join(new_restart_dir,"v.nc"))
362 p.dump(os.path.join(new_restart_dir,"p.nc"))
363 T.dump(os.path.join(new_restart_dir,"T.nc"))
364 if CREATE_TOPOGRAPHY: topography.dump(os.path.join(new_restart_dir,"topo.nc"))
365 removeRestartDirectory(old_restart_dir)
366 print "<%s> Restart files written to %s."%(time.asctime(),new_restart_dir)
367 print "<%s> Calculation finalized after %s sec."%(time.asctime(),time.time() - t1)
368

  ViewVC Help
Powered by ViewVC 1.1.26