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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3981 - (show annotations)
Fri Sep 21 02:47:54 2012 UTC (7 years, 2 months ago) by jfenwick
File MIME type: text/x-python
File size: 2557 byte(s)
First pass of updating copyright notices
1 ##############################################################################
2 #
3 # Copyright (c) 2003-2012 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 __copyright__="""Copyright (c) 2003-2012 by University of Queensland
16 http://www.uq.edu.au
17 Primary Business: Queensland, Australia"""
18 __license__="""Licensed under the Open Software License version 3.0
19 http://www.opensource.org/licenses/osl-3.0.php"""
20 __url__="https://launchpad.net/escript-finley"
21
22 # $Id$
23
24 from esys.escript import *
25 from esys.escript import unitsSI as U
26 from esys.escript.pdetools import MaskFromBoundaryTag
27 from esys.finley import ReadMesh
28 from esys.weipa import saveVTK
29 from esys.escript.models import StokesProblemCartesian
30 from math import ceil
31 #
32 # Parameter
33 #
34 DIM=2
35 MESHFILE="sub.fly"
36 ETA=1.e22*U.Pa*U.sec
37 V_MAX=1.*U.cm/U.yr
38 ALPHA=30*U.DEG
39 STRIKE=10*U.DEG
40 DIP=30*U.DEG
41 N=1 # boudary layer control
42
43
44 g=9.81*U.m/U.sec**2
45 #
46 # derived values
47 #
48 dom=ReadMesh(MESHFILE)
49 DIM=dom.getDim()
50 bb=boundingBox(dom)
51 LX=bb[0][1]-bb[0][0]
52 if DIM == 3: LY=bb[1][1]-bb[1][0]
53 DEPTH=bb[DIM-1][1]-bb[DIM-1][0]
54
55 sc=StokesProblemCartesian(dom)
56 x = dom.getX()
57 #
58 v=Vector(0.,Solution(dom))
59 mask=Vector(0.,Solution(dom))
60 #
61 # in subduction zone:
62 #
63
64 if DIM==2:
65 S=numpy.array([0.,0.])
66 X0=[bb[0][0],bb[1][1]]
67 dd=[-cos(ALPHA),-sin(ALPHA)]
68 else:
69 S=numpy.array([sin(STRIKE),cos(STRIKE),0.])
70 X0=[bb[0][0],bb[1][0],bb[2][1]]
71 dd=[-cos(ALPHA),0.,-sin(ALPHA)]
72 r=sqrt(length(x-X0)**2-inner(X0-x,S)**2)
73 v=V_MAX*r*dd
74 mask=MaskFromBoundaryTag(dom,"subduction")*[ 1. for i in range(DIM) ]
75 #
76 # back of the domain
77 #
78 v=v*(1.-whereZero(x[0]-bb[0][1])*kronecker(DIM)[0])
79 mask+=whereZero(x[0]-bb[0][1])*kronecker(DIM)[0]
80 #
81 # bottom of the domain
82 #
83 v=v*(1.-((bb[DIM-1][1]-x[DIM-1])/DEPTH)**N*kronecker(DIM)[DIM-1])
84 mask+=whereZero(x[DIM-1]-bb[DIM-1][0])*kronecker(DIM)[DIM-1]
85 #
86 # faces of the domain:
87 #
88 if DIM==3:
89 v=v*(1.-(((x[1]-bb[1][0])*(bb[1][1]-x[1])/(0.5*LY)**2)**N*kronecker(DIM)[1]))
90 mask+=(whereZero(x[1]-bb[1][0])+whereZero(x[1]-bb[1][1]))*kronecker(DIM)[1]
91 sc.initialize(eta=ETA, fixed_u_mask= mask)
92 p=Scalar(0.,ReducedSolution(dom))
93 v,p=sc.solve(v,p, verbose=True)
94 saveVTK("u.vtu",velocity=v,pressure=p,m=mask)
95

  ViewVC Help
Powered by ViewVC 1.1.26