Skip to content

Commit

Permalink
CLSVOF lambda fix, scales properly with domain size (#790)
Browse files Browse the repository at this point in the history 
* 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
  • Loading branch information
@jhcollins @cekees
jhcollins authored and cekees committed Aug 28, 2018
1 parent 6fdcd24 commit 72b74d3
Show file tree
Hide file tree
Showing 4 changed files with 106 additions and 16 deletions.
81 changes: 70 additions & 11 deletions proteus/mprans/CLSVOF.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
{
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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<nSpace;I++)
nrm2_v+=df[I]*df[I];
cfl = nrm2_v*norm_factor_lagged/(epsFactHeaviside*lambdaFact*h*h*h);
cfl = std::sqrt(cfl);
}

inline
void exteriorNumericalAdvectiveFlux(const int& isDOFBoundary_u,
const int& isFluxBoundary_u,
Expand Down Expand Up @@ -395,6 +416,18 @@ namespace proteus
return d;
}

inline double smoothedNormalizedDirac(double eps, double u)
{
double d;
if (u > 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;
Expand Down Expand Up @@ -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,
Expand All @@ -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<nElements_global;eN++)
{
Expand Down Expand Up @@ -601,8 +637,15 @@ namespace proteus
mesh_velocity[1] = yt;
mesh_velocity[2] = zt;

double lambda = lambdaFact*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN])/degree_polynomial/norm_factor_lagged;
double epsHeaviside = epsFactHeaviside*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN])/degree_polynomial;
double hK=(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN])/degree_polynomial;
double epsHeaviside = epsFactHeaviside*hK;
double lambda = lambdaFact*hK/norm_factor_lagged;
if (LAMBDA_SCALING==1)
{
double delta = 1.0;//fmax(smoothedNormalizedDirac(2*epsHeaviside,un),1E-6);
lambda = lambdaFact*VelMax*delta;
}

double Sn = smoothedSign(epsHeaviside,un);
double Snp1 = smoothedSign(epsHeaviside,u);
double Hnp1 = smoothedHeaviside(epsHeaviside,u);
Expand All @@ -622,6 +665,7 @@ namespace proteus
relative_velocity[I] = (velocity[eN_k_nSpace+I]-MOVING_DOMAIN*mesh_velocity[I]);
relative_velocity_old[I] = (velocity_old[eN_k_nSpace+I]-MOVING_DOMAIN*mesh_velocity[I]);
fnp1[I] = relative_velocity[I]*Snp1; //implicit advection via BDF
//fnp1[I] = relative_velocity_old[I]*Sn; // explicit advection
fnHalf[I] = 0.5*(relative_velocity[I]*Snp1+relative_velocity_old[I]*Sn); //implicit advection via CN
grad_unHalf[I] = 0.5*(grad_u[I]+grad_un[I]);
}
Expand All @@ -630,6 +674,12 @@ namespace proteus
// CALCULATE CELL BASED CFL //
//////////////////////////////
calculateCFL(elementDiameter[eN]/degree_polynomial,relative_velocity,cfl[eN_k]);
//calculateNonlinearCFL(elementDiameter[eN]/degree_polynomial,
// relative_velocity,
// norm_factor_lagged,
// epsFactHeaviside,
// lambdaFact,
// cfl[eN_k]);

/////////////////////
// TIME DERIVATIVE //
Expand All @@ -639,9 +689,10 @@ namespace proteus
// CALCULATE min, max and mean distance
if (eN<nElements_owned) // locally owned?
{
min_distance[0] = fmin(min_distance[0],u);
max_distance[0] = fmax(max_distance[0],u);
mean_distance[0] += u*dV;
min_distance[0] = fmin(min_distance[0],fabs(u));
max_distance[0] = fmax(max_distance[0],fabs(u));
mean_distance[0] += fabs(u)*dV;
volume_domain[0] += dV;
}

//double norm_grad_un = 0;
Expand Down Expand Up @@ -939,7 +990,8 @@ namespace proteus
double epsFactDirac,
double lambdaFact,
// normalization factor
double norm_factor_lagged)
double norm_factor_lagged,
double VelMax)
{
double timeCoeff=1.0;
if (timeOrder==2)
Expand Down Expand Up @@ -1008,10 +1060,17 @@ namespace proteus
mesh_velocity[1] = yt;
mesh_velocity[2] = zt;

double lambda = lambdaFact*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN])/degree_polynomial/norm_factor_lagged;
double epsDirac = epsFactDirac*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN])/degree_polynomial;
double epsHeaviside = epsFactHeaviside*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN])/degree_polynomial;
double dSnp1 = smoothedDerivativeSign(epsDirac,u); //derivative of smoothed sign
double hK=(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN])/degree_polynomial;
double epsHeaviside = epsFactHeaviside*hK;
double lambda = lambdaFact*hK/norm_factor_lagged;
if (LAMBDA_SCALING==1)
{
double delta = 1.0;//fmax(smoothedNormalizedDirac(2*epsHeaviside,un),1E-6);
lambda = lambdaFact*VelMax*delta;
}

double epsDirac = epsFactDirac*hK;
double dSnp1 = smoothedDerivativeSign(epsDirac,u); //derivative of smoothed sign

for (int I=0;I<nSpace;I++)
{
Expand Down
25 changes: 23 additions & 2 deletions proteus/mprans/CLSVOF.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -92,6 +92,7 @@ def attachModels(self,modelList):
if modelList[self.flowModelIndex].ebq.has_key(('f',0)):
self.ebq_v = modelList[self.flowModelIndex].ebq[('f',0)]
self.q_v_old = numpy.copy(self.q_v)
self.q_v_tStar = numpy.copy(self.q_v)
#VRANS
self.flowCoefficients = modelList[self.flowModelIndex].coefficients
if hasattr(self.flowCoefficients,'q_porosity'):
Expand Down Expand Up @@ -146,16 +147,28 @@ def preStep(self,t,firstStep=False):
if (self.computeMetrics > 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) #
if self.model.hasVelocityFieldAsFunction:
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

Expand Down Expand Up @@ -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')
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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,
Expand All @@ -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
Expand Down Expand Up @@ -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:
Expand Down
15 changes: 12 additions & 3 deletions proteus/mprans/cCLSVOF.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -426,7 +431,9 @@ cdef class cCLSVOF_base:
<double*> min_distance.data,
<double*> max_distance.data,
<double*> mean_distance.data,
<double*> volume_domain.data,
norm_factor_lagged,
VelMax,
<double*> projected_qx_tn.data,
<double*> projected_qy_tn.data,
<double*> projected_qz_tn.data,
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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,
Expand Down
1 change: 1 addition & 0 deletions proteus/tests/CLSVOF/pure_level_set/clsvof.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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')

0 comments on commit 72b74d3

Please sign in to comment.