From 9f213e065e9c4668492452e7b409fb904532389b Mon Sep 17 00:00:00 2001
From: Lars von Wolff <lars.von-wolff@ians.uni-stuttgart.de>
Date: Wed, 2 Jun 2021 15:38:26 +0200
Subject: [PATCH] Added the option to choose between shapes for the initial
values. Currently implemented: Square and Circle
---
dune/phasefield/base_parameters.hh | 3 ++-
dune/phasefield/porenetwork/pn_initial.hh | 28 ++++++++++++++++++-----
src/pntest/pn1p.ini | 7 +++---
src/pntest/pntest.cc | 7 +++++-
4 files changed, 34 insertions(+), 11 deletions(-)
diff --git a/dune/phasefield/base_parameters.hh b/dune/phasefield/base_parameters.hh
index 950a9da..33d0dbe 100644
--- a/dune/phasefield/base_parameters.hh
+++ b/dune/phasefield/base_parameters.hh
@@ -256,8 +256,9 @@ class Params_Initial
{
public:
const Dune::ParameterTree ptree;
+ int shape;
Params_Initial(const Dune::ParameterTree& tree) :
- ptree(tree)
+ ptree(tree), shape(ptree.get("shape",0))
{
}
};
diff --git a/dune/phasefield/porenetwork/pn_initial.hh b/dune/phasefield/porenetwork/pn_initial.hh
index df0c621..2f8c77d 100644
--- a/dune/phasefield/porenetwork/pn_initial.hh
+++ b/dune/phasefield/porenetwork/pn_initial.hh
@@ -1,5 +1,8 @@
#pragma once
+#define INITIAL_SHAPE_SQUARE 0
+#define INITIAL_SHAPE_CIRCLE 1
+
template<typename GV, typename RF, typename Params>
class PhiInitial_Quartersquare :
public Dune::PDELab::AnalyticGridFunctionBase<
@@ -8,6 +11,7 @@ class PhiInitial_Quartersquare :
const Params& param;
RF eps;
RF r;
+ int myshape;
public:
typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1> Traits;
@@ -19,15 +23,27 @@ public:
eps = params.eps;
r = params.initial.ptree.get("Radius",0.2);
std::cout << "EPS" << eps;
+ myshape = params.initial.shape;
}
inline void evaluateGlobal(const typename Traits::DomainType & x, typename Traits::RangeType & y) const
{
- RF rf1 = param.dw.eval_shape((x[0]-r)*1/eps);
- RF rf2 = param.dw.eval_shape((x[1]-r)*1/eps);
- //RF rf1 = 0.5*(1+tanh(1/eps/sqrt(2)*(x[0]-r)));
- //RF rf2 = 0.5*(1+tanh(1/eps/sqrt(2)*(x[1]-r)));
-
- y=rf1*rf2;
+ if(myshape == INITIAL_SHAPE_SQUARE)
+ {
+ RF rf1 = param.dw.eval_shape((x[0]-r)*1/eps);
+ RF rf2 = param.dw.eval_shape((x[1]-r)*1/eps);
+ //RF rf1 = 0.5*(1+tanh(1/eps/sqrt(2)*(x[0]-r)));
+ //RF rf2 = 0.5*(1+tanh(1/eps/sqrt(2)*(x[1]-r)));
+ y=rf1*rf2;
+ }
+ else if(myshape == INITIAL_SHAPE_CIRCLE)
+ {
+ RF rad = sqrt((x[0]-1)*(x[0]-1) + (x[1]-1)*(x[1]-1));
+ y = param.dw.eval_shape(((1-r)-rad)*1/eps);
+ }
+ else
+ {
+ std::cout << "WARNING: pn_initial.hh cannot find shape" << myshape << std::endl;
+ }
}
};
diff --git a/src/pntest/pn1p.ini b/src/pntest/pn1p.ini
index e695424..d78203a 100644
--- a/src/pntest/pn1p.ini
+++ b/src/pntest/pn1p.ini
@@ -29,10 +29,11 @@ Sigma3=1
dtMax = 0.002
[Initial]
+shape = 0
# FlowVyy= 10# pConst/(Domain_width * mu)
-pConst=0.1
-uConst=0.5
-uInflow=0.3
+#pConst=0.1
+#uConst=0.5
+#uInflow=0.3
[Refinement]
etaRefine = 0.2
diff --git a/src/pntest/pntest.cc b/src/pntest/pntest.cc
index 1d66b2e..e078868 100644
--- a/src/pntest/pntest.cc
+++ b/src/pntest/pntest.cc
@@ -48,10 +48,15 @@ int main(int argc, char** argv)
> Parameters;
Parameters param(ptree);
+
//Set input/output filenamens
std::string GridFilename = ptree.get("domain.filename","square.msh");
std::string OutputFilename = ptree.get("output.filename","output");
+ //Set initial Shape
+ //param.initial.shape = INITIAL_SHAPE_SQUARE;
+ //param.initial.shape = INITIAL_SHAPE_CIRCLE;
+
//Create Simulation
Pn_1PFlow<Parameters> pn_1pFlow(GridFilename,OutputFilename,param, 10e-6); //Last Parameter is initial pore Radius in m
@@ -59,7 +64,7 @@ int main(int argc, char** argv)
//Example of usecase for Pn_1PFlow
std::cout << "Kf = " << pn_1pFlow.getKf() << " phiSolid = " << pn_1pFlow.getPhiSolid() << std::endl;
pn_1pFlow.writeVTK(0.0);
- pn_1pFlow.grow(2e-6); //Growth is measured in m
+ pn_1pFlow.grow(0.2e-6); //Growth is measured in m
pn_1pFlow.writeVTK(1.0);
std::cout << "Kf = " << pn_1pFlow.getKf() << " phiSolid = " << pn_1pFlow.getPhiSolid() << std::endl;
std::cout << "Total Flux = - Kf/viscosity * gradient p" << std::endl;
--
GitLab