escript/py_src/mountains.py


########################################################
#
# Copyright (c) 2003-2009 by University of Queensland
# Earth Systems Science Computational Center (ESSCC)
# http://www.uq.edu.au/esscc
#
# Primary Business: Queensland, Australia
# Licensed under the Open Software License version 3.0
# http://www.opensource.org/licenses/osl-3.0.php
#
########################################################

__copyright__="""Copyright (c) 2003-2009 by University of Queensland
Earth Systems Science Computational Center (ESSCC)
http://www.uq.edu.au/esscc
Primary Business: Queensland, Australia"""
__license__="""Licensed under the Open Software License version 3.0
http://www.opensource.org/licenses/osl-3.0.php"""
__url__="http://www.uq.edu.au/esscc/escript-finley"

from esys.escript import *
from esys.escript.linearPDEs import LinearPDE
from esys import finley
import sys

class Mountains:
  """
  The Mountains class is defined by the following equations:
  
  (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.
  
  (2) Integration of topography PDE using Taylor-Galerkin upwinding to stabilize the advection terms
      H^(t+dt)=H^t+dt*w_3+w_hat*dt*[(div(w_hat*H^t)+w_3)+(dt/2)+H^t],
      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 (?) 
      

  """
  def __init__(self,domain,v,eps=0.01,z=1):
    """
    Sets up the level set method.

    @param domain: the domain where the mountains is used
    @param eps: the parameter for (1)
    @param z: the height of the box
    """
    self.__domain = domain
    self.__eps = eps
    self.__PDE_W = LinearPDE(domain)
    self.__PDE_W.setSymmetryOn()
    self.__PDE_H = LinearPDE(domain)
    self.__PDE_H.setSymmetryOn()
    self.__PDE_H.setValue(D=1.0)
    self.__PDE_H.setSolverMethod(solver=LinearPDE.LUMPING)
    self.__h = inf(domain.getSize())
    self.__v=v
    self.__z=z
    self.__x=domain.getX()
    self.__DIM=domain.getDim()
    self.__fixed_w_mask=whereZero(self.__z-self.__x[self.__DIM-1])+whereZero(self.__x[self.__DIM-1])
    self.__fixed_H_mask=whereZero(self.__x[self.__DIM-1])
    self.__A=kronecker(domain)
    self.__A[self.__DIM-1,self.__DIM-1]=1/self.__eps
    self.__PDE_W.setValue(A=self.__A, q=self.__fixed_w_mask)
    self.__PDE_H.setValue(q=self.__fixed_H_mask)
    self.__H_t=Scalar(0.0, Solution(domain))
    self.__dt=0.

  def update(self,u=None,H_t=None,dt=None):
      """
      Sets a new W and updates the H function.

      @param dt: time step forward
      """
      if H_t==None:
        H_t=self.__H_t
        
      if u==None:
         u=self.__v
         
      for d in range(self.__DIM):
        self.__PDE_W.setValue(r=u[d]*self.__x[self.__DIM-1]/self.__z)
        u[d]=self.__PDE_W.getSolution()

      if dt==None:
        dt=0.5*self.__h/sup(u)
        self.__dt=dt
      else:
        self.__dt=dt
 
      w=u
      w_tilda=1.*w
      w_tilda[self.__DIM-1]=0

      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)
      self.__H_t=self.__PDE_H.getSolution()
      
      self.__v=w

      return w,self.__H_t

  def getH(self):
      """
      Returns the mesh size.
      """
      return self.__h

  def getDt(self):
      """
      Returns the mesh size.
      """
      return self.__dt
    

  def getDomain(self):
      """
      Returns the domain.
      """
      return self.__domain

  def setTolerance(self,tolerance=1e-3):
    self.__PDE_W.setTolerance(tolerance)
    self.__PDE_H.setTolerance(tolerance)


1
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+dtw_3+w_hatdt[(div(w_hatH^t)+w_3)+(dt/2)+H^t],
35	where w_hat=w[1,1,0], dt<0.5d/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_tildaH_t)+w[self.__DIM-1])dt/2+H_t)w_tildadt, 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