From 72b74d33b3f8cb60f96e21b39824b71a9434feaa Mon Sep 17 00:00:00 2001 From: jhcollins <37124212+jhcollins@users.noreply.github.com> Date: Tue, 28 Aug 2018 11:24:42 -0500 Subject: [PATCH] CLSVOF lambda fix, scales properly with domain size (#790) * CLSVOF lambda fix, scales properly with domain size * change name of parameter in pyx file * Fix CLSVOF pure level set tests * fix tests of CLSVOF with RANS3PF * Fix the test for CLSVOF with RANS2P * Made the changes to lambda term in CLSVOF.h * added both lambda scaling methods for use in testing; lambda_scaling 0 is default with abs(), lambda_scaling 1 is (1-dirac) * Add more changes about scaling * Fix tests for CLSVOF --- proteus/mprans/CLSVOF.h | 81 ++++++++++++++++--- proteus/mprans/CLSVOF.py | 25 +++++- proteus/mprans/cCLSVOF.pyx | 15 +++- proteus/tests/CLSVOF/pure_level_set/clsvof.py | 1 + 4 files changed, 106 insertions(+), 16 deletions(-) diff --git a/proteus/mprans/CLSVOF.h b/proteus/mprans/CLSVOF.h index 78ecbff860..82c14a6031 100644 --- a/proteus/mprans/CLSVOF.h +++ b/proteus/mprans/CLSVOF.h @@ -8,6 +8,7 @@ // True characteristic functions #define heaviside(z) (z>0 ? 1. : (z<0 ? 0. : 0.5)) #define sign(z) (z>0 ? 1. : (z<0 ? -1. : 0.)) +#define LAMBDA_SCALING 0 namespace proteus { @@ -92,7 +93,9 @@ namespace proteus double* min_distance, double* max_distance, double* mean_distance, + double* volume_domain, double norm_factor_lagged, + double VelMax, // normal reconstruction double* projected_qx_tn, double* projected_qy_tn, @@ -165,7 +168,8 @@ namespace proteus double epsFactDirac, double lambdaFact, // normalization factor - double norm_factor_lagged)=0; + double norm_factor_lagged, + double VelMax)=0; virtual void calculateMetricsAtEOS( //EOS=End Of Simulation double* mesh_trial_ref, double* mesh_grad_trial_ref, @@ -325,6 +329,23 @@ namespace proteus cfl = nrm_v/h; } + inline + void calculateNonlinearCFL(const double& elementDiameter, + const double df[nSpace], + const double norm_factor_lagged, + const double epsFactHeaviside, + const double lambdaFact, + double& cfl) + { + double h,nrm2_v; + h = elementDiameter; + nrm2_v=0.0; + for(int I=0;I eps) + d=0.0; + else if (u < -eps) + d=0.0; + else + d = 0.5*(1.0 + cos(M_PI*u/eps)); + return d; + } + inline double smoothedSign(double eps, double u) { double H; @@ -497,7 +530,9 @@ namespace proteus double* min_distance, double* max_distance, double* mean_distance, + double* volume_domain, double norm_factor_lagged, + double VelMax, // normal reconstruction double* projected_qx_tn, double* projected_qy_tn, @@ -513,6 +548,7 @@ namespace proteus min_distance[0] = 1E10; max_distance[0] = -1E10; mean_distance[0] = 0.; + volume_domain[0] = 0.; for(int eN=0;eN 0 and firstStep==True): # Overwritten if spin up step is taken self.model.u0_dof[:] = self.model.u[0].dof + # + if (firstStep==True): + r=1 + else: + r=self.model.timeIntegration.dt/self.dt_old + # + self.dt_old=self.model.timeIntegration.dt + self.q_v_tStar[:] = (1+r)*self.q_v - r*self.q_v_old + # SAVE OLD VELOCITY # self.q_v_old[:] = self.q_v # COMPUTE NEW VELOCITY (if given by user) # @@ -153,9 +163,12 @@ def preStep(self,t,firstStep=False): self.model.updateVelocityFieldAsFunction() if (firstStep==True): self.q_v_old[:] = self.q_v + self.q_v_tStar[:] = self.q_v # GET NORMALS RECONSTRUCTION; i.e, compute qx_tn, qy_tn. This is done within CLSVOF solver # SAVE OLD SOLUTION (and VELOCITY) # self.model.u_dof_old[:] = self.model.u[0].dof + self.VelMax = max(self.q_v.max(),1E-6) + copyInstructions = {} return copyInstructions @@ -614,7 +627,9 @@ def __init__(self, self.min_distance = 0. self.max_distance = 0. self.mean_distance = 0. + self.volume_domain = 0. self.norm_factor_lagged = 1.0 + self.VelMax = 1.0 self.weighted_lumped_mass_matrix = numpy.zeros(self.u[0].dof.shape,'d') # rhs for normal reconstruction self.rhs_qx = numpy.zeros(self.u[0].dof.shape,'d') @@ -1042,6 +1057,7 @@ def getResidual(self,u,r): min_distance = numpy.zeros(1) max_distance = numpy.zeros(1) mean_distance = numpy.zeros(1) + volume_domain = numpy.zeros(1) self.clsvof.calculateResidual(#element self.timeIntegration.dt, @@ -1119,7 +1135,9 @@ def getResidual(self,u,r): min_distance, max_distance, mean_distance, + volume_domain, self.norm_factor_lagged, + self.VelMax, # normal reconstruction self.projected_qx_tn, self.projected_qy_tn, @@ -1137,7 +1155,9 @@ def getResidual(self,u,r): self.min_distance = -globalMax(-min_distance[0]) self.max_distance = globalMax(max_distance[0]) self.mean_distance = globalSum(mean_distance[0]) - + self.volume_domain = globalSum(volume_domain[0]) + self.mean_distance /= self.volume_domain + if self.forceStrongConditions:# for dofN,g in self.dirichletConditionsForceDOF.DOFBoundaryConditionsDict.iteritems(): r[dofN] = 0 @@ -1216,7 +1236,8 @@ def getJacobian(self,jacobian): self.coefficients.epsFactHeaviside, self.coefficients.epsFactDirac, self.coefficients.lambdaFact, - self.norm_factor_lagged) + self.norm_factor_lagged, + self.VelMax) #Load the Dirichlet conditions directly into residual if self.forceStrongConditions: diff --git a/proteus/mprans/cCLSVOF.pyx b/proteus/mprans/cCLSVOF.pyx index 154d27d406..384234cac7 100644 --- a/proteus/mprans/cCLSVOF.pyx +++ b/proteus/mprans/cCLSVOF.pyx @@ -73,7 +73,9 @@ cdef extern from "CLSVOF.h" namespace "proteus": double* min_distance, double* max_distance, double* mean_distance, + double* volume_domain, double norm_factor_lagged, + double VelMax, double* projected_qx_tn, double* projected_qy_tn, double* projected_qz_tn, @@ -135,7 +137,8 @@ cdef extern from "CLSVOF.h" namespace "proteus": double epsFactHeaviside, double epsFactDirac, double lambdaFact, - double norm_factor_lagged) + double norm_factor_lagged, + double VelMax) void calculateMetricsAtEOS(double* mesh_trial_ref, double* mesh_grad_trial_ref, double* mesh_dof, @@ -349,7 +352,9 @@ cdef class cCLSVOF_base: numpy.ndarray min_distance, numpy.ndarray max_distance, numpy.ndarray mean_distance, + numpy.ndarray volume_domain, double norm_factor_lagged, + double VelMax, numpy.ndarray projected_qx_tn, numpy.ndarray projected_qy_tn, numpy.ndarray projected_qz_tn, @@ -426,7 +431,9 @@ cdef class cCLSVOF_base: min_distance.data, max_distance.data, mean_distance.data, + volume_domain.data, norm_factor_lagged, + VelMax, projected_qx_tn.data, projected_qy_tn.data, projected_qz_tn.data, @@ -490,7 +497,8 @@ cdef class cCLSVOF_base: double epsFactHeaviside, double epsFactDirac, double lambdaFact, - norm_factor_lagged): + norm_factor_lagged, + VelMax): cdef numpy.ndarray rowptr,colind,globalJacobian_a (rowptr,colind,globalJacobian_a) = globalJacobian.getCSRrepresentation() self.thisptr.calculateJacobian(dt, @@ -546,7 +554,8 @@ cdef class cCLSVOF_base: epsFactHeaviside, epsFactDirac, lambdaFact, - norm_factor_lagged) + norm_factor_lagged, + VelMax) def calculateMetricsAtEOS(self, numpy.ndarray mesh_trial_ref, numpy.ndarray mesh_grad_trial_ref, diff --git a/proteus/tests/CLSVOF/pure_level_set/clsvof.py b/proteus/tests/CLSVOF/pure_level_set/clsvof.py index ceb9c138bc..c9b0367118 100644 --- a/proteus/tests/CLSVOF/pure_level_set/clsvof.py +++ b/proteus/tests/CLSVOF/pure_level_set/clsvof.py @@ -107,3 +107,4 @@ def attachModels(self,modelList): self.q_v = np.zeros(self.model.q[('grad(u)',0)].shape,'d') self.ebqe_v = np.zeros(self.model.ebqe[('grad(u)',0)].shape,'d') self.q_v_old = np.zeros(self.model.q[('grad(u)',0)].shape,'d') + self.q_v_tStar = np.zeros(self.model.q[('grad(u)',0)].shape,'d')