/[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 4657 - (show annotations)
Thu Feb 6 06:12:20 2014 UTC (5 years, 4 months ago) by jfenwick
File MIME type: text/x-python
File size: 14514 byte(s)
I changed some files.
Updated copyright notices, added GeoComp.



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

  ViewVC Help
Powered by ViewVC 1.1.26