Skip to content

GitLab

Commit 69838ebe authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

BaseProblems updated and some Warnings removed

parent c00e1666
Hide whitespace changes
Inline Side-by-side
Showing with 132 additions and 119 deletions
extensions/base_problems/CahnHilliardNavierStokes_TwoPhase_RB.cc
View file @ 69838ebe
......@@ -38,7 +38,7 @@ CahnHilliardNavierStokes_TwoPhase_RB::CahnHilliardNavierStokes_TwoPhase_RB(const
Parameters::get(name + "->epsilon", eps); // interface width
// type of double well: 0= [0,1], 1= [-1,1]
Parameters::get(name + "->double-well type", doubleWell);
// Parameters::get(name + "->double-well type", doubleWell);
// parameters for navier-stokes
Initfile::get(name + "->viscosity1", viscosity1);
......@@ -52,12 +52,14 @@ CahnHilliardNavierStokes_TwoPhase_RB::CahnHilliardNavierStokes_TwoPhase_RB(const
// Parameters for CH-Coupling
// 0... u*grad(c), 1... -div(u*c)
Initfile::get(name + "->use conservation form", useConservationForm);
// Initfile::get(name + "->use conservation form", useConservationForm);
// Parameters for NS-Coupling
// cahn-hiliard-force: sigma*mu*grad(c)
Initfile::get(name + "->sigma", sigma); // surface tension
surfaceTension = sigma*3.0/(2.0*sqrt(2.0)) / eps;
Initfile::get(name + "->surfaceTension", surfaceTension); // surface tension
surfaceTensionFactor = 2.0*sqrt(2.0)*eps/(3.0*sigma);
force.set(0.0);
......@@ -69,16 +71,22 @@ CahnHilliardNavierStokes_TwoPhase_RB::CahnHilliardNavierStokes_TwoPhase_RB(const
minusEpsInv = -epsInv;
epsSqr = sqr(eps);
minusEpsSqr = -epsSqr;
MSG("Simulation parameters:\n");
MSG("======================\n");
MSG("gamma=%e\n",gamma);
MSG("epsilon=%e\n",eps);
MSG("viscosity1=%e\n",viscosity1);
MSG("viscosity2=%e\n",viscosity2);
MSG("density1=%e\n",density1);
MSG("density2=%e\n",density2);
MSG("sigma=%e\n",sigma);
}
CahnHilliardNavierStokes_TwoPhase_RB::~CahnHilliardNavierStokes_TwoPhase_RB()
{ FUNCNAME("CahnHilliardNavierStokes_TwoPhase_RB::~CahnHilliardNavierStokes_TwoPhase_RB()");
if (reinit != NULL) {
delete reinit;
reinit = NULL;
}
if (velocity != NULL) {
delete velocity;
velocity = NULL;
......@@ -92,8 +100,6 @@ void CahnHilliardNavierStokes_TwoPhase_RB::initData()
{ FUNCNAME("CahnHilliardNavierStokes_TwoPhase_RB::initData()");
dim = getMesh()->getDim();
// create instance redistancing class
reinit = new HL_SignedDistTraverse("reinit", dim);
if (velocity == NULL)
velocity = new DOFVector<WorldVector<double> >(getFeSpace(2), "velocity");
......@@ -129,52 +135,102 @@ void CahnHilliardNavierStokes_TwoPhase_RB::fillOperators()
WorldVector<DOFVector<double>*> stage_velocity;
WorldVector<DOFVector<double>*> un_velocity;
for (size_t i = 0; i < dow; ++i) {
stage_velocity[i] = prob->getStageSolution(2+i);
un_velocity[i] = prob->getUnVec(2+i);
stage_velocity[i] = prob->getStageSolution(2+i); // intermediateSolution
un_velocity[i] = prob->getUnVec(2+i); // oldSolution
}
// mu + eps^2*laplace(c) - M'(c)
// ----------------------------------------------------------------------
/// < mu , psi >
Operator *opChMu = new Operator(prob->getFeSpace(0),prob->getFeSpace(0));
opChMu->addZeroOrderTerm(new Simple_ZOT(surfaceTensionFactor));
opChMu->setUhOld(prob->getStageSolution(1));
prob->addVectorOperator(*opChMu, 0);
prob->addMatrixOperator(*opChMu, 0, 1, &minus1, &minus1);
/// < -eps^2*grad(c) , grad(psi) >
Operator *opChLC = new Operator(prob->getFeSpace(1),prob->getFeSpace(0));
opChLC->addSecondOrderTerm(new Simple_SOT(minusEpsSqr));
opChLC->setUhOld(prob->getStageSolution(0));
prob->addVectorOperator(*opChLC, 0);
prob->addMatrixOperator(*opChLC, 0, 0, &minus1, &minus1);
/// < -M'(c) , psi >
// if (doubleWell == 0) {
// // jacobian operators
// Operator *opChMImpl = new Operator(prob->getFeSpace(1),prob->getFeSpace(0));
// opChMImpl->addZeroOrderTerm(new VecAtQP_ZOT(
// prob->getUnVec(0),
// new DoubleWell0Diff(-1.0))); // < (3c^2-3c+0.5)*c' , psi >
// prob->addMatrixOperator(*opChMImpl, 1, 0, &minus1, &minus1);
// // rhs operators
// Operator *opChMExpl = new Operator(prob->getFeSpace(1),prob->getFeSpace(0));
// opChMExpl->addZeroOrderTerm(new VecAtQP_ZOT(
// prob->getStageSolution(0),
// new DoubleWell0(-1.0))); // < (c^3-3/2*c^2+0.5) , psi >
// prob->addVectorOperator(*opChMExpl, 1);
// }
// else if (doubleWell == 1) {
// jacobian operators
Operator *opChMImpl = new Operator(prob->getFeSpace(1),prob->getFeSpace(0));
opChMImpl->addZeroOrderTerm(new VecAtQP_ZOT(
prob->getUnVec(0),
new DoubleWell1Diff(-1.0))); // < (3c^2-1)*c' , psi >
prob->addMatrixOperator(*opChMImpl, 0, 0, &minus1, &minus1);
// rhs operators
Operator *opChMExpl = new Operator(prob->getFeSpace(1),prob->getFeSpace(0));
opChMExpl->addZeroOrderTerm(new VecAtQP_ZOT(
prob->getStageSolution(0),
new DoubleWell1(-1.0))); // < (c^3-c) , psi >
prob->addVectorOperator(*opChMExpl, 0);
// }
// dt(c) = laplace(mu) - u*grad(c)
// -----------------------------------
/// < dt(c) , psi >
prob->addTimeOperator(0, 0);
/// d_t(c) + u*grad(c)
prob->addTimeOperator(1,0);
// div(M(c)grad(mu)), with M(c)=gamma/4*(c^2-1)^2
Operator *opChLM = new Operator(prob->getFeSpace(0),prob->getFeSpace(1));
if (useMobility) {
if (doubleWell == 0) {
// jacobian operators
opChLM->addTerm(new VecAtQP_SOT(
prob->getUnVec(0),
new MobilityCH0(gamma)));
prob->addMatrixOperator(*opChLM, 0, 1);
Operator *opChLM2 = new Operator(prob->getFeSpace(0),prob->getFeSpace(0));
opChLM2->addTerm(new VecGrad_FOT(
prob->getUnVec(0),
prob->getUnVec(1),
new MobilityCH0Diff(gamma)), GRD_PSI);
prob->addMatrixOperator(*opChLM2, 0, 0);
// rhs operators
Operator *opChLMe = new Operator(prob->getFeSpace(0),prob->getFeSpace(1));
opChLMe->addTerm(new VecAtQP_SOT(
prob->getStageSolution(0),
new MobilityCH0(gamma)));
opChLMe->setUhOld(prob->getStageSolution(1));
prob->addVectorOperator(*opChLMe, 0, &minus1, &minus1);
} else if (doubleWell == 1) {
// if (doubleWell == 0) {
// // jacobian operators
// opChLM->addTerm(new VecAtQP_SOT(
// prob->getUnVec(0),
// new MobilityCH0(gamma)));
// prob->addMatrixOperator(*opChLM, 0, 1);
//
// Operator *opChLM2 = new Operator(prob->getFeSpace(0),prob->getFeSpace(0));
// opChLM2->addTerm(new VecGrad_FOT(
// prob->getUnVec(0),
// prob->getUnVec(1),
// new MobilityCH0Diff(gamma)), GRD_PSI);
// prob->addMatrixOperator(*opChLM2, 0, 0);
//
// // rhs operators
// Operator *opChLMe = new Operator(prob->getFeSpace(0),prob->getFeSpace(1));
// opChLMe->addTerm(new VecAtQP_SOT(
// prob->getStageSolution(0),
// new MobilityCH0(gamma)));
// opChLMe->setUhOld(prob->getStageSolution(1));
// prob->addVectorOperator(*opChLMe, 0, &minus1, &minus1);
// } else if (doubleWell == 1) {
// jacobian operators
opChLM->addTerm(new VecAtQP_SOT(
prob->getUnVec(0),
new MobilityCH1(gamma)));
prob->addMatrixOperator(*opChLM, 0, 1);
prob->addMatrixOperator(*opChLM, 1, 1);
Operator *opChLM2 = new Operator(prob->getFeSpace(0),prob->getFeSpace(0));
opChLM2->addTerm(new VecGrad_FOT(
prob->getUnVec(0),
prob->getUnVec(1),
new MobilityCH1Diff(gamma)), GRD_PSI);
prob->addMatrixOperator(*opChLM2, 0, 0);
prob->addMatrixOperator(*opChLM2, 1, 0);
// rhs operators
Operator *opChLMe = new Operator(prob->getFeSpace(0),prob->getFeSpace(1));
......@@ -182,13 +238,13 @@ void CahnHilliardNavierStokes_TwoPhase_RB::fillOperators()
prob->getStageSolution(0),
new MobilityCH1(gamma)));
opChLMe->setUhOld(prob->getStageSolution(1));
prob->addVectorOperator(*opChLMe, 0, &minus1, &minus1);
}
prob->addVectorOperator(*opChLMe, 1, &minus1, &minus1);
// }
} else {
opChLM->addTerm(new Simple_SOT(gamma));
opChLM->setUhOld(prob->getStageSolution(1));
prob->addVectorOperator(*opChLM, 0, &minus1, &minus1);
prob->addMatrixOperator(*opChLM, 0, 1); /// < phi*grad(mu) , grad(psi) >
prob->addVectorOperator(*opChLM, 1, &minus1, &minus1);
prob->addMatrixOperator(*opChLM, 1, 1); /// < phi*grad(mu) , grad(psi) >
}
/// < u_s * grad(c_s), psi >
......@@ -198,7 +254,7 @@ void CahnHilliardNavierStokes_TwoPhase_RB::fillOperators()
else
uGradC->addTerm(new WorldVec_FOT(stage_velocity), GRD_PHI);
uGradC->setUhOld(prob->getStageSolution(0));
prob->addVectorOperator(*uGradC, 0, &minus1, &minus1);
prob->addVectorOperator(*uGradC, 1, &minus1, &minus1);
/// < u_old * grad(c), psi >
Operator *JuGradC1 = new Operator(prob->getFeSpace(0), prob->getFeSpace(0));
......@@ -206,7 +262,7 @@ void CahnHilliardNavierStokes_TwoPhase_RB::fillOperators()
JuGradC1->addTerm(new WorldVec_FOT(un_velocity, -1.0), GRD_PSI);
else
JuGradC1->addTerm(new WorldVec_FOT(un_velocity), GRD_PHI);
prob->addMatrixOperator(*JuGradC1, 0, 0);
prob->addMatrixOperator(*JuGradC1, 1, 0);
/// < u * grad(c_old), psi >
for (size_t i = 0; i < dow; i++) {
......@@ -216,56 +272,9 @@ void CahnHilliardNavierStokes_TwoPhase_RB::fillOperators()
else {
JuGradC2->addTerm(new PartialDerivative_ZOT(prob->getUnVec(0), i));
}
prob->addMatrixOperator(*JuGradC2, 0, 2+i);
}
// mu + eps^2*laplace(c) - M'(c)
// ----------------------------------------------------------------------
/// < mu , psi >
Operator *opChMu = new Operator(prob->getFeSpace(0),prob->getFeSpace(0));
opChMu->addZeroOrderTerm(new Simple_ZOT);
opChMu->setUhOld(prob->getStageSolution(1));
prob->addVectorOperator(*opChMu, 1);
prob->addMatrixOperator(*opChMu, 1, 1, &minus1, &minus1);
/// < -eps^2*grad(c) , grad(psi) >
Operator *opChLC = new Operator(prob->getFeSpace(1),prob->getFeSpace(0));
opChLC->addSecondOrderTerm(new Simple_SOT(minusEpsSqr));
opChLC->setUhOld(prob->getStageSolution(0));
prob->addVectorOperator(*opChLC, 1);
prob->addMatrixOperator(*opChLC, 1, 0, &minus1, &minus1);
/// < -M'(c) , psi >
if (doubleWell == 0) {
// jacobian operators
Operator *opChMImpl = new Operator(prob->getFeSpace(1),prob->getFeSpace(0));
opChMImpl->addZeroOrderTerm(new VecAtQP_ZOT(
prob->getUnVec(0),
new DoubleWell0Diff(-1.0))); // < (3c^2-3c+0.5)*c' , psi >
prob->addMatrixOperator(*opChMImpl, 1, 0, &minus1, &minus1);
// rhs operators
Operator *opChMExpl = new Operator(prob->getFeSpace(1),prob->getFeSpace(0));
opChMExpl->addZeroOrderTerm(new VecAtQP_ZOT(
prob->getStageSolution(0),
new DoubleWell0(-1.0))); // < (c^3-3/2*c^2+0.5) , psi >
prob->addVectorOperator(*opChMExpl, 1);
}
else if (doubleWell == 1) {
// jacobian operators
Operator *opChMImpl = new Operator(prob->getFeSpace(1),prob->getFeSpace(0));
opChMImpl->addZeroOrderTerm(new VecAtQP_ZOT(
prob->getUnVec(0),
new DoubleWell1Diff(-1.0))); // < (3c^2-1)*c' , psi >
prob->addMatrixOperator(*opChMImpl, 1, 0, &minus1, &minus1);
// rhs operators
Operator *opChMExpl = new Operator(prob->getFeSpace(1),prob->getFeSpace(0));
opChMExpl->addZeroOrderTerm(new VecAtQP_ZOT(
prob->getStageSolution(0),
new DoubleWell1(-1.0))); // < (c^3-c) , psi >
prob->addVectorOperator(*opChMExpl, 1);
prob->addMatrixOperator(*JuGradC2, 1, 2+i);
}
// Navier-Stokes part
// -----------------------------------------------------------------------------------
bool variableDensity = abs(density1 - density2) > DBL_TOL;
......@@ -319,26 +328,25 @@ void CahnHilliardNavierStokes_TwoPhase_RB::fillOperators()
// forces by Cahn-Hilliard terms
// -----------------------------
if (surfaceTension > DBL_TOL) {
/// < mu_old * grad(c) , psi >
Operator *JmuGradC = new Operator(getFeSpace(2+i),getFeSpace(0));
JmuGradC->addTerm(new VecAndPartialDerivative_FOT(
prob->getUnVec(1), i, -1.0), GRD_PHI);
prob->addMatrixOperator(*JmuGradC, 2+i, 0, &surfaceTension, &surfaceTension);
/// < mu * grad(c_old) , psi >
Operator *JmuGradC2 = new Operator(getFeSpace(2+i),getFeSpace(1));
JmuGradC2->addTerm(new PartialDerivative_ZOT(
prob->getUnVec(0), i, -1.0));
prob->addMatrixOperator(*JmuGradC2, 2+i, 1, &surfaceTension, &surfaceTension);
/// < mu_s * grad(c_s) , psi >
Operator *opMuGradC = new Operator(getFeSpace(2+i),getFeSpace(1));
opMuGradC->addTerm(new VecAndPartialDerivative_ZOT(
prob->getStageSolution(1),
prob->getStageSolution(0), i));
prob->addVectorOperator(*opMuGradC, 2+i, &surfaceTension, &surfaceTension);
}
/// < mu_old * grad(c) , psi >
Operator *JmuGradC = new Operator(getFeSpace(2+i),getFeSpace(0));
JmuGradC->addTerm(new VecAndPartialDerivative_FOT(
prob->getUnVec(1), i, -1.0), GRD_PHI);
prob->addMatrixOperator(*JmuGradC, 2+i, 0);
/// < mu * grad(c_old) , psi >
Operator *JmuGradC2 = new Operator(getFeSpace(2+i),getFeSpace(1));
JmuGradC2->addTerm(new PartialDerivative_ZOT(
prob->getUnVec(0), i, -1.0));
prob->addMatrixOperator(*JmuGradC2, 2+i, 1);
/// < mu_s * grad(c_s) , psi >
Operator *opMuGradC = new Operator(getFeSpace(2+i),getFeSpace(1));
opMuGradC->addTerm(new VecAndPartialDerivative_ZOT(
prob->getStageSolution(1),
prob->getStageSolution(0), i));
prob->addVectorOperator(*opMuGradC, 2+i);
// ---------------------------------------
if (variableDensity) {
......@@ -433,7 +441,7 @@ void CahnHilliardNavierStokes_TwoPhase_RB::addStressTerm(int i)
opLaplaceUijExpl->addTerm(new VecAtQP_IJ_SOT(
prob->getStageSolution(0),
new LinearInterpolation(viscosity1, viscosity2),
1-i, 1-j));
j, i));
opLaplaceUijExpl->setUhOld(prob->getStageSolution(2+j));
prob->addVectorOperator(*opLaplaceUijExpl, 2+i, &minus1, &minus1);
......@@ -442,7 +450,7 @@ void CahnHilliardNavierStokes_TwoPhase_RB::addStressTerm(int i)
opLaplaceUijImpl->addTerm(new VecAtQP_IJ_SOT(
prob->getUnVec(0),
new LinearInterpolation(viscosity1, viscosity2),
1-i, 1-j));
j, i));
prob->addMatrixOperator(*opLaplaceUijImpl, 2+i, 2+j);
/// < c * nu'(c_old)*d_j(u_old_i) , d_i(psi) >
......@@ -450,7 +458,7 @@ void CahnHilliardNavierStokes_TwoPhase_RB::addStressTerm(int i)
opLaplaceUijImpl2->addTerm(new VecAndPartialDerivativeIJ_FOT(
prob->getUnVec(0),
prob->getUnVec(2+j),
1-i, 1-j,
j, i,
new LinearInterpolationDiffFactor(viscosity1, viscosity2)), GRD_PSI);
prob->addMatrixOperator(*opLaplaceUijImpl2, 2+i, 0);
}
......
extensions/base_problems/CahnHilliardNavierStokes_TwoPhase_RB.h
View file @ 69838ebe
......@@ -5,7 +5,6 @@
#include "AMDiS.h"
#include "BaseProblem_RB.h"
#include "HL_SignedDistTraverse.h"
#include "chns.h"
using namespace AMDiS;
......@@ -38,11 +37,16 @@ public: // public methods
virtual void fillBoundaryConditions() {}
// virtual void addNonlinearTerm(int i);
virtual void addStressTerm(int i);
virtual void addMaterialDerivativeConst(int i, int j,double factor = 1.0); // rho = const
virtual void addMaterialDerivativeNonConst(int i, int j); // rho = rho(x)
virtual void addLaplaceTermConst(int i); // viscosity = const
virtual void addLaplaceTermNonConst(int i); // viscosity = viscosity(x)
protected: // protected variables
HL_SignedDistTraverse *reinit;
DOFVector<WorldVector<double> >* velocity;
void calcVelocity()
......@@ -56,13 +60,13 @@ protected: // protected variables
}
bool useMobility;
bool useConservationForm;
// bool useConservationForm;
unsigned dim;
int laplaceType;
int nonLinTerm;
int doubleWell;
// int doubleWell;
// navierStokes parameters
double oldTimestep;
......@@ -83,6 +87,7 @@ protected: // protected variables
double sigma; // coupling parameter to calculate the surface tension
double surfaceTension;// := sigma/epsilon
double surfaceTensionFactor;// := epsilon*2*sqrt(2)/(3*sigma)
WorldVector<double> force;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment