correcting mean mixing ratio limiter

gusto/core/kernels.py

@@ -10,7 +10,7 @@
 """
 
 from firedrake import dx
-from firedrake.parloops import par_loop, READ, WRITE, MIN, MAX, op2
+from firedrake.parloops import par_loop, READ, WRITE, INC, MIN, MAX, op2
 import numpy as np
 
 
@@ -132,19 +132,20 @@ def __init__(self, V):
         domain = "{{[i]: 0 <= i < {nDOFs_base}}}".format(**shapes)
 
         instrs = ("""
+                  <float64> min_value1 = 0.0
+                  <float64> min_value2 = 0.0
                   <float64> min_value = 0.0
+                  <float64> new_lamda = 0.0
                   for i
-                      min_value = fmin(mX_field[i*4], mX_field[i*4+1])
-                      min_value = fmin(min_value, mX_field[i*4+2])
-                      min_value = fmin(min_value, mX_field[i*4+3])
+                      min_value1 = fmin(mX_field[i*4], mX_field[i*4+1])
+                      min_value2 = fmin(mX_field[i*4+2], mX_field[i*4+3])
+                      min_value = fmin(min_value1, min_value2)
                       if min_value < 0.0
-                          lamda[i] = -min_value/(mean_field[i] - min_value)
-                      else
-                          lamda[i] = 0.0
+                          lamda[i] = fmax(lamda[i],-min_value/(mean_field[i] - min_value))
                       end
                   end
                   """)
-
+        #lamda[i] = fmax(lamda[i],-min_value/(mean_field[i] - min_value))
         self._kernel = (domain, instrs)
 
     def apply(self, lamda, mX_field, mean_field):
@@ -156,7 +157,7 @@ def apply(self, lamda, mX_field, mean_field):
                 lives in the continuous target space.
         """
         par_loop(self._kernel, dx,
-                 {"lamda": (lamda, WRITE),
+                 {"lamda": (lamda, INC),
                   "mX_field": (mX_field, READ),
                   "mean_field": (mean_field, READ)})
 

gusto/spatial_methods/augmentation.py

@@ -307,11 +307,12 @@ def __init__(
             self.mX_idxs.append(mX_idx)
 
             # Define a limiter
-            self.limiters.append(MeanLimiter(eqns.spaces[mX_idx]))
+            #self.limiters.append(MeanLimiter(eqns.spaces[mX_idx]))
             #self.sublimiters.update({mX_name: MeanLimiter(eqns.spaces[mX_idx])})
 
         #self.limiter = MixedFSLimiter(self.eqn_orig, self.sublimiters)
         #self.limiter = MixedFSLimiter(sublimiters)
+        self.limiters = MeanLimiter(mX_spaces)
 
         # Contruct projectors for computing the mean mixing ratios
         self.mX_ins = [Function(mX_spaces[i]) for i in range(self.mX_num)]
@@ -451,7 +452,10 @@ def pre_apply(self, x_in):
 
         #for idx, field in enumerate(self.eqn.field_names):
         #    self.x_in.subfunctions[idx].assign(x_in.subfunctions[idx])
+        print('in pre-apply')
         for idx in range(self.idx_orig):
+            print(np.min(x_in.subfunctions[idx].dat.data))
+            print('\n')
             self.x_in.subfunctions[idx].assign(x_in.subfunctions[idx])
 
         # Set the mean mixing ratio to be zero, just because
@@ -469,16 +473,18 @@ def post_apply(self, x_out):
         Args:
             x_out (:class:`Function`): The output fields
         """
-
+        print('in post apply')
         for idx in range(self.idx_orig):
+            print(np.min(self.x_out.subfunctions[idx].dat.data))
+            print('\n')
             x_out.subfunctions[idx].assign(self.x_out.subfunctions[idx])
 
     def update(self, x_in_mixed):
         """
         Compute the mean mixing ratio field by projecting the mixing 
         ratio from DG1 into DG0.
 
-        SHouldn't this be a conservative projection??!!
+        To DO: Shouldn't this be a conservative projection??!!
         This requires density fields...
 
         Args:
@@ -487,10 +493,37 @@ def update(self, x_in_mixed):
         print('in update')
         for i in range(self.mX_num):
             self.mX_ins[i].assign(x_in_mixed.subfunctions[self.mX_idxs[i]])
+            print('\n min of mX field:')
+            print(np.min(self.mX_ins[i].dat.data))
             self.compute_means[i].project()
             self.x_in.subfunctions[self.mean_idxs[i]].assign(self.mean_outs[i])
 
     def limit(self, x_in_mixed):
+        # Ensure non-negativity by applying the blended limiter
+        mX_pre = []
+        means = []
+        print('limiting within the augmentation')
+        for i in range(self.mX_num):
+            mX_pre.append(x_in_mixed.subfunctions[self.mX_idxs[i]])
+            means.append(x_in_mixed.subfunctions[self.mean_idxs[i]])
+
+            print('\n min of mX field:')
+            print(np.min(mX_pre[i].dat.data))
+            print('\n min of mean field:')
+            print(np.min(means[i].dat.data))
+
+        self.limiters.apply(mX_pre, means)
+
+        print('\n After applying blended limiter')
+
+        for i in range(self.mX_num):
+            x_in_mixed.subfunctions[self.mX_idxs[i]].assign(mX_pre[i])
+            print('\n min of mX field:')
+            print(np.min(mX_pre[i].dat.data))
+            print('\n min of mean field:')
+            print(np.min(means[i].dat.data))
+
+    def limit_old(self, x_in_mixed):
         # Ensure non-negativity by applying the blended limiter
         for i in range(self.mX_num):
             print('limiting within the augmentation')

gusto/spatial_methods/limiters.py

@@ -269,7 +269,7 @@ class MeanLimiter(object):
     factor is given by the DG0 function lamda.
     """
 
-    def __init__(self, space):
+    def __init__(self, spaces):
         """
         Args:
             space: The mixed function space for the equation set
@@ -279,16 +279,19 @@ def __init__(self, space):
             ValueError: If the space is not appropriate for the limiter.
         """
 
-        self.space = space 
-        mesh = space.mesh()
+        #self.space = space 
+        #mesh = space.mesh()
 
         # check that space is DG1
-        degree = space.ufl_element().degree()
-        if (space.ufl_element().sobolev_space.name != 'L2'
-            or ((type(degree) is tuple and np.any([deg != 1 for deg in degree]))
-                and degree != 1)):
-            raise NotImplementedError('MeanLimiter only implemented for mixing' +
-                                      'ratios in the DG1 space')
+        for space in spaces:
+            degree = space.ufl_element().degree()
+            if (space.ufl_element().sobolev_space.name != 'L2'
+                or ((type(degree) is tuple and np.any([deg != 1 for deg in degree]))
+                    and degree != 1)):
+                raise NotImplementedError('MeanLimiter only implemented for mixing' +
+                                          'ratios in the DG1 space')
+        self.space = spaces[0] 
+        mesh = self.space.mesh()
 
         # Create equispaced DG1 space needed for limiting
         if space.extruded:
@@ -301,17 +304,22 @@ def __init__(self, space):
             DG1_element = FiniteElement("DG", cell, 1, variant="equispaced")
 
         DG1_equispaced = FunctionSpace(mesh, DG1_element)
-        print(DG1_equispaced.finat_element.space_dimension())
-
         DG0 = FunctionSpace(mesh, 'DG', 0)
 
-        self.lamda = Function(DG1_equispaced)
+        self.lamda = Function(DG0)
+        #self.minus_lamda = Function(DG0)
         self.mX_field = Function(DG1_equispaced)
         self.mean_field = Function(DG0)
 
+        print(DG1_equispaced.finat_element.space_dimension())
+        print(DG0.finat_element.space_dimension())
+        print(self.space.finat_element.space_dimension())
+        print(self.space.mesh().topological_dimension())
+
         self._kernel = MeanMixingRatioWeights(self.space)
+        self._clip_zero_kernel = ClipZero(self.space)
 
-    def apply(self, mX_field, mean_field):
+    def apply(self, mX_fields, mean_fields):
         """
         The application of the limiter to the field.
 
@@ -322,14 +330,34 @@ def apply(self, mX_field, mean_field):
              AssertionError: If the field is not in the correct space.
          """
 
-        #self.field_old.interpolate(mX_field)
-        #self.mean_field.interpolate(mean_field)
+        # Set the weights, lamda, to zero
+        self.lamda.interpolate(Constant(0.0))
+
+        for i in range(len(mX_fields)):
+            # Update the weights, lamda
+            self._kernel.apply(self.lamda, mX_fields[i], mean_fields[i])
+            #print(self.lamda.dat.data)
+
+        #self.minus_lamda.interpolate(Constant(1.0) - self.lamda)
+
+        #As a hack for now, clip zero when required
+
+        for i in range(len(mX_fields)):
+
+            #mX_fields[i].interpolate(self.minus_lamda*mX_fields[i] + self.lamda*mean_fields[i])
+            mX_fields[i].interpolate((Constant(1.0) - self.lamda)*mX_fields[i] + self.lamda*mean_fields[i])
+            #mX_fields[i].interpolate(self.one_m_lamda*mX_fields[i] + self.lamda*mean_fields[i])
+            #mX_fields[i].interpolate(self.lamda*mean_fields[i])
 
+            #As a hack for now, clip zero when required
+            self._clip_zero_kernel.apply(mX_fields[i], mX_fields[i])
+
+
         # Compute the weights, lamda:
-        self._kernel.apply(self.lamda, mX_field, mean_field)
+        #self._kernel.apply(self.lamda, mX_field, mean_field)
 
-        print('applying limiter')
+        #print('applying limiter')
         #print(self.lamda.dat.data)
 
         # Compute the blended field
-        mX_field.interpolate((Constant(1.0) - self.lamda)*mX_field + self.lamda*mean_field)
+        #mX_field.interpolate((Constant(1.0) - self.lamda)*mX_field + self.lamda*mean_field)