# Diff of /trunk/escript/py_src/mountains.py

revision 2562 by gross, Tue Jun 30 05:49:22 2009 UTC revision 2563 by gross, Tue Jul 28 03:50:45 2009 UTC
# Line 34  class Mountains: Line 34  class Mountains:
34        H^(t+dt)=H^t+dt*w_3+w_hat*dt*[(div(w_hat*H^t)+w_3)+(dt/2)+H^t],        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 (?)        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    def __init__(self,domain,v,eps=0.01,z=1):    def __init__(self,domain,eps=0.01):
39      """      """
40      Sets up the level set method.      Sets up the level set method.
41
42      @param domain: the domain where the mountains is used      @param domain: the domain where the mountains is used
43      @param eps: the parameter for (1)      @param eps: the smoothing parameter for (1)
@param z: the height of the box
44      """      """
45        if eps<=0:
46            raise ValueError("Smmoting parameter eps must be positive.")
47      self.__domain = domain      self.__domain = domain
48      self.__eps = eps      self.__eps = eps
49        self.__DIM=domain.getDim()
50        z=domain.getX()[self.__DIM-1]
51
52      self.__PDE_W = LinearPDE(domain)      self.__PDE_W = LinearPDE(domain)
53      self.__PDE_W.setSymmetryOn()      self.__PDE_W.setSymmetryOn()
54        A=kronecker(domain)
55        A[self.__DIM-1,self.__DIM-1]=1/self.__eps
56        self.__PDE_W.setValue(A=A, q=whereZero(sup(z)-z)+whereZero(inf(z)-z))
57
58      self.__PDE_H = LinearPDE(domain)      self.__PDE_H = LinearPDE(domain)
59      self.__PDE_H.setSymmetryOn()      self.__PDE_H.setSymmetryOn()
60      self.__PDE_H.setValue(D=1.0)      self.__PDE_H.setValue(D=1.0)
61      self.__PDE_H.getSolverOptions().setSolverMethod(SolverOptions.LUMPING)      self.__PDE_H.getSolverOptions().setSolverMethod(SolverOptions.LUMPING)
62      self.__h = inf(domain.getSize())      self.__PDE_H.setValue(q=whereZero(inf(z)-z))
self.__v=v
self.__z=z
self.__x=domain.getX()
self.__DIM=domain.getDim()
self.__A=kronecker(domain)
self.__A[self.__DIM-1,self.__DIM-1]=1/self.__eps
self.__H_t=Scalar(0.0, Solution(domain))
self.__dt=0.
63
64    def update(self,u=None,H_t=None,dt=None):      self.setVelocity()
65        self.setTopography()
66      def getSolverOptionsForSmooting(self):
67         """
68         returns the solver options for the smoothing/extrapolation
69         """
70         return self.__PDE_W.getSolverOptions()
71
72      def getSolverOptionsForUpdate(self):
73         """
74         returns the solver options for the topograthy update
75         """
76         return self.__PDE_H.getSolverOptions()
77      def getDomain(self):
78        """        """
79        Sets a new W and updates the H function.        Returns the domain.

@param dt: time step forward
80        """        """
81        if H_t==None:        return self.__domain
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
82
83        return w,self.__H_t    def setVelocity(self,v=None):
84          """
85          set a new velocity. v is define on the entire domain but only the surface values are used.
86
87    def getH(self):        @param v: velocity field. If None zero is used.
88          @type v: vector
89        """        """
90        Returns the mesh size.        self.__dt=None
91          self.__v=Vector(0.,Solution(self.getDomain()))
92          if not v == None:
93            z=self.getDomain().getX()[self.__DIM-1]
94            z_min=inf(z)
95            z_max=sup(z)
96            f=(z-z_min)/(z_max-z_min)
97            for d in range(self.__DIM):
98               self.__PDE_W.setValue(r=v[d]*f)
99               self.__v[d]=self.__PDE_W.getSolution()
100      def getVelocity(self):
101        """        """
102        return self.__h        returns the smoothed/extrapolated velocity
103          @rtype: vector L{Data}
104          """
105          return self.__v
106
107      def setTopography(self,H=None):
108        """
109        set the topography to H where H defines the vertical displacement. H is defined for the entire domain.
110
111        @param H: the topography.  If None zero is used.
112        @type H: scalar
113        """
114
115    def getDt(self):      if H==None:
116           self.__H=Scalar(0.0, Solution(self.getDomain()))
117        else:
118           self.__H=interpolate(H, Solution(self.getDomain()))
119
120      def getTopography(self):
121         """
122         returns the current topography.
123         @rtype: scalar L{Data}
124         """
125         return self.__H
126
127      def getSafeTimeStepSize(self):
128        """        """
129        Returns the time step value.        Returns the time step value.
130
131          @rtype: C{float}
132        """        """
133          if self.__dt == None:
134               self.__dt=0.5*inf(self.getDomain().getSize()/length(self.getVelocity()))
135        return self.__dt        return self.__dt
136          def update(self,dt=None):

def getDomain(self):
137        """        """
138        Returns the domain.        Sets a new W and updates the H function.
139
140          @param dt: time step forward. If None the save time step size is used.
141          @type dt: positve C{float} which is less or equal than the safe time step size.
142
143        """        """
144        return self.__domain        if dt == None:
145                dt = self.getSafeTimeStepSize()
146          if dt<=0:
147               raise ValueError("Time step size must be positive.")
148          if dt>self.getSafeTimeStepSize():
149               raise ValueError("Time step must be less than safe time step size = %e."%self.getSafeTimeStepSize())
150
151          H=self.getTopography()
152          w=self.getVelocity()
153          w_tilda=1.*w
154          w_tilda[self.__DIM-1]=0
155          w_z=w[self.__DIM-1]
156
157          self.__PDE_H.setValue(X=((div(w_tilda*H)+w_z)*dt/2+H)*w_tilda*dt, Y=w_z*dt+H, r=H)
158          self.setTopography(self.__PDE_H.getSolution())
159
160          return self.getTopography()
161
def setTolerance(self,tolerance=1e-3):
self.__PDE_W.getSolverOptions().setTolerance(tolerance)
self.__PDE_H.getSolverOptions().setTolerance(tolerance)

Legend:
 Removed from v.2562 changed lines Added in v.2563