/[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 6651 - (show annotations)
Wed Feb 7 02:12:08 2018 UTC (20 months, 1 week ago) by jfenwick
File MIME type: text/x-python
File size: 14552 byte(s)
Make everyone sad by touching all the files

Copyright dates update

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

  ViewVC Help
Powered by ViewVC 1.1.26