Skip to content
Snippets Groups Projects
Commit ece2ecfa authored by Lars von Wolff's avatar Lars von Wolff
Browse files

Started Work of Bachelor Thesis

parent fc46149a
No related branches found
No related tags found
No related merge requests found
...@@ -2,6 +2,8 @@ ...@@ -2,6 +2,8 @@
#define INITIAL_SHAPE_SQUARE 0 #define INITIAL_SHAPE_SQUARE 0
#define INITIAL_SHAPE_CIRCLE 1 #define INITIAL_SHAPE_CIRCLE 1
#define INITIAL_SHAPE_ELLIPSE 2
#define INITIAL_SHAPE_TWOCIRCLES 3
template<typename GV, typename RF, typename Params> template<typename GV, typename RF, typename Params>
class PhiInitial_Quartersquare : class PhiInitial_Quartersquare :
...@@ -41,6 +43,28 @@ public: ...@@ -41,6 +43,28 @@ public:
RF rad = sqrt((x[0]-1)*(x[0]-1) + (x[1]-1)*(x[1]-1)); 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); y = param.dw.eval_shape(((1-r)-rad)*1/eps);
} }
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);
}
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
RF rad = sqrt((x[0]-x01)*(x[0]-x01) + (x[1]-x11)*(x[1]-x11));
RF y1 = param.dw.eval_shape(((1-r)-(rad)/radius1)*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)-(rad2)/radius2)*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 else
{ {
std::cout << "WARNING: pn_initial.hh cannot find shape" << myshape << std::endl; std::cout << "WARNING: pn_initial.hh cannot find shape" << myshape << std::endl;
......
#add_executable(pntest pntest.cc) add_executable(pntest pntest.cc)
add_executable(pn2ptest pn2ptest.cc) #add_executable(pn2ptest pn2ptest.cc)
dune_symlink_to_source_files(FILES grids pn1p.ini pn2p.ini) dune_symlink_to_source_files(FILES grids pn1p.ini pn2p.ini)
...@@ -3,6 +3,7 @@ ...@@ -3,6 +3,7 @@
#endif #endif
#include<iostream> #include<iostream>
#include<fstream>
// dune includes // dune includes
#include<dune/common/parallel/mpihelper.hh> #include<dune/common/parallel/mpihelper.hh>
#include<dune/common/parametertreeparser.hh> #include<dune/common/parametertreeparser.hh>
...@@ -55,27 +56,37 @@ int main(int argc, char** argv) ...@@ -55,27 +56,37 @@ int main(int argc, char** argv)
//Set initial Shape //Set initial Shape
//param.initial.shape = INITIAL_SHAPE_SQUARE; //param.initial.shape = INITIAL_SHAPE_SQUARE;
param.initial.shape = INITIAL_SHAPE_CIRCLE; //param.initial.shape = INITIAL_SHAPE_CIRCLE;
//param.initial.shape = INITIAL_SHAPE_ELLIPSE;
param.initial.shape = INITIAL_SHAPE_TWOCIRCLES;
//Create Simulation //Create Simulation
Pn_1PFlow<Parameters> pn_1pFlow(GridFilename,OutputFilename,param, 10e-6); //Last Parameter is initial pore Radius in m Pn_1PFlow<Parameters> pn_1pFlow(GridFilename,OutputFilename,param, 10e-6); //Last Parameter is initial pore Radius in m
//Example of usecase for Pn_1PFlow //Example of usecase for Pn_1PFlow
std::cout << "Kf = " << pn_1pFlow.getKf() << " phiSolid = " << pn_1pFlow.getPhiSolid() << 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 = " << pn_1pFlow.getPoreArea() << " m^2" << std::endl;
std::cout << "Pore Area Fraction (compared to inital values) = " << pn_1pFlow.getPoreAreaFraction() << std::endl; std::cout << "Pore Area Fraction (compared to inital values) = " << pn_1pFlow.getPoreAreaFraction() << std::endl;
pn_1pFlow.writeVTK(0.0); pn_1pFlow.writeVTK(0.0);
std::string filename("test.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.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.growArea((10e-6)*(10e-6)*0.1); // Growth of Area in m^2
pn_1pFlow.writeVTK(1.0); //pn_1pFlow.writeVTK(1.0);
pn_1pFlow.writeVTK(i);
std::cout << "Kf = " << pn_1pFlow.getKf() << " phiSolid = " << pn_1pFlow.getPhiSolid() << std::endl; std::cout << "Kf = " << pn_1pFlow.getKf() << " phiSolid = " << pn_1pFlow.getPhiSolid() << std::endl;
std::cout << "Total Flux = - Kf/viscosity * gradient p" << std::endl; std::cout << "Total Flux = - Kf/viscosity * gradient p" << std::endl;
std::cout << "Kf has unit m^4" << 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 = " << pn_1pFlow.getPoreArea() << " m^2" << std::endl;
std::cout << "Pore Area Fraction (compared to inital values) = " << pn_1pFlow.getPoreAreaFraction() << 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;
}
} }
catch (Dune::Exception &e){ catch (Dune::Exception &e){
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment