diff --git a/gusto/core/kernels.py b/gusto/core/kernels.py index 422b7456..44ad8457 100644 --- a/gusto/core/kernels.py +++ b/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 = (""" + min_value1 = 0.0 + min_value2 = 0.0 min_value = 0.0 + 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)}) diff --git a/gusto/spatial_methods/augmentation.py b/gusto/spatial_methods/augmentation.py index e271061f..a301fd30 100644 --- a/gusto/spatial_methods/augmentation.py +++ b/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,8 +473,10 @@ 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): @@ -478,7 +484,7 @@ 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') diff --git a/gusto/spatial_methods/limiters.py b/gusto/spatial_methods/limiters.py index c0431154..823bb30b 100644 --- a/gusto/spatial_methods/limiters.py +++ b/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)