| 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 |
|
Parent Directory
| 
Revision Log