From 4baf4abeb8143d82f47cdbf65b16b4e75e595af7 Mon Sep 17 00:00:00 2001
From: Simon Praetorius <simon.praetorius@tu-dresden.de>
Date: Wed, 15 Aug 2012 10:18:44 +0000
Subject: [PATCH] base problems
---
.../base_problems/CahnHilliardNavierStokes.cc | 74 ++-
.../base_problems/CahnHilliardNavierStokes.h | 5 +-
.../CahnHilliardNavierStokes_RB.cc | 369 ++++++++++++++
.../CahnHilliardNavierStokes_RB.h | 272 ++++++++++
.../CahnHilliardNavierStokes_TwoPhase_RB.cc | 477 ++++++++++++++++++
.../CahnHilliardNavierStokes_TwoPhase_RB.h | 326 ++++++++++++
extensions/base_problems/DiffuseDomainFsi.cc | 184 ++++---
extensions/base_problems/DiffuseDomainFsi.h | 9 +-
extensions/base_problems/LinearElasticity.cc | 44 +-
extensions/base_problems/LinearElasticity.h | 2 +
.../base_problems/NavierStokes_TaylorHood.cc | 10 +-
11 files changed, 1666 insertions(+), 106 deletions(-)
create mode 100644 extensions/base_problems/CahnHilliardNavierStokes_RB.cc
create mode 100644 extensions/base_problems/CahnHilliardNavierStokes_RB.h
create mode 100644 extensions/base_problems/CahnHilliardNavierStokes_TwoPhase_RB.cc
create mode 100644 extensions/base_problems/CahnHilliardNavierStokes_TwoPhase_RB.h
diff --git a/extensions/base_problems/CahnHilliardNavierStokes.cc b/extensions/base_problems/CahnHilliardNavierStokes.cc
index 6ef2b8ce..c194ad5a 100644
--- a/extensions/base_problems/CahnHilliardNavierStokes.cc
+++ b/extensions/base_problems/CahnHilliardNavierStokes.cc
@@ -10,7 +10,7 @@ CahnHilliardNavierStokes::CahnHilliardNavierStokes(const std::string &name_) :
super(name_),
useMobility(false),
useConservationForm(false),
- doubleWell(0),
+ doubleWell(1),
gamma(1.0),
eps(0.1),
minusEps(-0.1),
@@ -19,7 +19,7 @@ CahnHilliardNavierStokes::CahnHilliardNavierStokes(const std::string &name_) :
epsSqr(0.01),
minusEpsSqr(-0.01),
laplaceType(0),
- nonLinTerm(2),
+ nonLinTerm(3),
oldTimestep(0.0),
viscosity(1.0),
density(1.0),
@@ -31,8 +31,8 @@ CahnHilliardNavierStokes::CahnHilliardNavierStokes(const std::string &name_) :
fileWriter(NULL)
{
// parameters for CH
- Parameters::get(name + "->use mobility", useMobility); // mobility
- Parameters::get(name + "->gamma", gamma); // mobility
+ Parameters::get(name + "->use mobility", useMobility); // dt(c) = div(M(c)*grad(c)) + ...
+ Parameters::get(name + "->gamma", gamma); // M(c) = gamma * Double-Well
Parameters::get(name + "->epsilon", eps); // interface width
// type of double well: 0= [0,1], 1= [-1,1]
@@ -41,15 +41,23 @@ CahnHilliardNavierStokes::CahnHilliardNavierStokes(const std::string &name_) :
// parameters for navier-stokes
Initfile::get(name + "->viscosity", viscosity);
Initfile::get(name + "->density", density);
- // type of laplace operator: 0... div(nu*grad(u)), 1... div(0.5*nu*(grad(u)+grad(u)^T))
+
+ // type of laplace operator: 0... div(nu*grad(u)),
+ // 1... div(0.5*nu*(grad(u)+grad(u)^T))
Initfile::get(name + "->laplace operator", laplaceType);
- // type of non-linear term: 0... u^old*grad(u_i^old), 1... u'*grad(u_i^old), 2... u^old*grad(u'_i)
+
+ // type of non-linear term: 0... u^old*grad(u_i^old),
+ // 1... u*grad(u_i^old),
+ // 2... u^old*grad(u_i)
+ // 3... u*grad(u_i^old) + u^old*grad(u_i) - u^old*grad(u_i^old)
Initfile::get(name + "->non-linear term", nonLinTerm);
// Parameters for CH-Coupling
+ // 0... u*grad(c), 1... -div(u*c)
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;
@@ -97,6 +105,22 @@ void CahnHilliardNavierStokes::initData()
}
+void CahnHilliardNavierStokes::transferInitialSolution(AdaptInfo *adaptInfo)
+{ FUNCNAME("CahnHilliardNavierStokes::transferInitialSolution()");
+
+ calcVelocity();
+
+ for (int i = 0; i < 2+dow; i++)
+ prob->setExactSolution(prob->getSolution()->getDOFVector(i), i);
+
+ fileWriter->writeFiles(adaptInfo, false);
+ writeFiles(adaptInfo, false);
+
+ // initial parameters for detecting mesh changes
+ oldMeshChangeIdx= getMesh()->getChangeIndex();
+}
+
+
void CahnHilliardNavierStokes::fillOperators()
{ FUNCNAME("CahnHilliardNavierStokes::fillOperators()");
MSG("CahnHilliardNavierStokes::fillOperators()\n");
@@ -143,9 +167,9 @@ void CahnHilliardNavierStokes::fillOperators()
uGradC->addTerm(new WorldVec_FOT(vel, -1.0), GRD_PSI);
else
uGradC->addTerm(new WorldVec_FOT(vel), GRD_PHI);
+ prob->addMatrixOperator(*uGradC, 0, 0);
uGradC->setUhOld(prob->getSolution()->getDOFVector(0));
prob->addVectorOperator(*uGradC, 0);
- prob->addMatrixOperator(*uGradC, 0, 0);
/// < u * grad(c_old) , psi >
for (size_t i = 0; i < dow; i++) {
@@ -208,27 +232,28 @@ void CahnHilliardNavierStokes::fillOperators()
prob->addMatrixOperator(*opTime, 2+i, 2+i, getInvTau(), getInvTau());
/// < (1/tau)*u_i^old , psi >
Operator *opTimeOld = new Operator(getFeSpace(2+i), getFeSpace(2+i));
- opTime->addTerm(new Phase_ZOT(prob->getSolution()->getDOFVector(2+i), density));
+ opTimeOld->addTerm(new Phase_ZOT(prob->getSolution()->getDOFVector(2+i), density));
prob->addVectorOperator(*opTimeOld, 2+i, getInvTau(), getInvTau());
/// < u^old*grad(u_i^old) , psi >
- Operator *opUGradU0 = new Operator(getFeSpace(2+i), getFeSpace(2+i));
- opUGradU0->addTerm(new WorldVector_FOT(velocity, -density), GRD_PHI);
- opUGradU0->setUhOld(prob->getSolution()->getDOFVector(2+i));
- if (nonLinTerm == 0) {
+ if (nonLinTerm == 0 || nonLinTerm == 3) {
+ Operator *opUGradU0 = new Operator(getFeSpace(2+i), getFeSpace(2+i));
+ opUGradU0->addTerm(new WorldVec_FOT(vel, (nonLinTerm == 0 ? -1.0 : 1.0)*density),
+ GRD_PHI);
+ opUGradU0->setUhOld(prob->getSolution()->getDOFVector(2+i));
prob->addVectorOperator(*opUGradU0, 2+i);
}
-
- if (nonLinTerm == 1) {
- /// < u'*grad(u_i^old) , psi >
+ if (nonLinTerm == 1 || nonLinTerm == 3) {
+ /// < u*grad(u_i^old) , psi >
for (size_t j = 0; j < dow; ++j) {
Operator *opUGradU1 = new Operator(getFeSpace(2+i),getFeSpace(2+i));
opUGradU1->addTerm(new PartialDerivative_ZOT(
prob->getSolution()->getDOFVector(2+i), j, density));
prob->addMatrixOperator(*opUGradU1, 2+i, 2+j);
}
- } else if (nonLinTerm == 2) {
- /// < u^old*grad(u'_i) , psi >
+ }
+ if (nonLinTerm == 2 || nonLinTerm == 3) {
+ /// < u^old*grad(u_i) , psi >
for(size_t j = 0; j < dow; ++j) {
Operator *opUGradU2 = new Operator(getFeSpace(2+i),getFeSpace(2+i));
opUGradU2->addTerm(new VecAndPartialDerivative_FOT(
@@ -258,16 +283,21 @@ void CahnHilliardNavierStokes::fillOperators()
prob->addVectorOperator(*opForce, 2+i);
}
- /// forces by Cahn-Hilliard terms
+ // forces by Cahn-Hilliard terms
+ // -----------------------------
+ /// < mu_old * grad(c) , psi >
Operator *opMuGradC = new Operator(getFeSpace(2+i),getFeSpace(0));
opMuGradC->addTerm(new VecAndPartialDerivative_FOT(
prob->getSolution()->getDOFVector(1), i), GRD_PHI);
prob->addMatrixOperator(*opMuGradC, 2+i, 0, &surfaceTension, &surfaceTension);
+ /// < mu * grad(c_old) , psi >
Operator *opMuGradC2 = new Operator(getFeSpace(2+i),getFeSpace(1));
opMuGradC2->addTerm(new PartialDerivative_ZOT(
prob->getSolution()->getDOFVector(0), i));
prob->addMatrixOperator(*opMuGradC2, 2+i, 1, &surfaceTension, &surfaceTension);
+
+ /// < mu_old * grad(c_old) , psi >
opMuGradC2->setUhOld(prob->getSolution()->getDOFVector(1));
prob->addVectorOperator(*opMuGradC2, 2+i, &surfaceTension, &surfaceTension);
@@ -300,3 +330,11 @@ void CahnHilliardNavierStokes::addLaplaceTerm(int i)
opLaplaceUi->addTerm(new Simple_SOT(viscosity));
prob->addMatrixOperator(*opLaplaceUi, 2+i, 2+i);
}
+
+void CahnHilliardNavierStokes::closeTimestep(AdaptInfo *adaptInfo)
+{ FUNCNAME("CahnHilliardNavierStokes::closeTimestep()");
+
+ calcVelocity();
+ fileWriter->writeFiles(adaptInfo, false);
+ writeFiles(adaptInfo, false);
+}
diff --git a/extensions/base_problems/CahnHilliardNavierStokes.h b/extensions/base_problems/CahnHilliardNavierStokes.h
index 56686fdd..54d746e1 100644
--- a/extensions/base_problems/CahnHilliardNavierStokes.h
+++ b/extensions/base_problems/CahnHilliardNavierStokes.h
@@ -23,6 +23,9 @@ public: // public methods
virtual void initData();
+ virtual void transferInitialSolution(AdaptInfo *adaptInfo);
+ virtual void closeTimestep(AdaptInfo *adaptInfo);
+
double getEpsilon() { return eps; }
int getDoubleWellType() { return doubleWell; }
@@ -31,8 +34,6 @@ public: // public methods
return velocity;
}
-protected: // protected methods
-
virtual void fillOperators();
virtual void fillBoundaryConditions() {}
virtual void addLaplaceTerm(int i);
diff --git a/extensions/base_problems/CahnHilliardNavierStokes_RB.cc b/extensions/base_problems/CahnHilliardNavierStokes_RB.cc
new file mode 100644
index 00000000..70a1b7be
--- /dev/null
+++ b/extensions/base_problems/CahnHilliardNavierStokes_RB.cc
@@ -0,0 +1,369 @@
+#include "CahnHilliardNavierStokes_RB.h"
+#include "Views.h"
+#include "SignedDistFunctors.h"
+#include "PhaseFieldConvert.h"
+#include "POperators.h"
+
+using namespace AMDiS;
+
+CahnHilliardNavierStokes_RB::CahnHilliardNavierStokes_RB(const std::string &name_) :
+ super(name_),
+ useMobility(false),
+ useConservationForm(false),
+ doubleWell(1),
+ gamma(1.0),
+ eps(0.1),
+ minusEps(-0.1),
+ epsInv(10.0),
+ minusEpsInv(-10.0),
+ epsSqr(0.01),
+ minusEpsSqr(-0.01),
+ laplaceType(0),
+ nonLinTerm(3),
+ oldTimestep(0.0),
+ viscosity(1.0),
+ density(1.0),
+ c(0.0),
+ sigma(1.0),
+ surfaceTension(1.0),
+ reinit(NULL),
+ velocity(NULL),
+ fileWriter(NULL)
+{
+ // parameters for CH
+ Parameters::get(name + "->use mobility", useMobility); // dt(c) = div(M(c)*grad(c)) + ...
+ Parameters::get(name + "->gamma", gamma); // M(c) = gamma * Double-Well
+ 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 for navier-stokes
+ Initfile::get(name + "->viscosity", viscosity);
+ Initfile::get(name + "->density", density);
+
+ // type of laplace operator: 0... div(nu*grad(u)),
+ // 1... div(0.5*nu*(grad(u)+grad(u)^T))
+ Initfile::get(name + "->laplace operator", laplaceType);
+
+ // type of non-linear term: 0... u^old*grad(u_i^old),
+ // 1... u*grad(u_i^old),
+ // 2... u^old*grad(u_i)
+ // 3... u*grad(u_i^old) + u^old*grad(u_i) - u^old*grad(u_i^old)
+ Initfile::get(name + "->non-linear term", nonLinTerm);
+
+ // Parameters for CH-Coupling
+ // 0... u*grad(c), 1... -div(u*c)
+ 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;
+
+
+ force.set(0.0);
+ Initfile::get(name + "->force", force);
+
+ // transformation of the parameters
+ minusEps = -eps;
+ epsInv = 1.0/eps;
+ minusEpsInv = -epsInv;
+ epsSqr = sqr(eps);
+ minusEpsSqr = -epsSqr;
+}
+
+
+CahnHilliardNavierStokes_RB::~CahnHilliardNavierStokes_RB()
+{ FUNCNAME("CahnHilliardNavierStokes_RB::~CahnHilliardNavierStokes_RB()");
+
+ if (reinit != NULL) {
+ delete reinit;
+ reinit = NULL;
+ }
+ if (velocity != NULL) {
+ delete velocity;
+ velocity = NULL;
+ }
+ delete fileWriter;
+ fileWriter = NULL;
+}
+
+
+void CahnHilliardNavierStokes_RB::initData()
+{ FUNCNAME("CahnHilliardNavierStokes_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");
+ fileWriter = new FileVectorWriter(name + "->velocity->output", getFeSpace(2)->getMesh(), velocity);
+
+ super::initData();
+}
+
+
+void CahnHilliardNavierStokes_RB::transferInitialSolution(AdaptInfo *adaptInfo)
+{ FUNCNAME("CahnHilliardNavierStokes_RB::transferInitialSolution()");
+
+ calcVelocity();
+
+ for (int i = 0; i < 2+dow; i++)
+ prob->setExactSolution(prob->getSolution()->getDOFVector(i), i);
+
+ fileWriter->writeFiles(adaptInfo, false);
+ writeFiles(adaptInfo, false);
+
+ // initial parameters for detecting mesh changes
+ oldMeshChangeIdx= getMesh()->getChangeIndex();
+}
+
+
+void CahnHilliardNavierStokes_RB::fillOperators()
+{ FUNCNAME("CahnHilliardNavierStokes_RB::fillOperators()");
+ MSG("CahnHilliardNavierStokes_RB::fillOperators()\n");
+
+ // variable order:
+ // (c, mu, u0, u1 [, u2], p)
+
+ 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);
+ }
+
+ // dt(c) = laplace(mu) - u*grad(c)
+ // -----------------------------------
+ /// < dt(c) , psi >
+ prob->addTimeOperator(0, 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)));
+ opChLM->addTerm(new VecGrad_FOT(
+ prob->getUnVec(0),
+ prob->getUnVec(1),
+ new MobilityCH0Diff(gamma)), GRD_PSI);
+ prob->addMatrixOperator(*opChLM, 0, 1);
+ // rhs operators
+ Operator *opChLM2 = new Operator(prob->getFeSpace(0),prob->getFeSpace(1));
+ opChLM2->addTerm(new VecAtQP_SOT(
+ prob->getStageSolution(0),
+ new MobilityCH0(gamma)));
+ opChLM2->setUhOld(prob->getStageSolution(1));
+ prob->addVectorOperator(*opChLM2, 0, &minus1, &minus1);
+ } else {
+ // jacobian operators
+ opChLM->addTerm(new VecAtQP_SOT(
+ prob->getUnVec(0),
+ new MobilityCH1(gamma)));
+ opChLM->addTerm(new VecGrad_FOT(
+ prob->getUnVec(0),
+ prob->getUnVec(1),
+ new MobilityCH1Diff(gamma)), GRD_PSI);
+ prob->addMatrixOperator(*opChLM, 0, 1);
+ // rhs operators
+ Operator *opChLM2 = new Operator(prob->getFeSpace(0),prob->getFeSpace(1));
+ opChLM2->addTerm(new VecAtQP_SOT(
+ prob->getStageSolution(0),
+ new MobilityCH1(gamma)));
+ opChLM2->setUhOld(prob->getStageSolution(1));
+ prob->addVectorOperator(*opChLM2, 0, &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) >
+ }
+
+ /// < u_s * grad(c_s), psi >
+ Operator *uGradC = new Operator(prob->getFeSpace(0), prob->getFeSpace(0));
+ if (useConservationForm)
+ uGradC->addTerm(new WorldVec_FOT(stage_velocity, -1.0), GRD_PSI);
+ else
+ uGradC->addTerm(new WorldVec_FOT(stage_velocity), GRD_PHI);
+ uGradC->setUhOld(prob->getStageSolution(0));
+ prob->addVectorOperator(*uGradC, 0, &minus1, &minus1);
+
+ /// < u_old * grad(c), psi >
+ Operator *JuGradC1 = new Operator(prob->getFeSpace(0), prob->getFeSpace(0));
+ if (useConservationForm)
+ 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);
+
+ /// < u * grad(c_old), psi >
+ for (size_t i = 0; i < dow; i++) {
+ Operator *JuGradC2 = new Operator(prob->getFeSpace(0), prob->getFeSpace(0));
+ if (useConservationForm)
+ JuGradC2->addTerm(new VecAndPartialDerivative_FOT(prob->getUnVec(0), i, -1.0), GRD_PSI);
+ 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);
+ }
+
+ // Navier-Stokes part
+ // -----------------------------------------------------------------------------------
+ for (size_t i = 0; i < dow; ++i) {
+ /// < d_t(u_i) , psi >
+ prob->addTimeOperator(2+i, 2+i);
+
+ /// < u_s*grad(u_s_i) , psi >
+ Operator *opUGradU = new Operator(getFeSpace(2+i), getFeSpace(2+i));
+ opUGradU->addTerm(new WorldVec_FOT(stage_velocity, -1.0), GRD_PHI);
+ opUGradU->setUhOld(stage_velocity[i]);
+ prob->addVectorOperator(*opUGradU, 2+i);
+
+ /// < u_old*grad(u_i) , psi >
+ Operator *opUGradV = new Operator(getFeSpace(2+i), getFeSpace(2+i));
+ opUGradV->addTerm(new WorldVec_FOT(un_velocity, -1.0), GRD_PHI);
+ prob->addMatrixOperator(*opUGradV, 2+i, 2+i, &minus1, &minus1);
+
+ /// < u*grad(u_old_i) , psi >
+ for (size_t j = 0; j < dow; ++j) {
+ Operator *opVGradU = new Operator(getFeSpace(2+i), getFeSpace(2+j));
+ opVGradU->addTerm(new PartialDerivative_ZOT(un_velocity[i], j, -1.0));
+ prob->addMatrixOperator(*opVGradU, 2+i, 2+j, &minus1, &minus1);
+ }
+
+ /// Diffusion-Operator (including Stress-Tensor for space-dependent viscosity)
+ addLaplaceTerm(i);
+
+ /// < p , d_i(psi) >
+ Operator *opGradP = new Operator(getFeSpace(2+i),getFeSpace(2+dow));
+ opGradP->addTerm(new PartialDerivative_FOT(i), GRD_PSI);
+ opGradP->setUhOld(prob->getStageSolution(2+dow));
+ prob->addVectorOperator(*opGradP, 2+i);
+ prob->addMatrixOperator(*opGradP, 2+i, 2+dow, &minus1, &minus1);
+
+ /// external force, i.e. gravitational force
+ if (norm(force) > DBL_TOL) {
+ Operator *opForce = new Operator(getFeSpace(2+i), getFeSpace(2+i));
+ opForce->addTerm(new Simple_ZOT(force[i]));
+ prob->addVectorOperator(*opForce, 2+i);
+ }
+
+
+ // forces by Cahn-Hilliard terms
+ // -----------------------------
+ /// < 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);
+ }
+
+ // div(u) = 0
+ // ----------
+ for (size_t i = 0; i < dow; ++i) {
+ /// < d_i(u_i) , psi >
+ Operator *opDivU = new Operator(getFeSpace(2+dow),getFeSpace(2+i));
+ opDivU->addTerm(new PartialDerivative_FOT(i), GRD_PHI);
+ prob->addMatrixOperator(*opDivU, 2+dow, 2+i, &minus1, &minus1);
+ /// < d_i(u_s_i) , psi >
+ opDivU->setUhOld(stage_velocity[i]);
+ prob->addVectorOperator(*opDivU, 2+dow);
+ }
+}
+
+
+void CahnHilliardNavierStokes_RB::addLaplaceTerm(int i)
+{ FUNCNAME("CahnHilliardNavierStokes_RB::addLaplaceTerm()");
+
+ /// < nu*[grad(u)+grad(u)^t] , grad(psi) >
+ if (laplaceType == 1) {
+ for (unsigned j = 0; j < dow; ++j) {
+ Operator *opLaplaceUi1 = new Operator(getFeSpace(2+i), getFeSpace(2+j));
+ opLaplaceUi1->addTerm(new MatrixIJ_SOT(1-i, 1-j, viscosity));
+ prob->addMatrixOperator(*opLaplaceUi1, 2+i, 2+j);
+ opLaplaceUi1->setUhOld(prob->getStageSolution(2+j));
+ prob->addVectorOperator(*opLaplaceUi1, 2+i, &minus1, &minus1);
+ }
+ }
+
+ /// < nu*grad(u'_i) , grad(psi) >
+ Operator *opLaplaceUi = new Operator(getFeSpace(2+i), getFeSpace(2+i));
+ opLaplaceUi->addTerm(new Simple_SOT(viscosity));
+ prob->addMatrixOperator(*opLaplaceUi, 2+i, 2+i);
+ opLaplaceUi->setUhOld(prob->getStageSolution(2+i));
+ prob->addVectorOperator(*opLaplaceUi, 2+i, &minus1, &minus1);
+}
+
+void CahnHilliardNavierStokes_RB::closeTimestep(AdaptInfo *adaptInfo)
+{ FUNCNAME("CahnHilliardNavierStokes_RB::closeTimestep()");
+
+ calcVelocity();
+ fileWriter->writeFiles(adaptInfo, false);
+ writeFiles(adaptInfo, false);
+}
diff --git a/extensions/base_problems/CahnHilliardNavierStokes_RB.h b/extensions/base_problems/CahnHilliardNavierStokes_RB.h
new file mode 100644
index 00000000..bcab1017
--- /dev/null
+++ b/extensions/base_problems/CahnHilliardNavierStokes_RB.h
@@ -0,0 +1,272 @@
+/** \file CahnHilliardNavierStokes_RB.h */
+
+#ifndef CAHN_HILLIARD_H
+#define CAHN_HILLIARD_H
+
+#include "AMDiS.h"
+#include "BaseProblem_RB.h"
+#include "HL_SignedDistTraverse.h"
+
+using namespace AMDiS;
+
+class CahnHilliardNavierStokes_RB : public BaseProblem_RB
+{
+public: // definition of types
+
+ typedef BaseProblem_RB super;
+
+public: // public methods
+
+ CahnHilliardNavierStokes_RB(const std::string &name_);
+ ~CahnHilliardNavierStokes_RB();
+
+ virtual void initData();
+
+ virtual void transferInitialSolution(AdaptInfo *adaptInfo);
+ virtual void closeTimestep(AdaptInfo *adaptInfo);
+
+ double getEpsilon() { return eps; }
+ int getDoubleWellType() { return doubleWell; }
+
+ DOFVector<WorldVector<double> >* getVelocity()
+ {
+ return velocity;
+ }
+
+ virtual void fillOperators();
+ virtual void fillBoundaryConditions() {}
+ virtual void addLaplaceTerm(int i);
+
+protected: // protected variables
+
+ HL_SignedDistTraverse *reinit;
+
+ DOFVector<WorldVector<double> >* velocity;
+
+ void calcVelocity()
+ {
+ if (dow == 1)
+ transformDOF(prob->getSolution()->getDOFVector(2), velocity, new AMDiS::Vec1WorldVec<double>);
+ else if (dow == 2)
+ transformDOF(prob->getSolution()->getDOFVector(2), prob->getSolution()->getDOFVector(3), velocity, new AMDiS::Vec2WorldVec<double>);
+ else if (dow == 3)
+ transformDOF(prob->getSolution()->getDOFVector(2), prob->getSolution()->getDOFVector(3), prob->getSolution()->getDOFVector(4), velocity, new AMDiS::Vec3WorldVec<double>);
+ }
+
+ bool useMobility;
+ bool useConservationForm;
+
+ unsigned dim;
+
+ int laplaceType;
+ int nonLinTerm;
+ int doubleWell;
+
+ // navierStokes parameters
+ double oldTimestep;
+ double viscosity;
+ double density;
+ double c;
+
+ // Cahn-Hilliard parameters
+ double gamma;
+ double eps;
+ double minusEps;
+ double epsInv;
+ double minusEpsInv;
+ double epsSqr;
+ double minusEpsSqr;
+
+ double sigma; // coupling parameter to calculate the surface tension
+ double surfaceTension;// := sigma/epsilon
+
+ WorldVector<double> force;
+
+ FileVectorWriter *fileWriter;
+};
+
+
+/** \brief
+ * Abstract function for Cahn-Hilliard mobility
+ */
+class MobilityCH0 : public AbstractFunction<double,double>
+{
+ public:
+ MobilityCH0(double gamma_=1.0) :
+ AbstractFunction<double,double>(4),
+ gamma(gamma_),
+ delta(1.e-6) { }
+ double operator()(const double &ch) const
+ {
+ double phase = std::max(0.0, std::min(1.0, ch));
+ double mobility = 0.25*sqr(phase)*sqr(phase-1.0);
+ return gamma * std::max(mobility, delta);
+ }
+
+ protected:
+ double gamma;
+ double delta;
+};
+
+class MobilityCH1 : public AbstractFunction<double,double>
+{
+ public:
+ MobilityCH1(double gamma_=1.0) :
+ AbstractFunction<double,double>(4),
+ gamma(gamma_),
+ delta(1.e-6) { }
+ double operator()(const double &ch) const
+ {
+ double phase = std::max(-1.0, std::min(1.0, ch));
+ double mobility = 0.25*sqr(sqr(phase)-1.0);
+ return gamma * std::max(mobility, delta);
+ }
+
+ protected:
+ double gamma;
+ double delta;
+};
+
+class ScaleFunctor : public AbstractFunction<double,double>
+{
+ public:
+ ScaleFunctor(double factor_=1.0) :
+ AbstractFunction<double,double>(1),
+ factor(factor_) { }
+ double operator()(const double &ch) const
+ {
+ return factor * ch;
+ }
+
+ protected:
+ double factor;
+};
+
+class Pow2Functor : public AbstractFunction<double,double>
+{
+ public:
+ Pow2Functor(double factor_=1.0) :
+ AbstractFunction<double,double>(2),
+ factor(factor_) { }
+ double operator()(const double &ch) const
+ {
+ return factor * sqr(ch);
+ }
+
+ protected:
+ double factor;
+};
+
+class Pow3Functor : public AbstractFunction<double,double>
+{
+ public:
+ Pow3Functor(double factor_=1.0) :
+ AbstractFunction<double,double>(3),
+ factor(factor_) { }
+ double operator()(const double &ch) const
+ {
+ return factor * sqr(ch) * ch;
+ }
+
+ protected:
+ double factor;
+};
+
+
+class MobilityCH0Diff : public BinaryAbstractFunction<WorldVector<double>,double,WorldVector<double> >
+{
+ public:
+ MobilityCH0Diff(double factor_=1.0) :
+ BinaryAbstractFunction<WorldVector<double>,double,WorldVector<double> >(4),
+ factor(factor_) { }
+ WorldVector<double> operator()(const double &c, const WorldVector<double>& grdMu) const
+ {
+ double phase = std::max(0.0, std::min(1.0, c));
+ return (factor * (sqr(phase)*phase - 1.5*sqr(phase) + 0.5*phase)) * grdMu;
+ }
+
+ protected:
+ double factor;
+};
+
+class MobilityCH1Diff : public BinaryAbstractFunction<WorldVector<double>,double,WorldVector<double> >
+{
+ public:
+ MobilityCH1Diff(double factor_=1.0) :
+ BinaryAbstractFunction<WorldVector<double>,double,WorldVector<double> >(4),
+ factor(factor_) { }
+ WorldVector<double> operator()(const double &c, const WorldVector<double>& grdMu) const
+ {
+ double phase = std::max(-1.0, std::min(1.0, c));
+ return (factor * (sqr(phase)*phase - phase)) * grdMu;
+ }
+
+ protected:
+ double factor;
+};
+
+class DoubleWell0 : public AbstractFunction<double,double>
+{
+ public:
+ DoubleWell0(double factor_=1.0) :
+ AbstractFunction<double,double>(4),
+ factor(factor_) { }
+
+ double operator()(const double &c) const
+ {
+ return factor * (sqr(c)*c - 1.5*sqr(c) + 0.5*c);
+ }
+
+ protected:
+ double factor;
+};
+
+class DoubleWell1 : public AbstractFunction<double,double>
+{
+ public:
+ DoubleWell1(double factor_=1.0) :
+ AbstractFunction<double,double>(4),
+ factor(factor_) { }
+
+ double operator()(const double &c) const
+ {
+ return factor * (sqr(c)*c - c);
+ }
+
+ protected:
+ double factor;
+};
+
+class DoubleWell0Diff : public AbstractFunction<double,double>
+{
+ public:
+ DoubleWell0Diff(double factor_=1.0) :
+ AbstractFunction<double,double>(4),
+ factor(factor_) { }
+
+ double operator()(const double &c) const
+ {
+ return factor * (3.0*sqr(c) - 3.0*c + 0.5);
+ }
+
+ protected:
+ double factor;
+};
+
+class DoubleWell1Diff : public AbstractFunction<double,double>
+{
+ public:
+ DoubleWell1Diff(double factor_=1.0) :
+ AbstractFunction<double,double>(4),
+ factor(factor_) { }
+
+ double operator()(const double &c) const
+ {
+ return factor * (3.0*sqr(c) - 1.0);
+ }
+
+ protected:
+ double factor;
+};
+
+#endif // CAHN_HILLIARD_H
diff --git a/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase_RB.cc b/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase_RB.cc
new file mode 100644
index 00000000..38b0df01
--- /dev/null
+++ b/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase_RB.cc
@@ -0,0 +1,477 @@
+#include "CahnHilliardNavierStokes_TwoPhase_RB.h"
+#include "Views.h"
+#include "SignedDistFunctors.h"
+#include "PhaseFieldConvert.h"
+#include "POperators.h"
+
+using namespace AMDiS;
+
+CahnHilliardNavierStokes_TwoPhase_RB::CahnHilliardNavierStokes_TwoPhase_RB(const std::string &name_) :
+ super(name_),
+ useMobility(false),
+ useConservationForm(false),
+ doubleWell(1),
+ gamma(1.0),
+ eps(0.1),
+ minusEps(-0.1),
+ epsInv(10.0),
+ minusEpsInv(-10.0),
+ epsSqr(0.01),
+ minusEpsSqr(-0.01),
+ laplaceType(0),
+ nonLinTerm(3),
+ oldTimestep(0.0),
+ viscosity1(1.0),
+ viscosity2(1.0),
+ density1(1.0),
+ density2(1.0),
+ c(0.0),
+ sigma(1.0),
+ surfaceTension(1.0),
+ reinit(NULL),
+ velocity(NULL),
+ fileWriter(NULL)
+{
+ // parameters for CH
+ Parameters::get(name + "->use mobility", useMobility); // dt(c) = div(M(c)*grad(c)) + ...
+ Parameters::get(name + "->gamma", gamma); // M(c) = gamma * Double-Well
+ 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 for navier-stokes
+ Initfile::get(name + "->viscosity1", viscosity1);
+ Initfile::get(name + "->viscosity2", viscosity2);
+ Initfile::get(name + "->density1", density1);
+ Initfile::get(name + "->density2", density2);
+
+ // type of laplace operator: 0... div(nu*grad(u)),
+ // 1... div(0.5*nu*(grad(u)+grad(u)^T))
+ Initfile::get(name + "->laplace operator", laplaceType);
+
+ // Parameters for CH-Coupling
+ // 0... u*grad(c), 1... -div(u*c)
+ 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;
+
+
+ force.set(0.0);
+ Initfile::get(name + "->force", force);
+
+ // transformation of the parameters
+ minusEps = -eps;
+ epsInv = 1.0/eps;
+ minusEpsInv = -epsInv;
+ epsSqr = sqr(eps);
+ minusEpsSqr = -epsSqr;
+}
+
+
+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;
+ }
+ delete fileWriter;
+ fileWriter = NULL;
+}
+
+
+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");
+ fileWriter = new FileVectorWriter(name + "->velocity->output", getFeSpace(2)->getMesh(), velocity);
+
+ super::initData();
+}
+
+
+void CahnHilliardNavierStokes_TwoPhase_RB::transferInitialSolution(AdaptInfo *adaptInfo)
+{ FUNCNAME("CahnHilliardNavierStokes_TwoPhase_RB::transferInitialSolution()");
+
+ calcVelocity();
+
+ for (int i = 0; i < 2+dow; i++)
+ prob->setExactSolution(prob->getSolution()->getDOFVector(i), i);
+
+ fileWriter->writeFiles(adaptInfo, false);
+ writeFiles(adaptInfo, false);
+
+ // initial parameters for detecting mesh changes
+ oldMeshChangeIdx= getMesh()->getChangeIndex();
+}
+
+
+void CahnHilliardNavierStokes_TwoPhase_RB::fillOperators()
+{ FUNCNAME("CahnHilliardNavierStokes_TwoPhase_RB::fillOperators()");
+ MSG("CahnHilliardNavierStokes_TwoPhase_RB::fillOperators()\n");
+
+ // variable order:
+ // (c, mu, u0, u1 [, u2], p)
+
+ 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);
+ }
+
+ // dt(c) = laplace(mu) - u*grad(c)
+ // -----------------------------------
+ /// < dt(c) , psi >
+ prob->addTimeOperator(0, 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)));
+ opChLM->addTerm(new VecGrad_FOT(
+ prob->getUnVec(0),
+ prob->getUnVec(1),
+ new MobilityCH0Diff(gamma)), GRD_PSI);
+ prob->addMatrixOperator(*opChLM, 0, 1);
+ // rhs operators
+ Operator *opChLM2 = new Operator(prob->getFeSpace(0),prob->getFeSpace(1));
+ opChLM2->addTerm(new VecAtQP_SOT(
+ prob->getStageSolution(0),
+ new MobilityCH0(gamma)));
+ opChLM2->setUhOld(prob->getStageSolution(1));
+ prob->addVectorOperator(*opChLM2, 0, &minus1, &minus1);
+ } else {
+ // jacobian operators
+ opChLM->addTerm(new VecAtQP_SOT(
+ prob->getUnVec(0),
+ new MobilityCH1(gamma)));
+ opChLM->addTerm(new VecGrad_FOT(
+ prob->getUnVec(0),
+ prob->getUnVec(1),
+ new MobilityCH1Diff(gamma)), GRD_PSI);
+ prob->addMatrixOperator(*opChLM, 0, 1);
+ // rhs operators
+ Operator *opChLM2 = new Operator(prob->getFeSpace(0),prob->getFeSpace(1));
+ opChLM2->addTerm(new VecAtQP_SOT(
+ prob->getStageSolution(0),
+ new MobilityCH1(gamma)));
+ opChLM2->setUhOld(prob->getStageSolution(1));
+ prob->addVectorOperator(*opChLM2, 0, &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) >
+ }
+
+ /// < u_s * grad(c_s), psi >
+ Operator *uGradC = new Operator(prob->getFeSpace(0), prob->getFeSpace(0));
+ if (useConservationForm)
+ uGradC->addTerm(new WorldVec_FOT(stage_velocity, -1.0), GRD_PSI);
+ else
+ uGradC->addTerm(new WorldVec_FOT(stage_velocity), GRD_PHI);
+ uGradC->setUhOld(prob->getStageSolution(0));
+ prob->addVectorOperator(*uGradC, 0, &minus1, &minus1);
+
+ /// < u_old * grad(c), psi >
+ Operator *JuGradC1 = new Operator(prob->getFeSpace(0), prob->getFeSpace(0));
+ if (useConservationForm)
+ 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);
+
+ /// < u * grad(c_old), psi >
+ for (size_t i = 0; i < dow; i++) {
+ Operator *JuGradC2 = new Operator(prob->getFeSpace(0), prob->getFeSpace(0));
+ if (useConservationForm)
+ JuGradC2->addTerm(new VecAndPartialDerivative_FOT(prob->getUnVec(0), i, -1.0), GRD_PSI);
+ 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);
+ }
+
+ // Navier-Stokes part
+ // -----------------------------------------------------------------------------------
+ for (size_t i = 0; i < dow; ++i) {
+ /// < d_t(u_i) , psi >
+ prob->addTimeOperator(3+dow+i, 2+i);
+
+ /// < z , psi >
+ Operator *opZ = new Operator(getFeSpace(3+dow+i), getFeSpace(3+dow+i));
+ opZ->addTerm(new Simple_ZOT);
+ prob->addMatrixOperator(*opZ, 3+dow+i, 3+dow+i, &minus1, &minus1);
+ opZ->setUhOld(prob->getStageSolution(3+dow+i));
+ prob->addVectorOperator(*opZ, 3+dow+i);
+
+
+ /// < u_s*grad(u_s_i) , psi >
+ Operator *opUGradU = new Operator(getFeSpace(2+i), getFeSpace(2+i));
+ opUGradU->addTerm(new WorldVec_FOT(stage_velocity, -1.0), GRD_PHI);
+ opUGradU->setUhOld(stage_velocity[i]);
+ prob->addVectorOperator(*opUGradU, 3+dow+i);
+
+ /// < u_old*grad(u_i) , psi >
+ Operator *opUGradV = new Operator(getFeSpace(2+i), getFeSpace(2+i));
+ opUGradV->addTerm(new WorldVec_FOT(un_velocity, -1.0), GRD_PHI);
+ prob->addMatrixOperator(*opUGradV, 3+dow+i, 2+i, &minus1, &minus1);
+
+ /// < u*grad(u_old_i) , psi >
+ for (size_t j = 0; j < dow; ++j) {
+ Operator *opVGradU = new Operator(getFeSpace(2+i), getFeSpace(2+j));
+ opVGradU->addTerm(new PartialDerivative_ZOT(un_velocity[i], j, -1.0));
+ prob->addMatrixOperator(*opVGradU, 3+dow+i, 2+j, &minus1, &minus1);
+ }
+ }
+
+ for (size_t i = 0; i < dow; ++i) {
+
+ /// Diffusion-Operator (including Stress-Tensor for space-dependent viscosity)
+ addStressTerm(i);
+
+ /// external force, i.e. gravitational force
+ if (norm(force) > DBL_TOL) {
+ Operator *opForce = new Operator(getFeSpace(2+i), getFeSpace(2+i));
+ opForce->addTerm(new Simple_ZOT(force[i]));
+ prob->addVectorOperator(*opForce, 2+i);
+ }
+
+ // forces by Cahn-Hilliard terms
+ // -----------------------------
+ /// < 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);
+
+ // ---------------------------------------
+
+ /// < rho(c_s)*z_s , psi >
+ Operator *opDensityExpl = new Operator(prob->getFeSpace(2+i),prob->getFeSpace(3+dow+1));
+ opDensityExpl->addZeroOrderTerm(new VecAtQP_ZOT(
+ prob->getStageSolution(0),
+ new LinearInterpolation(density1, density2, -1.0)));
+ opDensityExpl->setUhOld(prob->getStageSolution(3+dow+i));
+ prob->addVectorOperator(*opDensityExpl, 2+i);
+
+ /// < rho(c_old)*z , psi >
+ Operator *opDensityImpl = new Operator(prob->getFeSpace(1),prob->getFeSpace(0));
+ opDensityImpl->addZeroOrderTerm(new VecAtQP_ZOT(
+ prob->getUnVec(0),
+ new LinearInterpolation(density1, density2, -1.0)));
+ prob->addMatrixOperator(*opDensityImpl, 2+i, 3+dow+i, &minus1, &minus1);
+
+ /// < rho'(c_old)*z_old * c , psi >
+ Operator *opDensityImpl2 = new Operator(prob->getFeSpace(1),prob->getFeSpace(0));
+ opDensityImpl2->addZeroOrderTerm(new Vec2AtQP_ZOT(
+ prob->getUnVec(0),
+ prob->getUnVec(3+dow+i),
+ new LinearInterpolationDiffFactor(density1, density2, -1.0)));
+ prob->addMatrixOperator(*opDensityImpl2, 2+i, 0, &minus1, &minus1);
+ }
+
+ // div(u) = 0
+ // ----------
+ for (size_t i = 0; i < dow; ++i) {
+ /// < d_i(u_i) , psi >
+ Operator *opDivU = new Operator(getFeSpace(2+dow),getFeSpace(2+i));
+ opDivU->addTerm(new PartialDerivative_FOT(i), GRD_PHI);
+ prob->addMatrixOperator(*opDivU, 2+dow, 2+i, &minus1, &minus1);
+ /// < d_i(u_s_i) , psi >
+ opDivU->setUhOld(stage_velocity[i]);
+ prob->addVectorOperator(*opDivU, 2+dow);
+ }
+}
+
+
+// void CahnHilliardNavierStokes_TwoPhase_RB::addNonlinearTerm(int i)
+// { FUNCNAME("CahnHilliardNavierStokes_TwoPhase_RB::addNonlinearTerm()");
+//
+// WorldVector<DOFVector<double>*> stage_velocity;
+// WorldVector<DOFVector<double>*> un_velocity;
+// for (size_t j = 0; j < dow; ++j) {
+// stage_velocity[j] = prob->getStageSolution(2+j);
+// un_velocity[j] = prob->getUnVec(2+j);
+// }
+//
+// /// < rho(c_s)*u_s*grad(u_s_i) , psi >
+// Operator *opUGradU = new Operator(getFeSpace(2+i), getFeSpace(2+i));
+// opUGradU->addTerm(new VecAndWorldVec_FOT(
+// prob->getStageSolution(0),
+// stage_velocity,
+// new LinearInterpolation(density1, density2, -1.0)), GRD_PHI);
+// opUGradU->setUhOld(stage_velocity[i]);
+// prob->addVectorOperator(*opUGradU, 2+i);
+//
+// /// < rho(c_old)*u_old*grad(u_i) , psi >
+// Operator *opUGradV = new Operator(getFeSpace(2+i), getFeSpace(2+i));
+// opUGradV->addTerm(new VecAndWorldVec_FOT(
+// prob->getUnVec(0),
+// un_velocity,
+// new LinearInterpolation(density1, density2, -1.0)), GRD_PHI);
+// prob->addMatrixOperator(*opUGradV, 2+i, 2+i, &minus1, &minus1);
+//
+// /// < c * rho'(c_old)*u_old*grad(u_old_i) , psi >
+// for (size_t j = 0; j < dow; ++j) {
+// Operator *opVGradU = new Operator(getFeSpace(2+i), getFeSpace(2+j));
+// opVGradU->addTerm(new Vec2AndPartialDerivative_ZOT(
+// prob->getUnVec(0),
+// un_velocity[i],
+// un_velocity[j],
+// j,
+// new LinearInterpolationDiffFactor(density1, density2, -1.0)));
+// prob->addMatrixOperator(*opVGradU, 2+i, 0, &minus1, &minus1);
+// }
+// }
+
+
+void CahnHilliardNavierStokes_TwoPhase_RB::addStressTerm(int i)
+{ FUNCNAME("CahnHilliardNavierStokes_TwoPhase_RB::addStressTerm()");
+
+ /// < nu*[grad(u)+grad(u)^t] , grad(psi) >
+ if (laplaceType == 1) {
+ for (unsigned j = 0; j < dow; ++j) {
+ /// < nu(c_s)*d_j(u_s_i) , d_i(psi) >
+ Operator *opLaplaceUijExpl = new Operator(getFeSpace(2+i), getFeSpace(2+j));
+ opLaplaceUijExpl->addTerm(new VecAtQP_IJ_SOT(
+ prob->getStageSolution(0),
+ new LinearInterpolation(viscosity1, viscosity2),
+ 1-i, 1-j));
+ opLaplaceUijExpl->setUhOld(prob->getStageSolution(2+j));
+ prob->addVectorOperator(*opLaplaceUijExpl, 2+i, &minus1, &minus1);
+
+ /// < nu(c_old)*d_j(u_i) , d_i(psi) >
+ Operator *opLaplaceUijImpl = new Operator(getFeSpace(2+i), getFeSpace(2+j));
+ opLaplaceUijImpl->addTerm(new VecAtQP_IJ_SOT(
+ prob->getUnVec(0),
+ new LinearInterpolation(viscosity1, viscosity2),
+ 1-i, 1-j));
+ prob->addMatrixOperator(*opLaplaceUijImpl, 2+i, 2+j);
+
+ /// < c * nu'(c_old)*d_j(u_old_i) , d_i(psi) >
+ Operator *opLaplaceUijImpl2 = new Operator(getFeSpace(2+i), getFeSpace(2+j));
+ opLaplaceUijImpl2->addTerm(new VecAndPartialDerivativeIJ_FOT(
+ prob->getUnVec(0),
+ prob->getUnVec(2+j),
+ 1-i, 1-j,
+ new LinearInterpolationDiffFactor(viscosity1, viscosity2)), GRD_PSI);
+ prob->addMatrixOperator(*opLaplaceUijImpl2, 2+i, 0);
+ }
+ }
+
+ /// < nu*grad(u'_i) , grad(psi) >
+ Operator *opLaplaceUExpl = new Operator(getFeSpace(2+i), getFeSpace(2+i));
+ opLaplaceUExpl->addTerm(new VecAtQP_SOT(
+ prob->getStageSolution(0),
+ new LinearInterpolation(viscosity1, viscosity2)));
+ opLaplaceUExpl->setUhOld(prob->getStageSolution(2+i));
+ prob->addVectorOperator(*opLaplaceUExpl, 2+i, &minus1, &minus1);
+
+ Operator *opLaplaceUImpl = new Operator(getFeSpace(2+i), getFeSpace(2+i));
+ opLaplaceUImpl->addTerm(new VecAtQP_SOT(
+ prob->getUnVec(0),
+ new LinearInterpolation(viscosity1, viscosity2)));
+ prob->addMatrixOperator(*opLaplaceUImpl, 2+i, 2+i);
+
+ Operator *opLaplaceUImpl2 = new Operator(getFeSpace(2+i), getFeSpace(2+i));
+ opLaplaceUImpl2->addTerm(new VecGrad_FOT(
+ prob->getUnVec(0),
+ prob->getUnVec(2+i),
+ new LinearInterpolationDiffVector(viscosity1, viscosity2)), GRD_PSI);
+ prob->addMatrixOperator(*opLaplaceUImpl2, 2+i, 0);
+
+ /// < p , d_i(psi) >
+ Operator *opGradP = new Operator(getFeSpace(2+i),getFeSpace(2+dow));
+ opGradP->addTerm(new PartialDerivative_FOT(i), GRD_PSI);
+ opGradP->setUhOld(prob->getStageSolution(2+dow));
+ prob->addVectorOperator(*opGradP, 2+i);
+ prob->addMatrixOperator(*opGradP, 2+i, 2+dow, &minus1, &minus1);
+}
+
+void CahnHilliardNavierStokes_TwoPhase_RB::closeTimestep(AdaptInfo *adaptInfo)
+{ FUNCNAME("CahnHilliardNavierStokes_TwoPhase_RB::closeTimestep()");
+
+ calcVelocity();
+ fileWriter->writeFiles(adaptInfo, false);
+ writeFiles(adaptInfo, false);
+}
diff --git a/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase_RB.h b/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase_RB.h
new file mode 100644
index 00000000..a4ccf450
--- /dev/null
+++ b/extensions/base_problems/CahnHilliardNavierStokes_TwoPhase_RB.h
@@ -0,0 +1,326 @@
+/** \file CahnHilliardNavierStokes_TwoPhase_RB.h */
+
+#ifndef CAHN_HILLIARD_H
+#define CAHN_HILLIARD_H
+
+#include "AMDiS.h"
+#include "BaseProblem_RB.h"
+#include "HL_SignedDistTraverse.h"
+
+using namespace AMDiS;
+
+class CahnHilliardNavierStokes_TwoPhase_RB : public BaseProblem_RB
+{
+public: // definition of types
+
+ typedef BaseProblem_RB super;
+
+public: // public methods
+
+ CahnHilliardNavierStokes_TwoPhase_RB(const std::string &name_);
+ ~CahnHilliardNavierStokes_TwoPhase_RB();
+
+ virtual void initData();
+
+ virtual void transferInitialSolution(AdaptInfo *adaptInfo);
+ virtual void closeTimestep(AdaptInfo *adaptInfo);
+
+ double getEpsilon() { return eps; }
+ int getDoubleWellType() { return doubleWell; }
+
+ DOFVector<WorldVector<double> >* getVelocity()
+ {
+ return velocity;
+ }
+
+ virtual void fillOperators();
+ virtual void fillBoundaryConditions() {}
+// virtual void addNonlinearTerm(int i);
+ virtual void addStressTerm(int i);
+
+protected: // protected variables
+
+ HL_SignedDistTraverse *reinit;
+
+ DOFVector<WorldVector<double> >* velocity;
+
+ void calcVelocity()
+ {
+ if (dow == 1)
+ transformDOF(prob->getSolution()->getDOFVector(2), velocity, new AMDiS::Vec1WorldVec<double>);
+ else if (dow == 2)
+ transformDOF(prob->getSolution()->getDOFVector(2), prob->getSolution()->getDOFVector(3), velocity, new AMDiS::Vec2WorldVec<double>);
+ else if (dow == 3)
+ transformDOF(prob->getSolution()->getDOFVector(2), prob->getSolution()->getDOFVector(3), prob->getSolution()->getDOFVector(4), velocity, new AMDiS::Vec3WorldVec<double>);
+ }
+
+ bool useMobility;
+ bool useConservationForm;
+
+ unsigned dim;
+
+ int laplaceType;
+ int nonLinTerm;
+ int doubleWell;
+
+ // navierStokes parameters
+ double oldTimestep;
+ double viscosity1;
+ double viscosity2;
+ double density1;
+ double density2;
+ double c;
+
+ // Cahn-Hilliard parameters
+ double gamma;
+ double eps;
+ double minusEps;
+ double epsInv;
+ double minusEpsInv;
+ double epsSqr;
+ double minusEpsSqr;
+
+ double sigma; // coupling parameter to calculate the surface tension
+ double surfaceTension;// := sigma/epsilon
+
+ WorldVector<double> force;
+
+ FileVectorWriter *fileWriter;
+};
+
+
+/** \brief
+ * Abstract function for Cahn-Hilliard mobility
+ */
+class MobilityCH0 : public AbstractFunction<double,double>
+{
+ public:
+ MobilityCH0(double gamma_=1.0) :
+ AbstractFunction<double,double>(4),
+ gamma(gamma_),
+ delta(1.e-6) { }
+ double operator()(const double &ch) const
+ {
+ double phase = std::max(0.0, std::min(1.0, ch));
+ double mobility = 0.25*sqr(phase)*sqr(phase-1.0);
+ return gamma * std::max(mobility, delta);
+ }
+
+ protected:
+ double gamma;
+ double delta;
+};
+
+class MobilityCH1 : public AbstractFunction<double,double>
+{
+ public:
+ MobilityCH1(double gamma_=1.0) :
+ AbstractFunction<double,double>(4),
+ gamma(gamma_),
+ delta(1.e-6) { }
+ double operator()(const double &ch) const
+ {
+ double phase = std::max(-1.0, std::min(1.0, ch));
+ double mobility = 0.25*sqr(sqr(phase)-1.0);
+ return gamma * std::max(mobility, delta);
+ }
+
+ protected:
+ double gamma;
+ double delta;
+};
+
+class ScaleFunctor : public AbstractFunction<double,double>
+{
+ public:
+ ScaleFunctor(double factor_=1.0) :
+ AbstractFunction<double,double>(1),
+ factor(factor_) { }
+ double operator()(const double &ch) const
+ {
+ return factor * ch;
+ }
+
+ protected:
+ double factor;
+};
+
+class Pow2Functor : public AbstractFunction<double,double>
+{
+ public:
+ Pow2Functor(double factor_=1.0) :
+ AbstractFunction<double,double>(2),
+ factor(factor_) { }
+ double operator()(const double &ch) const
+ {
+ return factor * sqr(ch);
+ }
+
+ protected:
+ double factor;
+};
+
+class Pow3Functor : public AbstractFunction<double,double>
+{
+ public:
+ Pow3Functor(double factor_=1.0) :
+ AbstractFunction<double,double>(3),
+ factor(factor_) { }
+ double operator()(const double &ch) const
+ {
+ return factor * sqr(ch) * ch;
+ }
+
+ protected:
+ double factor;
+};
+
+
+class MobilityCH0Diff : public BinaryAbstractFunction<WorldVector<double>,double,WorldVector<double> >
+{
+ public:
+ MobilityCH0Diff(double factor_=1.0) :
+ BinaryAbstractFunction<WorldVector<double>,double,WorldVector<double> >(4),
+ factor(factor_) { }
+ WorldVector<double> operator()(const double &c, const WorldVector<double>& grdMu) const
+ {
+ double phase = std::max(0.0, std::min(1.0, c));
+ return (factor * (sqr(phase)*phase - 1.5*sqr(phase) + 0.5*phase)) * grdMu;
+ }
+
+ protected:
+ double factor;
+};
+
+class MobilityCH1Diff : public BinaryAbstractFunction<WorldVector<double>,double,WorldVector<double> >
+{
+ public:
+ MobilityCH1Diff(double factor_=1.0) :
+ BinaryAbstractFunction<WorldVector<double>,double,WorldVector<double> >(4),
+ factor(factor_) { }
+ WorldVector<double> operator()(const double &c, const WorldVector<double>& grdMu) const
+ {
+ double phase = std::max(-1.0, std::min(1.0, c));
+ return (factor * (sqr(phase)*phase - phase)) * grdMu;
+ }
+
+ protected:
+ double factor;
+};
+
+class DoubleWell0 : public AbstractFunction<double,double>
+{
+ public:
+ DoubleWell0(double factor_=1.0) :
+ AbstractFunction<double,double>(4),
+ factor(factor_) { }
+
+ double operator()(const double &c) const
+ {
+ return factor * (sqr(c)*c - 1.5*sqr(c) + 0.5*c);
+ }
+
+ protected:
+ double factor;
+};
+
+class DoubleWell1 : public AbstractFunction<double,double>
+{
+ public:
+ DoubleWell1(double factor_=1.0) :
+ AbstractFunction<double,double>(4),
+ factor(factor_) { }
+
+ double operator()(const double &c) const
+ {
+ return factor * (sqr(c)*c - c);
+ }
+
+ protected:
+ double factor;
+};
+
+class DoubleWell0Diff : public AbstractFunction<double,double>
+{
+ public:
+ DoubleWell0Diff(double factor_=1.0) :
+ AbstractFunction<double,double>(4),
+ factor(factor_) { }
+
+ double operator()(const double &c) const
+ {
+ return factor * (3.0*sqr(c) - 3.0*c + 0.5);
+ }
+
+ protected:
+ double factor;
+};
+
+class DoubleWell1Diff : public AbstractFunction<double,double>
+{
+ public:
+ DoubleWell1Diff(double factor_=1.0) :
+ AbstractFunction<double,double>(4),
+ factor(factor_) { }
+
+ double operator()(const double &c) const
+ {
+ return factor * (3.0*sqr(c) - 1.0);
+ }
+
+ protected:
+ double factor;
+};
+
+struct LinearInterpolation : AbstractFunction<double,double>
+{
+ LinearInterpolation(double v0_, double v1_, double factor_ = 1.0) :
+ AbstractFunction<double,double>(4),
+ v0(v0_), v1(v1_), factor(factor_) { }
+
+ double operator()(const double &c) const
+ {
+ return factor * (v0*(c+1.0) - v1*(c-1.0))/2.0;
+ }
+
+private:
+ double v0;
+ double v1;
+ double factor;
+};
+
+struct LinearInterpolationDiffFactor : BinaryAbstractFunction<double,double,double>
+{
+ LinearInterpolationDiffFactor(double v0_, double v1_, double factor_ = 1.0) :
+ BinaryAbstractFunction<double,double,double>(4),
+ v0(v0_), v1(v1_), factor(factor_) { }
+
+ double operator()(const double &c, const double& z) const
+ {
+ return factor * z * (v0 - v1)/2.0;
+ }
+
+private:
+ double v0;
+ double v1;
+ double factor;
+};
+
+struct LinearInterpolationDiffVector : BinaryAbstractFunction<WorldVector<double>,double,WorldVector<double> >
+{
+ LinearInterpolationDiffVector(double v0_, double v1_, double factor_ = 1.0) :
+ BinaryAbstractFunction<WorldVector<double>,double,WorldVector<double> >(4),
+ v0(v0_), v1(v1_), factor(factor_) { }
+
+ WorldVector<double> operator()(const double &c, const WorldVector<double>& grdU) const
+ {
+ return (factor * (v0 - v1)/2.0) * grdU;
+ }
+
+private:
+ double v0;
+ double v1;
+ double factor;
+};
+
+#endif // CAHN_HILLIARD_H
diff --git a/extensions/base_problems/DiffuseDomainFsi.cc b/extensions/base_problems/DiffuseDomainFsi.cc
index 0b034703..247c8508 100644
--- a/extensions/base_problems/DiffuseDomainFsi.cc
+++ b/extensions/base_problems/DiffuseDomainFsi.cc
@@ -6,6 +6,7 @@ using namespace AMDiS;
DiffuseDomainFsi::DiffuseDomainFsi(const std::string &name_) :
super(name_),
phase(NULL),
+ bcPhase(NULL),
phaseInv(NULL),
velocity(NULL),
displacement(NULL),
@@ -17,22 +18,31 @@ DiffuseDomainFsi::DiffuseDomainFsi(const std::string &name_) :
density(1.0),
lambda(1.0),
rho(1.0),
+ nonLinTerm(3),
fileWriterPhase(NULL),
fileWriterVelocity(NULL),
fileWriterDisplacement(NULL)
{
- Parameters::get(name + "->beta", beta);
Parameters::get(name + "->epsilon", epsilon);
- Parameters::get(name + "->alpha", alpha);
// parameters for navier-stokes
Initfile::get(name + "->viscosity", viscosity);
Initfile::get(name + "->density", density);
+
+ Parameters::get(name + "->beta", beta);
+ Parameters::get(name + "->alpha", alpha);
+
// parameters for linear elasticity
Parameters::get(name + "->mu", mu);
Parameters::get(name + "->lambda", lambda);
Parameters::get(name + "->rho", rho);
+
+ // type of non-linear term: 0... u^old*grad(u_i^old),
+ // 1... u*grad(u_i^old),
+ // 2... u^old*grad(u_i)
+ // 3... u*grad(u_i^old) + u^old*grad(u_i) - u^old*grad(u_i^old)
+ Parameters::get(name + "->non-linear term", nonLinTerm);
};
DiffuseDomainFsi::~DiffuseDomainFsi()
@@ -43,20 +53,26 @@ DiffuseDomainFsi::~DiffuseDomainFsi()
delete fileWriterVelocity;
if (fileWriterDisplacement != NULL)
delete fileWriterDisplacement;
+
+ delete phaseInv;
+ delete phaseInvDisturbed;
+ delete velocity;
+ delete displacement;
}
void DiffuseDomainFsi::initData()
{
super::initData();
- phaseInv = new DOFVector<double>(getFeSpace(dow+1), "phaseInv");
+ phaseInv = new DOFVector<double>(getFeSpace(0), "phaseInv");
+ phaseInvDisturbed = new DOFVector<double>(getFeSpace(0), "phaseInvDisturbed");
velocity = new DOFVector<WorldVector<double> >(getFeSpace(0), "velocity");
displacement = new DOFVector<WorldVector<double> >(getFeSpace(dow+1), "displacement");
dbc_factor = beta/pow(epsilon, alpha);
MSG_DBG("dbc_factor = %f\n", dbc_factor);
pressure_factor = -(lambda+(2.0/3.0)*mu);
- MSG_DBG("pressure_factor = %f\n", pressure_factor);
+ MSG("pressure_factor = %f\n", pressure_factor);
}
@@ -68,7 +84,6 @@ void DiffuseDomainFsi::transferInitialSolution(AdaptInfo *adaptInfo)
fileWriterDisplacement = new FileVectorWriter(name + "->displacement->output", getFeSpace()->getMesh(), displacement);
super::transferInitialSolution(adaptInfo);
- transformDOF(1.0, phase, phaseInv, new Subtract<double>);
calcVelocity();
calcDisplacement();
@@ -83,6 +98,15 @@ void DiffuseDomainFsi::transferInitialSolution(AdaptInfo *adaptInfo)
}
+void DiffuseDomainFsi::initTimestep(AdaptInfo *adaptInfo)
+{
+ super::initTimestep(adaptInfo);
+ transformDOF(1.0, phase, phaseInv, new Subtract<double>);
+ transformDOF(phaseInv, 1.e-7, phaseInvDisturbed, new Max<double>);
+ VtkWriter::writeFile(phaseInv, "phaseInv.vtu");
+}
+
+
void DiffuseDomainFsi::closeTimestep(AdaptInfo *adaptInfo)
{
calcVelocity();
@@ -108,47 +132,48 @@ void DiffuseDomainFsi::fillOperators()
/// < (1/tau)*u'_i , psi >
Operator *opTime = new Operator(prob->getFeSpace(i), prob->getFeSpace(i));
opTime->addTerm(new Simple_ZOT(density));
- prob->addMatrixOperator(*opTime, i, i, getInvTau());
+ prob->addMatrixOperator(*opTime, i, i, getInvTau(), getInvTau());
/// < (1/tau)*u_i^old , psi >
Operator *opTimeOld = new Operator(prob->getFeSpace(i), prob->getFeSpace(i));
- opTimeOld->addTerm(new VecAtQP_ZOT(prob->getSolution()->getDOFVector(i), new Factor<double>(density)));
- prob->addVectorOperator(*opTimeOld, i, getInvTau());
+ opTimeOld->addTerm(new Phase_ZOT(prob->getSolution()->getDOFVector(i), density));
+ prob->addVectorOperator(*opTimeOld, i, getInvTau(), getInvTau());
/// < (1/tau)*(1-phase)*eta_i , psi >
Operator *opTimeE = new Operator(getFeSpace(dow+1+i), getFeSpace(dow+1+i));
- opTimeE->addTerm(new VecAtQP_ZOT(phaseInv));
+ opTimeE->addTerm(new VecAtQP_ZOT(phaseInvDisturbed));
prob->addMatrixOperator(*opTimeE, dow+1+i, dow+1+i, getInvTau(), getInvTau());
/// < (1/tau)*(1-phase)*eta_i^old , psi >
Operator *opTimeEOld = new Operator(getFeSpace(i+dow));
- opTimeEOld->addTerm(new Vec2AtQP_ZOT(phaseInv, prob->getSolution()->getDOFVector(dow+1+i)));
+ opTimeEOld->addTerm(new Vec2AtQP_ZOT(phaseInvDisturbed, prob->getSolution()->getDOFVector(dow+1+i)));
prob->addVectorOperator(*opTimeEOld, dow+1+i, getInvTau(), getInvTau());
/// < (1-phase)*u , psi >
Operator *opU = new Operator(getFeSpace(dow+1+i), getFeSpace(i));
- opU->addTerm(new Phase_ZOT(phaseInv, -1.0));
- prob->addMatrixOperator(*opU, dow+1+i, dow+1+i);
+ opU->addTerm(new Phase_ZOT(phaseInvDisturbed, -1.0));
+ prob->addMatrixOperator(*opU, dow+1+i, i);
// nonlinear terms
// ===============
/// < u^old*grad(u_i^old) , psi >
- Operator *opUGradU0 = new Operator(getFeSpace(i), getFeSpace(i));
- opUGradU0->addTerm(new WorldVector_FOT(velocity, -density), GRD_PHI);
- opUGradU0->setUhOld(prob->getSolution()->getDOFVector(i));
- if (nonLinTerm == 0)
+ if (nonLinTerm == 0 || nonLinTerm == 3) {
+ Operator *opUGradU0 = new Operator(getFeSpace(i), getFeSpace(i));
+ opUGradU0->addTerm(new WorldVec_FOT(vel, (nonLinTerm == 0 ? -1.0 : 1.0)*density), GRD_PHI);
+ opUGradU0->setUhOld(prob->getSolution()->getDOFVector(i));
prob->addVectorOperator(*opUGradU0, i);
-
- if (nonLinTerm == 1) {
+ }
+ if (nonLinTerm == 1 || nonLinTerm == 3) {
/// < u'*grad(u_i^old) , psi >
- for (unsigned j = 0; j < dow; ++j) {
- Operator *opUGradU1 = new Operator(getFeSpace(i),getFeSpace(i));
+ for (size_t j = 0; j < dow; ++j) {
+ Operator *opUGradU1 = new Operator(getFeSpace(i),getFeSpace(j));
opUGradU1->addTerm(new PartialDerivative_ZOT(
prob->getSolution()->getDOFVector(i), j, density));
prob->addMatrixOperator(*opUGradU1, i, j);
}
- } else if (nonLinTerm == 2) {
+ }
+ if (nonLinTerm == 2 || nonLinTerm == 3) {
/// < u^old*grad(u'_i) , psi >
- for(unsigned j = 0; j < dow; ++j) {
+ for(size_t j = 0; j < dow; ++j) {
Operator *opUGradU2 = new Operator(getFeSpace(i),getFeSpace(i));
opUGradU2->addTerm(new VecAndPartialDerivative_FOT(
prob->getSolution()->getDOFVector(j), j, density), GRD_PHI);
@@ -158,32 +183,30 @@ void DiffuseDomainFsi::fillOperators()
// ====================
- /// < u^old*grad(u_i^old) , psi >
- Operator *opUGradE0 = new Operator(prob->getFeSpace(i),prob->getFeSpace(dow+1+i));
- opUGradE0->addTerm(new WorldVecPhase_FOT(phaseInv, vel), GRD_PHI);
- opUGradE0->setUhOld(prob->getSolution()->getDOFVector(dow+1+i));
- if (nonLinTerm == 0)
- prob->addVectorOperator(*opUGradE0, dow+1+i);
-
- if (nonLinTerm == 1) {
- /// < u'*grad(u_i^old) , psi >
- for (size_t j = 0; j < dow; ++j) {
- Operator *opUGradE1 = new Operator(prob->getFeSpace(i),prob->getFeSpace(dow+1+i));
- opUGradE1->addTerm(new VecAndPartialDerivative_ZOT(
- phaseInv,
- prob->getSolution()->getDOFVector(dow+1+i), j));
- prob->addMatrixOperator(*opUGradE1, dow+1+i, j);
- }
- } else if (nonLinTerm == 2) {
- /// < u^old*grad(u'_i) , psi >
- for(size_t j = 0; j < dow; ++j) {
- Operator *opUGradE2 = new Operator(prob->getFeSpace(i),prob->getFeSpace(i));
- opUGradE2->addTerm(new Vec2ProductPartial_FOT(
- phaseInv,
- prob->getSolution()->getDOFVector(j), j, density), GRD_PHI);
- prob->addMatrixOperator(*opUGradE2, dow+1+i, dow+1+i);
- }
- }
+ /// < (1-phase) * u^old*grad(eta_i^old) , psi >
+ Operator *opUGradE0 = new Operator(prob->getFeSpace(dow+1+i),prob->getFeSpace(dow+1+i));
+ opUGradE0->addTerm(new WorldVecPhase_FOT(phaseInvDisturbed, vel), GRD_PHI);
+ opUGradE0->setUhOld(prob->getSolution()->getDOFVector(dow+1+i));
+ prob->addVectorOperator(*opUGradE0, dow+1+i);
+
+ /// < (1-phase) * u*grad(eta_i^old) , psi >
+ for (size_t j = 0; j < dow; ++j) {
+ Operator *opUGradE1 = new Operator(prob->getFeSpace(dow+1+i),prob->getFeSpace(j));
+ opUGradE1->addTerm(new VecAndPartialDerivative_ZOT(
+ phaseInvDisturbed,
+ prob->getSolution()->getDOFVector(dow+1+i), j));
+ prob->addMatrixOperator(*opUGradE1, dow+1+i, j);
+ }
+
+ /// < (1-phase) * u^old*grad(eta_i) , psi >
+ for(size_t j = 0; j < dow; ++j) {
+ Operator *opUGradE2 = new Operator(prob->getFeSpace(dow+1+i),prob->getFeSpace(dow+1+i));
+ opUGradE2->addTerm(new Vec2ProductPartial_FOT(
+ phaseInvDisturbed,
+ prob->getSolution()->getDOFVector(j), j), GRD_PHI);
+ prob->addMatrixOperator(*opUGradE2, dow+1+i, dow+1+i);
+ }
+
/// Diffusion-Operator (including Stress-Tensor for space-dependent viscosity)
addLaplaceTerm(i);
@@ -207,39 +230,46 @@ void DiffuseDomainFsi::fillOperators()
opP->addTerm(new VecAtQP_ZOT(phaseInv));
prob->addMatrixOperator(*opP, dow, dow);
+ for (size_t i = 0; i < dow; ++i) {
+ /// dirichlet boundary condition
+ /// beta/eps^alpha (bcPhase)*(eta_i - 0)
+ Operator *opDbcPhaseInvU = new Operator(prob->getFeSpace(dow+1+i), prob->getFeSpace(dow+1+i));
+ opDbcPhaseInvU->addTerm(new Phase_ZOT(bcPhase, dbc_factor));
+ prob->addMatrixOperator(*opDbcPhaseInvU, dow+1+i, dow+1+i);
+ }
}
void DiffuseDomainFsi::addLaplaceTerm(int i)
{ FUNCNAME("DiffuseDomainFsi::addLaplaceTerm()");
- for (size_t j = 0; j < dow; ++j) {
- /// < nu*[grad(u)+grad(u)^t] , grad(psi) >
- Operator *opLaplaceUi1 = new Operator(prob->getFeSpace(i), prob->getFeSpace(j));
- opLaplaceUi1->addTerm(new MatrixIJPhase_SOT(phase, 1-i, 1-j, viscosity));
- prob->addMatrixOperator(*opLaplaceUi1, i, j);
-
- /// < (mu+lambda)*[grad(u)+grad(u)^t] , grad(psi) >
- Operator *opLaplaceEi1 = new Operator(prob->getFeSpace(dow+1+i), prob->getFeSpace(dow+1+j));
- opLaplaceEi1->addTerm(new MatrixIJPhase_SOT(phaseInv, 1-i, 1-j, mu+lambda));
- prob->addMatrixOperator(*opLaplaceEi1, dow+1+i, dow+1+j);
- }
-
- /// < nu*grad(u_i) , grad(psi) >
- Operator *opLaplaceUi = new Operator(prob->getFeSpace(i), prob->getFeSpace(i));
- opLaplaceUi->addTerm(new Phase_SOT(phase, viscosity));
- prob->addMatrixOperator(*opLaplaceUi, i, i);
-
- /// < mu*grad(eta_i) , grad(psi) >
- Operator *opLaplaceEi = new Operator(prob->getFeSpace(dow+1+i), prob->getFeSpace(dow+1+i));
- opLaplaceEi->addTerm(new Phase_SOT(phaseInv, mu));
- prob->addMatrixOperator(*opLaplaceEi, dow+1+i, dow+1+i);
-
- /// < p , d_i(psi) >
- Operator *opGradP = new Operator(prob->getFeSpace(i),prob->getFeSpace(dow));
- opGradP->addTerm(new VecAndPartialDerivative_FOT(phase,i, -1.0), GRD_PSI);
- prob->addMatrixOperator(*opGradP, i, dow);
- Operator *opGradP2 = new Operator(prob->getFeSpace(i),prob->getFeSpace(dow));
- opGradP2->addTerm(new PartialDerivative_ZOT(phase,i, -1.0));
- prob->addMatrixOperator(*opGradP2, i, dow);
+ for (size_t j = 0; j < dow; ++j) {
+ /// < nu*[grad(u)+grad(u)^t] , grad(psi) >
+ Operator *opLaplaceUi1 = new Operator(prob->getFeSpace(i), prob->getFeSpace(j));
+ opLaplaceUi1->addTerm(new MatrixIJPhase_SOT(phase, 1-i, 1-j, viscosity));
+ prob->addMatrixOperator(*opLaplaceUi1, i, j);
+
+ /// < (mu+lambda)*[grad(u)+grad(u)^t] , grad(psi) >
+ Operator *opLaplaceEi1 = new Operator(prob->getFeSpace(i), prob->getFeSpace(dow+1+j));
+ opLaplaceEi1->addTerm(new MatrixIJPhase_SOT(phaseInv, 1-i, 1-j, mu+lambda));
+ prob->addMatrixOperator(*opLaplaceEi1, i, dow+1+j);
+ }
+
+ /// < nu*grad(u_i) , grad(psi) >
+ Operator *opLaplaceUi = new Operator(prob->getFeSpace(i), prob->getFeSpace(i));
+ opLaplaceUi->addTerm(new Phase_SOT(phase, viscosity));
+ prob->addMatrixOperator(*opLaplaceUi, i, i);
+
+ /// < mu*grad(eta_i) , grad(psi) >
+ Operator *opLaplaceEi = new Operator(prob->getFeSpace(i), prob->getFeSpace(dow+1+i));
+ opLaplaceEi->addTerm(new Phase_SOT(phaseInv, mu));
+ prob->addMatrixOperator(*opLaplaceEi, i, dow+1+i);
+
+ /// < p , d_i(psi) >
+ Operator *opGradP = new Operator(prob->getFeSpace(i),prob->getFeSpace(dow));
+ opGradP->addTerm(new VecAndPartialDerivative_FOT(phase,i, -1.0), GRD_PSI);
+ prob->addMatrixOperator(*opGradP, i, dow);
+ Operator *opGradP2 = new Operator(prob->getFeSpace(i),prob->getFeSpace(dow));
+ opGradP2->addTerm(new PartialDerivative_ZOT(phase,i, -1.0));
+ prob->addMatrixOperator(*opGradP2, i, dow);
}
diff --git a/extensions/base_problems/DiffuseDomainFsi.h b/extensions/base_problems/DiffuseDomainFsi.h
index f54cb1a2..bcfdeb3a 100644
--- a/extensions/base_problems/DiffuseDomainFsi.h
+++ b/extensions/base_problems/DiffuseDomainFsi.h
@@ -32,6 +32,7 @@ public: // methods
virtual void initData();
virtual void transferInitialSolution(AdaptInfo *adaptInfo);
+ virtual void initTimestep(AdaptInfo *adaptInfo);
virtual void closeTimestep(AdaptInfo *adaptInfo);
// === getting/setting methods ===
@@ -61,6 +62,11 @@ public: // methods
phase = phase_;
}
+ void setBcPhase(DOFVector<double>* phase_)
+ {
+ bcPhase = phase_;
+ }
+
void setEpsilon(double epsilon_)
{
epsilon = epsilon_;
@@ -88,7 +94,6 @@ protected: // variables
double epsilon;
double alpha;
- int laplaceType;
int nonLinTerm;
double oldTimestep;
@@ -103,6 +108,8 @@ protected: // variables
DOFVector<double> *phase;
DOFVector<double> *phaseInv;
+ DOFVector<double> *phaseInvDisturbed;
+ DOFVector<double> *bcPhase;
DOFVector<WorldVector<double> >* velocity;
diff --git a/extensions/base_problems/LinearElasticity.cc b/extensions/base_problems/LinearElasticity.cc
index 05755a5b..2b3e8f3c 100644
--- a/extensions/base_problems/LinearElasticity.cc
+++ b/extensions/base_problems/LinearElasticity.cc
@@ -1,4 +1,5 @@
#include "LinearElasticity.h"
+#include "POperators.h"
using namespace std;
using namespace AMDiS;
@@ -8,6 +9,7 @@ LinearElasticity::LinearElasticity(const std::string &name_) :
mu(1.0),
lambda(1.0),
rho(1.0),
+ nonLinTerm(2),
displacement(NULL)
{
Parameters::get(name + "->mu", mu);
@@ -16,6 +18,9 @@ LinearElasticity::LinearElasticity(const std::string &name_) :
force.set(0.0);
Parameters::get(name + "->force", force);
+
+ // type of non-linear term: 0... u^old*grad(u_i^old), 1... u'*grad(u_i^old), 2... u^old*grad(u'_i)
+ Parameters::get(name + "->non-linear term", nonLinTerm);
}
@@ -69,6 +74,10 @@ void LinearElasticity::fillOperators()
{ FUNCNAME("LinearElasticity::fillOperators()");
MSG("LinearElasticity::fillOperators()\n");
+ WorldVector<DOFVector<double>* > velocity;
+ for (size_t i = 0; i < dow; i++)
+ velocity[i] = prob->getSolution()->getDOFVector(i+dow);
+
// fill operators for prob
for (size_t i = 0; i < dow; ++i) {
/// dt(u) - v = 0
@@ -86,7 +95,7 @@ void LinearElasticity::fillOperators()
}
// < v , psi >
Operator *opV = new Operator(getFeSpace(i+dow), getFeSpace(i+dow));
- opV->addTerm(new Simple_ZOT(abs(rho) > DBL_TOL ? -1.0/rho : -1.0));
+ opV->addTerm(new Simple_ZOT(-1.0));
prob->addMatrixOperator(*opV, i+dow, i+dow);
/// dt(v) - mu*div(grad(u)) - (mu+lambda)*grad(div(u)) = F
@@ -96,14 +105,43 @@ void LinearElasticity::fillOperators()
if (abs(rho) > DBL_TOL) {
// < (1/tau)*v_i , psi >
Operator *opTimeV = new Operator(getFeSpace(i), getFeSpace(i+dow));
- opTimeV->addTerm(new Simple_ZOT);
+ opTimeV->addTerm(new Simple_ZOT(rho));
prob->addMatrixOperator(*opTimeV, i, i+dow, getInvTau(), getInvTau());
// < (1/tau)*v_i^old , psi >
Operator *opTimeVOld = new Operator(getFeSpace(i));
- opTimeVOld->addTerm(new VecAtQP_ZOT(prob->getSolution()->getDOFVector(i+dow)));
+ opTimeVOld->addTerm(new Phase_ZOT(prob->getSolution()->getDOFVector(i+dow), rho));
prob->addVectorOperator(*opTimeVOld, i, getInvTau(), getInvTau());
}
+
+
+ /// < v^old*grad(v_i^old) , psi >
+ if (nonLinTerm == 0) {
+ Operator *opUGradU0 = new Operator(getFeSpace(i), getFeSpace(i+dow));
+ opUGradU0->addTerm(new WorldVec_FOT(velocity, -rho), GRD_PHI);
+ opUGradU0->setUhOld(prob->getSolution()->getDOFVector(i+dow));
+ prob->addVectorOperator(*opUGradU0, i);
+ }
+
+ if (nonLinTerm == 1) {
+ /// < v*grad(v_i^old) , psi >
+ for (size_t j = 0; j < dow; ++j) {
+ Operator *opUGradU1 = new Operator(getFeSpace(i),getFeSpace(j+dow));
+ opUGradU1->addTerm(new PartialDerivative_ZOT(
+ velocity[i], j, rho));
+ prob->addMatrixOperator(*opUGradU1, i, j);
+ }
+ } else if (nonLinTerm == 2) {
+ /// < v^old*grad(v_i) , psi >
+ for(size_t j = 0; j < dow; ++j) {
+ Operator *opUGradU2 = new Operator(getFeSpace(i),getFeSpace(i+dow));
+ opUGradU2->addTerm(new VecAndPartialDerivative_FOT(
+ velocity[j], j, rho), GRD_PHI);
+ prob->addMatrixOperator(*opUGradU2, i, i);
+ }
+ }
+
+
// < grad(u_i) , grad(psi) >
Operator *opDivGradU = new Operator(getFeSpace(i), getFeSpace(i));
opDivGradU->addTerm(new Simple_SOT(mu));
diff --git a/extensions/base_problems/LinearElasticity.h b/extensions/base_problems/LinearElasticity.h
index 7db74f0f..b387507c 100644
--- a/extensions/base_problems/LinearElasticity.h
+++ b/extensions/base_problems/LinearElasticity.h
@@ -45,6 +45,8 @@ protected: // variables
double mu;
double lambda;
double rho;
+
+ int nonLinTerm;
WorldVector<double> force;
DOFVector<WorldVector<double> >* displacement;
diff --git a/extensions/base_problems/NavierStokes_TaylorHood.cc b/extensions/base_problems/NavierStokes_TaylorHood.cc
index b8d676bf..fd5db021 100644
--- a/extensions/base_problems/NavierStokes_TaylorHood.cc
+++ b/extensions/base_problems/NavierStokes_TaylorHood.cc
@@ -172,11 +172,11 @@ void NavierStokes_TaylorHood::fillOperators()
}
}
- for (int j = 0; j < dow; ++j) {
- Operator *opNull = new Operator(getFeSpace(i), getFeSpace(j));
- opNull->addTerm(new Simple_ZOT(0.0));
- prob->addMatrixOperator(*opNull, i, j);
- }
+// for (int j = 0; j < dow; ++j) {
+// Operator *opNull = new Operator(getFeSpace(i), getFeSpace(j));
+// opNull->addTerm(new Simple_ZOT(0.0));
+// prob->addMatrixOperator(*opNull, i, j);
+// }
/// Diffusion-Operator (including Stress-Tensor for space-dependent viscosity
addLaplaceTerm(i);
--
GitLab