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

  ViewVC Help
Powered by ViewVC 1.1.26