From 81d60deecbc63a5a65de1ad7a9b3f36455ed3642 Mon Sep 17 00:00:00 2001
From: Lars von Wolff <lars.von-wolff@mathematik.uni-stuttgart.de>
Date: Fri, 26 Aug 2022 12:24:42 +0200
Subject: [PATCH] Final State of Bachelorthesis
---
dune/phasefield/base_parameters.hh | 2 +
dune/phasefield/porenetwork/pn_initial.hh | 98 ++++++++++++++++++++---
src/pntest/pntest.cc | 92 +++++++++++++++------
3 files changed, 157 insertions(+), 35 deletions(-)
diff --git a/dune/phasefield/base_parameters.hh b/dune/phasefield/base_parameters.hh
index 121e32c..fab0667 100644
--- a/dune/phasefield/base_parameters.hh
+++ b/dune/phasefield/base_parameters.hh
@@ -290,6 +290,8 @@ class Params_Initial
public:
const Dune::ParameterTree ptree;
int shape;
+ float a_test;
+ float b_test;
Params_Initial(const Dune::ParameterTree& 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 ee50fb2..6e7f3fd 100644
--- a/dune/phasefield/porenetwork/pn_initial.hh
+++ b/dune/phasefield/porenetwork/pn_initial.hh
@@ -1,9 +1,10 @@
#pragma once
-
+#include <math.h>
#define INITIAL_SHAPE_SQUARE 0
#define INITIAL_SHAPE_CIRCLE 1
#define INITIAL_SHAPE_ELLIPSE 2
#define INITIAL_SHAPE_TWOCIRCLES 3
+#define INITIAL_SHAPE_RECTANGLE 4
template<typename GV, typename RF, typename Params>
class PhiInitial_Quartersquare :
@@ -14,6 +15,8 @@ class PhiInitial_Quartersquare :
RF eps;
RF r;
int myshape;
+ float a;
+ float b;
public:
typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1> Traits;
@@ -26,6 +29,8 @@ public:
r = params.initial.ptree.get("Radius",0.2);
std::cout << "EPS" << eps;
myshape = params.initial.shape;
+ a = params.initial.a_test;
+ b = params.initial.b_test;
}
inline void evaluateGlobal(const typename Traits::DomainType & x, typename Traits::RangeType & y) const
@@ -45,25 +50,94 @@ public:
}
else if(myshape == INITIAL_SHAPE_ELLIPSE)
{
- RF rad = sqrt((x[0]-1)/1.2*(x[0]-1)/1.2 + (x[1]-1)/0.8*(x[1]-1)/0.8);
- y = param.dw.eval_shape(((1-r)-rad)*1/eps);
+ float xnew1 = (x[1]-1)/a+1;
+ RF rad = sqrt((x[0]-1)*(x[0]-1) + (xnew1-1)*(xnew1-1));
+ RF y1 = param.dw.eval_shape(((1-r)-rad)*1/eps); //
+ y = y1;
+ }
+ else if(myshape == INITIAL_SHAPE_RECTANGLE) {
+ RF rf1 = param.dw.eval_shape((x[0]-r)*1/eps);
+ RF rf2 = param.dw.eval_shape((x[1]-(1-0.8*a))*1/eps);
+ y=rf1*rf2;
}
else if(myshape == INITIAL_SHAPE_TWOCIRCLES)
{
- float radius1 = 1.75; //radius circle 1
- float radius2 = 1.75; //radius circle 2
- float x01 = 0.9; //xo circle 1
- float x11 = 1.5; //x1 circle 1
- float x02 = 1.5; //x0 circle 2
- float x12 = 1.2; //x1 circle 2
+ float b_;
+ //convert parameter b to distance of circle centers to corner b_
+ if(b > 10000/10001.0f) {
+ b_ = 10000;
+ }
+ else {
+ b_ = 1/float(1.0f - b) - 1.0f;
+ }
+ int method = 3;
+ if(method == 0) {
+ // Idee: Der Radius eines Kreises wird verkleinert, um Stauchung zu erreichen
+ // Problem: Dadurch wird auch in die andere Richtung gestaucht, besonders auffällig
+ // bei b=0
+ float radius1 = 1.25*b_ + 1; //radius circle 1, Faktor 1.25, weil ja mit 0.8 skaliert wird
+ float radius2 = 1.25*b_ + 1*a; //radius circle 2, nur der Teil, der tatsächlich im "Bild" ist, wird mit a skaliert
+ float x01 = 1; //xo circle 1
+ float x11 = b_+1; //x1 circle 1
+ float x02 = b_+1; //x0 circle 2
+ float x12 = 1; //x1 circle 2
RF rad = sqrt((x[0]-x01)*(x[0]-x01) + (x[1]-x11)*(x[1]-x11));
RF y1 = param.dw.eval_shape(((1-r)*radius1-rad)*1/eps);
RF rad2 = sqrt((x[0]-x02)*(x[0]-x02) + (x[1]-x12)*(x[1]-x12));
RF y2 = param.dw.eval_shape(((1-r)*radius2-rad2)*1/eps);
y = y1 * y2;
- //RF rad1 = sqrt((x[0]-2)*(x[0]-2) + (x[1]-1)*(x[1]-1));
- //RF rad2 = sqrt((x[0]-1)*(x[0]-1) + (x[1]-2)*(x[1]-2));
- //y = param.dw.eval_shape(((1-r)-rad2)*1/eps) * param.dw.eval_shape(((1-r)-rad1)*1/eps);
+ }
+ else if(method == 1){
+ //Idee: Einen Kreis als Ellipse formulieren, um Stauchung in zweite Richtung
+ //zu verhindern
+ //Problem: sehr breite Übergangsregion, besonders auffällig für b = 1
+ //float radius1 = 1.25*b_ + 1; //radius circle 1
+ float radius2 = 1.25*b_ + 1; //radius circle 2
+ float x01 = 1; //xo circle 1
+ float x11 = b_+1; //x1 circle 1
+ float x02 = b_+1; //x0 circle 2
+ float x12 = 1; //x1 circle 2
+ //Ellipse mit Radien 1.25b+a und 1.25b+1
+ RF rad = sqrt((x[0]-x01)/(1.25*b_+1)*(x[0]-x01)/(1.25*b_+1) + (x[1]-x11)/(1.25*b_ + a)*(x[1]-x11)/(1.25*b_ + a));
+ RF y1 = param.dw.eval_shape(((1-r)-rad)*1/eps); //
+ RF rad2 = sqrt((x[0]-x02)*(x[0]-x02) + (x[1]-x12)*(x[1]-x12));
+ RF y2 = param.dw.eval_shape(((1-r)*radius2-rad2)*1/eps);
+ y = y1 * y2;
+ }
+ else if ( method == 2) {
+ //rumbasteln/sonstige Tests
+ float radius = 1.25 * b_ + 1;
+ float x01 = 1; //xo circle 1
+ float x11 = b_+1; //x1 circle 1
+ float x02 = b_+1; //x0 circle 2
+ float x12 = 1; //x1 circle 2
+ RF rad = sqrt((x[0]-x01)*(x[0]-x01) *(a*a)+ (x[1]-x11)*(x[1]-x11));
+ RF y1 = param.dw.eval_shape(((1-r)*radius-rad)*1/eps);
+ RF rad2 = sqrt((x[0]-x02)*(x[0]-x02) + (x[1]-x12)*(x[1]-x12));
+ RF y2 = param.dw.eval_shape(((1-r)*radius-rad2)*1/eps);
+ y = y1 * y2;
+ }
+ else if(method == 3){
+ //Idee: Einen Kreis als Ellipse formulieren, um Stauchung in zweite Richtung
+ //zu verhindern
+ //Problem: sehr breite Übergangsregion, besonders auffällig für b = 1
+ //float radius1 = 1.25*b_ + 1; //radius circle 1
+ //float xnew0 = x[0];
+ float xnew1 = (x[1]-1)/a+1;
+
+ float radius2 = 1.25*b_ + 1; //radius circle 2
+ float x01 = 1; //xo circle 1
+ float x11 = b_+1; //x1 circle 1
+ float x02 = b_+1; //x0 circle 2
+ float x12 = 1; //x1 circle 2
+ //Ellipse mit Radien 1.25b+a und 1.25b+1
+ RF rad = sqrt((x[0]-x01)*(x[0]-x01) + (xnew1-x11)*(xnew1-x11));
+ RF y1 = param.dw.eval_shape(((1-r)*radius2-rad)*1/eps); //
+ RF rad2 = sqrt((x[0]-x02)*(x[0]-x02) + (xnew1-x12)*(xnew1-x12));
+ RF y2 = param.dw.eval_shape(((1-r)*radius2-rad2)*1/eps);
+ y = y1 * y2;
+ }
+
}
else
{
diff --git a/src/pntest/pntest.cc b/src/pntest/pntest.cc
index 5ae546e..6ee5c7f 100644
--- a/src/pntest/pntest.cc
+++ b/src/pntest/pntest.cc
@@ -53,41 +53,87 @@ int main(int argc, char** argv)
//Set input/output filenamens
std::string GridFilename = ptree.get("domain.filename","square.msh");
std::string OutputFilename = ptree.get("output.filename","output");
-
+ bool runAll = false;
//Set initial Shape
//param.initial.shape = INITIAL_SHAPE_SQUARE;
//param.initial.shape = INITIAL_SHAPE_CIRCLE;
//param.initial.shape = INITIAL_SHAPE_ELLIPSE;
param.initial.shape = INITIAL_SHAPE_TWOCIRCLES;
-
+ param.initial.a_test = 1;
+ param.initial.b_test = 1;
+ float a_values[20] = {0.05, 0.1,0.15, 0.2, 0.25, 0.3, 0.35, 0.4,0.45, 0.5,0.55, 0.6,0.65, 0.7,0.75, 0.8, 0.85,0.9, 0.95,1.0};
+ float b_values[21] = {0, 0.05, 0.1,0.15, 0.2, 0.25, 0.3, 0.35, 0.4,0.45, 0.5,0.55, 0.6,0.65, 0.7,0.75, 0.8, 0.85,0.9, 0.95,1.0};
//Create Simulation
- Pn_1PFlow<Parameters> pn_1pFlow(GridFilename,OutputFilename,param, 10e-6); //Last Parameter is initial pore Radius in m
+ //for (int i = 0; i < 11; ++i) {
+ //param.initial.b_test = b_values[i];
+ //Pn_1PFlow<Parameters> pn_1pFlow(GridFilename,OutputFilename,param, 10e-6); //Last Parameter is initial pore Radius in m
+
+ //Example of usecase for Pn_1PFlow
+ //std::cout << "b= " << param.initial.b_test << std::endl;
+ //std::cout << "Kf = " << pn_1pFlow.getKf() << " phiSolid = " << pn_1pFlow.getPhiSolid() << std::endl;
+ //std::cout << "Pore Area = " << pn_1pFlow.getPoreArea() << " m^2" << std::endl;
+ //std::cout << "Pore Area Fraction (compared to inital values) = " << pn_1pFlow.getPoreAreaFraction() << std::endl;
+ //pn_1pFlow.writeVTK(0.0);
+ //opening stream to write into csv file
+ //std::string filename("squ.csv");
+ //std::ofstream file_out;
+ //file_out.open(filename, std::ios_base::app);
- //Example of usecase for Pn_1PFlow
- std::cout << "Kf = " << pn_1pFlow.getKf() << " phiSolid = " << pn_1pFlow.getPhiSolid() << std::endl;
- std::cout << "Pore Area = " << pn_1pFlow.getPoreArea() << " m^2" << std::endl;
- std::cout << "Pore Area Fraction (compared to inital values) = " << pn_1pFlow.getPoreAreaFraction() << std::endl;
- pn_1pFlow.writeVTK(0.0);
- std::string filename("test.csv");
+ //loop until phi solid is > 0.85
+ //int j = 0;
+ //while(1 - pn_1pFlow.getPhiSolid() > 0.55) {
+ //j++;
+ //pn_1pFlow.grow(2e-6); //Growth is measured in m
+ //pn_1pFlow.growArea((10e-6)*(10e-6)*0.1); // Growth of Area in m^2
+ //pn_1pFlow.writeVTK(1.0);
+ //pn_1pFlow.writeVTK(j);
+ //std::cout << "Kf = " << pn_1pFlow.getKf() << " phiSolid = " << pn_1pFlow.getPhiSolid() << std::endl;
+ //std::cout << "Total Flux = - Kf/viscosity * gradient p" << std::endl;
+ //std::cout << "Kf has unit m^4" << std::endl;
+
+ //std::cout << "Pore Area = " << pn_1pFlow.getPoreArea() << " m^2" << std::endl;
+ //std::cout << "Pore Area Fraction (compared to inital values) = " << pn_1pFlow.getPoreAreaFraction() << std::endl;
+ //file_out << pn_1pFlow.getKf() << "," << 1 - pn_1pFlow.getPhiSolid() << "b=" << param.initial.b_test << std::endl;
+ //}
+
+ if(runAll) {
+ //opening stream to write into csv file
+ std::string filename("data_rectangles.csv");
std::ofstream file_out;
file_out.open(filename, std::ios_base::app);
- int i = 0;
- while(1 - pn_1pFlow.getPhiSolid() > 0.15) {
- i++;
- //pn_1pFlow.grow(2e-6); //Growth is measured in m
- pn_1pFlow.growArea((10e-6)*(10e-6)*0.1); // Growth of Area in m^2
- //pn_1pFlow.writeVTK(1.0);
- pn_1pFlow.writeVTK(i);
- std::cout << "Kf = " << pn_1pFlow.getKf() << " phiSolid = " << pn_1pFlow.getPhiSolid() << std::endl;
- std::cout << "Total Flux = - Kf/viscosity * gradient p" << std::endl;
- std::cout << "Kf has unit m^4" << std::endl;
+ int step;
+ for(auto const& a: a_values) {
+ for(auto const& b: b_values) {
+ param.initial.b_test = b;
+ param.initial.a_test = a;
+ std::cout << "a= " << a << ", b=" << b << std::endl;
+ Pn_1PFlow<Parameters> pn_1pFlow(GridFilename,OutputFilename,param, 10e-6); //Last Parameter is initial pore Radius in m
+ float kf_init = pn_1pFlow.getKf();
+ float phiSolid_init = pn_1pFlow.getPhiSolid();
+ float poreArea_init = pn_1pFlow.getPoreArea();
+ step = 0;
+ file_out << a << ","<< b <<","<< pn_1pFlow.getKf() << "," << pn_1pFlow.getPhiSolid() << "," << pn_1pFlow.getPoreArea() << "," << kf_init << ","<< phiSolid_init << ","<< poreArea_init << "," << pn_1pFlow.getPoreAreaFraction() << ","<< step << std::endl;
+ while(1 - pn_1pFlow.getPhiSolid() > 0.01) {
+ pn_1pFlow.grow(1e-7); //Growth is measured in m
+ step++;
+ file_out << a << ","<< b <<","<< pn_1pFlow.getKf() << "," << pn_1pFlow.getPhiSolid() << "," << pn_1pFlow.getPoreArea() << "," << kf_init << ","<< phiSolid_init << ","<< poreArea_init << "," << pn_1pFlow.getPoreAreaFraction() << ","<< step << std::endl;
- std::cout << "Pore Area = " << pn_1pFlow.getPoreArea() << " m^2" << std::endl;
- std::cout << "Pore Area Fraction (compared to inital values) = " << pn_1pFlow.getPoreAreaFraction() << std::endl;
- file_out << pn_1pFlow.getKf() << "," << 1 - pn_1pFlow.getPhiSolid() << std::endl;
+ }
+ }
+ }
+ }//if
+ else {
+ param.initial.shape = INITIAL_SHAPE_SQUARE;
+ Pn_1PFlow<Parameters> pn_1pFlow(GridFilename,OutputFilename,param, 10e-6); //Last Parameter is initial pore Radius in m
+ pn_1pFlow.writeVTK(0);
+ int i = 1;
+ while(1 - pn_1pFlow.getPhiSolid() > 0.01) {
+ pn_1pFlow.grow(1e-7); //Growth is measured in m
+ pn_1pFlow.writeVTK(i);
+ i++;
+ }
}
-
}
catch (Dune::Exception &e){
std::cerr << "Dune reported error: " << e << std::endl;
--
GitLab