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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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.__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.  
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

  ViewVC Help
Powered by ViewVC 1.1.26