/[escript]/trunk/doc/examples/cookbook/example10d.py
ViewVC logotype

Annotation of /trunk/doc/examples/cookbook/example10d.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3892 - (hide annotations)
Tue Apr 10 08:57:23 2012 UTC (6 years, 10 months ago) by jfenwick
File MIME type: text/x-python
File size: 4861 byte(s)
Merged changes across from the attempt2 branch.
This version builds and passes python2 tests.
It also passes most python3 tests.



1 ahallam 3405
2     ########################################################
3     #
4     # Copyright (c) 2009 by University of Queensland
5     # Earth Systems Science Computational Center (ESSCC)
6     # http://www.uq.edu.au/esscc
7     #
8     # Primary Business: Queensland, Australia
9     # Licensed under the Open Software License version 3.0
10     # http://www.opensource.org/licenses/osl-3.0.php
11     #
12     ########################################################
13    
14     __copyright__="""Copyright (c) 2009 by University of Queensland
15     Earth Systems Science Computational Center (ESSCC)
16     http://www.uq.edu.au/esscc
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     """
23     Author: Antony Hallam antony.hallam@uqconnect.edu.au
24     """
25    
26     ############################################################FILE HEADER
27     # example10a.py
28     # Model of gravitational Potential for an infinite strike horizontal
29     # cylinder. This example tests the boundary condition criterion.
30    
31     #######################################################EXTERNAL MODULES
32     # To solve the problem it is necessary to import the modules we require.
33     from esys.escript import * # This imports everything from the escript library
34     from esys.escript.unitsSI import *
35     from esys.escript.linearPDEs import LinearPDE # This defines LinearPDE as LinearPDE
36     from esys.finley import Rectangle # This imports the rectangle domain function from finley
37     from esys.weipa import saveVTK # This imports the VTK file saver from weipa
38     import os, sys #This package is necessary to handle saving our data.
39     from math import pi, sqrt, sin, cos
40    
41     from esys.escript.pdetools import Projector, Locator
42     from cblib import toRegGrid
43     from esys.finley import ReadGmsh
44    
45     import matplotlib
46     matplotlib.use('agg') #It's just here for automated testing
47    
48     import pylab as pl #Plotting package
49     import numpy as np
50    
51     ########################################################MPI WORLD CHECK
52     if getMPISizeWorld() > 1:
53     import sys
54 jfenwick 3892 print("This example will not run in an MPI world.")
55 ahallam 3405 sys.exit(0)
56    
57     #################################################ESTABLISHING VARIABLES
58     #PDE related
59     rho=200.0 #the density contrast of the source
60     G=6.67300*10E-11 #gravitational constant
61     R=100. #radius of the source body
62     ################################################ESTABLISHING PARAMETERS
63     #the folder to put our outputs in, leave blank "" for script path
64     save_path= os.path.join("data","example10")
65     mesh_path = "data/example10m"
66     #ensure the dir exists
67     mkDir(save_path)
68    
69     ####################################################DOMAIN CONSTRUCTION
70     #Geometric and material property related variables.
71     domain=ReadGmsh(os.path.join(mesh_path,'example10m.msh'),2) # create the domain
72     x=domain.getX()
73     #Domain related dimensions from the imported file.
74     mx = Lsup(x)*m #meters - model length
75     my = Lsup(x)*m #meters - model width
76     rholoc=[mx/2,my/2]
77    
78     mask=wherePositive(R-length(x-rholoc))
79     rhoe=rho*mask
80     kro=kronecker(domain)
81    
82 ahallam 3409 #Define the analytic solution.
83 ahallam 3405 def analytic_gz(x,z,R,drho):
84     G=6.67300*10E-11
85     return G*2*np.pi*R*R*drho*(z/(x*x+z*z))
86    
87 ahallam 3409 #Define the boundary conditions.
88 ahallam 3405 q=whereZero(x[1]-my)+whereZero(x[1])+whereZero(x[0])+whereZero(x[0]-mx)
89     ###############################################ESCRIPT PDE CONSTRUCTION
90     mypde=LinearPDE(domain)
91     mypde.setValue(A=kro,Y=4.*3.1415*G*rhoe,q=q,r=0.)
92     mypde.setSymmetryOn()
93     sol=mypde.getSolution()
94    
95 ahallam 3409 g_field=grad(sol) #The gravitational acceleration g.
96 ahallam 3405 g_fieldz=g_field*[0,1] #The vertical component of the g field.
97     gz=length(g_fieldz) #The magnitude of the vertical component.
98     # Save the output to file.
99     saveVTK(os.path.join(save_path,"ex10d.vtu"),\
100     grav_pot=sol,g_field=g_field,g_fieldz=g_fieldz,gz=gz,rho=rhoe)
101    
102     ################################################MODEL SIZE SAMPLING
103     smoother = Projector(domain) #Function smoother.
104     z=1000. #Distance of profile from source.
105     sol_angz=[] #Array for analytic gz
106     sol_anx=[] #Array for x
107     #calculate analytic gz and x location.
108     for x in range(int(-mx/2.),int(mx/2.),10):
109     sol_angz.append(analytic_gz(x,z,R,rho))
110     sol_anx.append(x+mx/2)
111     #save analytic solution
112     pl.savetxt(os.path.join(save_path,"ex10d_as.asc"),sol_angz)
113    
114     sol_escgz=[] #Array for escript solution for gz
115     #calculate the location of the profile in the domain
116     for i in range(0,len(sol_anx)):
117     sol_escgz.append([sol_anx[i],my/2.-z])
118    
119     rec=Locator(gz.getFunctionSpace(),sol_escgz) #location to record
120     psol=rec.getValue(smoother(gz)) #extract the solution
121     #save X and escript solution profile to file
122     pl.savetxt(os.path.join(save_path,"ex10d_%05d.asc"%mx),psol)
123     pl.savetxt(os.path.join(save_path,"ex10d_x%05d.asc"%mx),sol_anx)
124    
125     #plot the pofile and analytic solution using example10q.py
126    

  ViewVC Help
Powered by ViewVC 1.1.26