Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
iwr
amdis
Commits
4baf4abe
Commit
4baf4abe
authored
Aug 15, 2012
by
Praetorius, Simon
Browse files
base problems
parent
13a42661
Changes
11
Hide whitespace changes
Inline
Side-by-side
extensions/base_problems/CahnHilliardNavierStokes.cc
View file @
4baf4abe
...
...
@@ -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
));
opTime
Old
->
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
);
}
extensions/base_problems/CahnHilliardNavierStokes.h
View file @
4baf4abe
...
...
@@ -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
);
...
...
extensions/base_problems/CahnHilliardNavierStokes_RB.cc
0 → 100644
View file @
4baf4abe
#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
);
}
extensions/base_problems/CahnHilliardNavierStokes_RB.h
0 → 100644
View file @
4baf4abe
/** \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
;
}