/[escript]/trunk/escript/py_src/mountains.py
ViewVC logotype

Annotation of /trunk/escript/py_src/mountains.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2210 - (hide annotations)
Tue Jan 13 03:49:47 2009 UTC (10 years, 8 months ago) by artak
File MIME type: text/x-python
File size: 3605 byte(s)
First version of the mountains class.
1 artak 2210
2     ########################################################
3     #
4     # Copyright (c) 2003-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) 2003-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__="http://www.uq.edu.au/esscc/escript-finley"
21    
22     from esys.escript import *
23     from esys.escript.linearPDEs import LinearPDE
24     from esys import finley
25     import sys
26    
27     class Mountains:
28     """
29     The Mountains class is defined by the following equations:
30    
31     (1) w_i,aa+(1/eps)*w_i,33=0 where 0<eps<<1 and a=1,2 and w_i is the extension of the surface velocity where w_i(x_3=1)=v_i.
32    
33     (2) Integration of topography PDE using Taylor-Galerkin upwinding to stabilize the advection terms
34     H^(t+dt)=H^t+dt*w_3+w_hat*dt*[(div(w_hat*H^t)+w_3)+(dt/2)+H^t],
35     where w_hat=w*[1,1,0], dt<0.5*d/max(w_i), d is a characteristic element size; H(x_3=1)=lambda (?)
36    
37    
38     """
39     def __init__(self,domain,v,eps=0.01,z=1):
40     """
41     Sets up the level set method.
42    
43     @param domain: the domain where the mountains is used
44     @param eps: the parameter for (1)
45     @param z: the height of the box
46     """
47     self.__domain = domain
48     self.__eps = eps
49     self.__PDE_W = LinearPDE(domain)
50     self.__PDE_W.setSymmetryOn()
51     self.__PDE_H = LinearPDE(domain)
52     self.__PDE_H.setSymmetryOn()
53     self.__PDE_H.setValue(D=1.0)
54     self.__PDE_H.setSolverMethod(solver=LinearPDE.LUMPING)
55     self.__h = inf(domain.getSize())
56     self.__v=v
57     self.__z=z
58     self.__x=domain.getX()
59     self.__DIM=domain.getDim()
60     self.__fixed_w_mask=whereZero(self.__z-self.__x[self.__DIM-1])+whereZero(self.__x[self.__DIM-1])
61     self.__fixed_H_mask=whereZero(self.__x[self.__DIM-1])
62     self.__A=kronecker(domain)
63     self.__A[self.__DIM-1,self.__DIM-1]=1/self.__eps
64     self.__PDE_W.setValue(A=self.__A, q=self.__fixed_w_mask)
65     self.__PDE_H.setValue(q=self.__fixed_H_mask)
66     self.__H_t=Scalar(0.0, Solution(domain))
67     self.__dt=0.
68    
69     def update(self,u=None,H_t=None,dt=None):
70     """
71     Sets a new W and updates the H function.
72    
73     @param dt: time step forward
74     """
75     if H_t==None:
76     H_t=self.__H_t
77    
78     if u==None:
79     u=self.__v
80    
81     for d in range(self.__DIM):
82     self.__PDE_W.setValue(r=u[d]*self.__x[self.__DIM-1]/self.__z)
83     u[d]=self.__PDE_W.getSolution()
84    
85     if dt==None:
86     dt=0.5*self.__h/sup(u)
87     self.__dt=dt
88     else:
89     self.__dt=dt
90    
91     w=u
92     w_tilda=1.*w
93     w_tilda[self.__DIM-1]=0
94    
95     self.__PDE_H.setValue(X=((div(w_tilda*H_t)+w[self.__DIM-1])*dt/2+H_t)*w_tilda*dt, Y=w[self.__DIM-1]*dt+H_t)
96     self.__H_t=self.__PDE_H.getSolution()
97    
98     self.__v=w
99    
100     return w,self.__H_t
101    
102     def getH(self):
103     """
104     Returns the mesh size.
105     """
106     return self.__h
107    
108     def getDt(self):
109     """
110     Returns the mesh size.
111     """
112     return self.__dt
113    
114    
115     def getDomain(self):
116     """
117     Returns the domain.
118     """
119     return self.__domain
120    
121     def setTolerance(self,tolerance=1e-3):
122     self.__PDE_W.setTolerance(tolerance)
123     self.__PDE_H.setTolerance(tolerance)
124    
125    

  ViewVC Help
Powered by ViewVC 1.1.26