diff --git a/dune/phasefield/base_parameters.hh b/dune/phasefield/base_parameters.hh
index 47e93008ed5af59af34f174cad9380c8d72b2ad0..121e32c17dd1638c5138f11a7348f355860a655b 100644
--- a/dune/phasefield/base_parameters.hh
+++ b/dune/phasefield/base_parameters.hh
@@ -450,17 +450,22 @@ public:
Number pf1Forcing;
Number pf2Forcing;
+ Number PCS_Search_VelThreshold;
+ Number PCS_Search_MaxCHNSIter;
+
Params_pn_2pflow(const Dune::ParameterTree& ptree) :
Params_base<Number, Mobility>(ptree),
delta(ptree.get("Phasefield.delta",this->eps)),
- mu_ratio(ptree.get("Pnysics.mu_ratio",1)),
+ mu_ratio(ptree.get("Physics.mu_ratio",(Number)1)),
tw(ptree.sub("Phasefield")),
htw(ptree.sub("Phasefield")),
reaction(ptree.sub("Physics")),
vDis(ptree.sub("Physics")),
SurfaceTensionScaling(ptree.get("Physics.SurfaceTensionScaling",(Number)1)),
pf1Forcing(0),
- pf2Forcing(0)
+ pf2Forcing(0),
+ PCS_Search_VelThreshold(ptree.get("PCS_Search.VelThreshold",(Number)1e-3)),
+ PCS_Search_MaxCHNSIter(ptree.get("PCS_Search.MaxCHNSIter",(Number)100))
{
}
};
diff --git a/dune/phasefield/porenetwork/pn_2pflow.hh b/dune/phasefield/porenetwork/pn_2pflow.hh
index a2bcca9a144132849ac9c5dc611c53707e19b8be..9af50e0faef5a1a362423201548197589557f7d7 100644
--- a/dune/phasefield/porenetwork/pn_2pflow.hh
+++ b/dune/phasefield/porenetwork/pn_2pflow.hh
@@ -161,7 +161,9 @@ private:
std::unique_ptr<NSX_T_Newton> nsxtNewton;
std::unique_ptr<Dune::VTKSequenceWriter<GV>> vtkwriter;
- //------------------------------------------------------------------------------------
+//------------------------------------------------------------------------------------
+// Operators and Solving
+//------------------------------------------------------------------------------------
void initializeGridOperators()
@@ -238,16 +240,24 @@ private:
remesh();
}
+
+//------------------------------------------------------------------------------------
+//Integration and Helper Functions
+//------------------------------------------------------------------------------------
+
+
void updateKf()
{
param.pf1Forcing = 1;
param.pf2Forcing = 0;
applyNsTransmissibility();
updateKf_integration(true);
+ writeVTK(0);
param.pf1Forcing = 0;
param.pf2Forcing = 1;
applyNsTransmissibility();
updateKf_integration(false);
+ writeVTK(1);
}
void updateKf_integration(bool phase1)
@@ -318,13 +328,16 @@ private:
capillaryPressure = pressureFluid2-pressureFluid1;
}
-/* RF calculateGrowthSpeed()
+ RF calculateGrowthSpeed()
{
auto productFunction = Dune::PDELab::makeGridFunctionFromCallable (gv, [&](const auto& element, const auto& xlocal){
Dune::FieldVector<double,1> phi1;
dgf->phiDgf.evaluate(element, xlocal, phi1);
+ Dune::FieldVector<double,1> phi2;
+ dgf->phi2Dgf.evaluate(element, xlocal, phi2);
+ RF phi3 = 1 - phi1[0] - phi2[0];
- RF Rspeed = phi1[0] < 0.95 ? param.dw.eval_q(phi1[0]) : 0.0f;
+ RF Rspeed = phi1[0] < 0.95 ? param.tw.eval_q(phi1[0],phi2[0],phi3) : 0.0f;
RF R1 = -Rspeed/param.eps*param.reaction.eval(1);
return Dune::FieldVector<double,1>(-R1);
});
@@ -333,7 +346,36 @@ private:
return integral[0];
}
-*/
+ void findStationaryState()
+ {
+
+ writeVTK(0);
+ Dune::FieldVector<double,1> integral;
+ param.reaction.reactionRate = 0;
+ int iter = 0;
+ RF dtmax = param.pTree.get("Time.dtMax",0.01);
+ param.time.dt = dtmax;
+
+ do
+ {
+ iter++;
+ applyChPressureSaturation(1);
+ applyNsPressureSaturation(1);
+ auto productFunction = Dune::PDELab::makeGridFunctionFromCallable (gv, [&](const auto& element, const auto& xlocal){
+ Dune::FieldVector<double,2> v;
+ dgf->vDgf.evaluate(element, xlocal, v);
+ Dune::FieldVector<double,1> phi1;
+ dgf->phiDgf.evaluate(element, xlocal, phi1);
+ Dune::FieldVector<double,1> phi2;
+ dgf->phi2Dgf.evaluate(element, xlocal, phi2);
+ return Dune::FieldVector<double,1>((phi1+phi2)*std::sqrt(v[0]*v[0]+v[1]*v[1]));
+ });
+ Dune::PDELab::integrateGridFunction(productFunction,integral,2);
+ //writeVTK(iter+1);
+ }
+ while ((integral[0]>param.PCS_Search_VelThreshold)
+ && (iter < param.PCS_Search_MaxCHNSIter));
+ }
//------------------------------------------------------------------------------------
@@ -403,146 +445,234 @@ public:
initialPhiSolid = phiSolid;
updateKf();
*/
+ //findStationaryState();
+ updateKf();
+ std::cout << kf_11 << ", "
+ << kf_12 << ", "
+ << kf_21 << ", "
+ << kf_22 << std::endl;
+ }
- RF dtmax = param.pTree.get("Time.dtMax",0.01);
- param.time.dt = dtmax;
- int InitialGrowthTime = param.pTree.get("test.InitialGrowthTime",0);
- int FillTime = param.pTree.get("test.FillTime",100);
- int EmptyTime = param.pTree.get("test.EmptyTime",200);
+//------------------------------------------------------------------------------------
- param.reaction.reactionRate = 0;
- for(int j=0; j<100; j++)
+
+ void writeVTK(RF time)
+ {
+ vtkwriter->write(time,Dune::VTK::ascii);
+ }
+
+ void writeCSV(std::string filename, RF time = 0, int iteration = 0)
+ {
+ try
{
- applyChPressureSaturation(1);
- applyNsPressureSaturation(1);
+ std::ofstream myfile;
+ myfile.open(filename,std::ofstream::app);
+ myfile << time << ", "
+ << iteration << ", "
+ << phiFluid1 << ", "
+ << phiFluid2 << ", "
+ << phiSolid << ", "
+ << saturation << ", "
+ << pressureFluid1 << ", "
+ << pressureFluid2 << ", "
+ << capillaryPressure << ", "
+ << kf_11 << ", "
+ << kf_12 << ", "
+ << kf_21 << ", "
+ << kf_22 << std::endl;
+ myfile.close();
}
- for(int j=0; j<InitialGrowthTime; j++)
+ catch(...)
{
- param.reaction.reactionRate = 1;
- applyChReaction(1);
+ std::cout << "Error in write CSV File" << std::endl;
}
- param.reaction.reactionRate = 0;
- for(int j=0; j<100; j++)
+ }
+
+ void writeCSVHeader(std::string filename)
+ {
+ try
{
- applyChPressureSaturation(1);
- applyNsPressureSaturation(1);
+ std::ofstream myfile;
+ myfile.open(filename,std::ofstream::app);
+ myfile << "Time" << ", "
+ << "Iteration" << ", "
+ << "Area Fluid1" << ", "
+ << "Area Fluid2" << ", "
+ << "Area Solid" << ", "
+ << "Saturation" << ", "
+ << "Pressure Fluid1" << ", "
+ << "Pressure Fluid2" << ", "
+ << "Capillary Pressure" << ", "
+ << "kf_11" << ", "
+ << "kf_12" << ", "
+ << "kf_21" << ", "
+ << "kf_22" << std::endl;
+ myfile.close();
}
+ catch(...)
+ {
+ std::cout << "Error in write CSV File" << std::endl;
+ }
+ }
+
- for(int i=0; i<FillTime+EmptyTime; i++)
+//------------------------------------------------------------------------------------
+
+ void calculateFullPCS_Curve()
+ {
+ RF dtmax = param.pTree.get("Time.dtMax",0.01);
+ param.time.dt = dtmax;
+
+ int InitialGrowthIter = param.pTree.get("FullPCS_Curve.InitialGrowthIter",0);
+ int FillIter = param.pTree.get("FullPCS_Curve.FillIter",100);
+ int EmptyIter = param.pTree.get("FullPCS_Curve.EmptyIter",200);
+
+ findStationaryState();
+ for(int j=0; j<InitialGrowthIter; j++)
{
- if(i<FillTime)
+ param.reaction.reactionRate = 1;
+ applyChReaction(1);
+ }
+ findStationaryState();
+
+ writeCSVHeader("pressuresaturation.csv");
+ for(int i=0; i<FillIter+EmptyIter; i++)
+ {
+ if(i<FillIter)
{
param.reaction.reactionRate = 1;
+ applyChPressureSaturation(1);
}
else
{
param.reaction.reactionRate = -1;
- }
- applyChPressureSaturation(1);
- param.reaction.reactionRate = 0;
- for(int j=0; j<20; j++)
- {
applyChPressureSaturation(1);
- applyNsPressureSaturation(1);
- //writeVTK(i+((RF)j+1)/21);
}
+
+ findStationaryState();
+
writeVTK(i+1);
updateSaturation();
updateCapillaryPressure();
updateKf();
- try
- {
- std::ofstream myfile;
- myfile.open("pressuresaturation.csv",std::ofstream::app);
- myfile << i << ", " << phiFluid1 << ", " << phiFluid2 << ", " << phiSolid << ", " << saturation << ", "
- << pressureFluid1 << ", " << pressureFluid2 << ", " << capillaryPressure << ", "
- << kf_11 << ", " << kf_12 << ", " << kf_21 << ", " << kf_22 << ", " << std::endl;
- myfile.close();
- }
- catch(...)
- {
- std::cout << "Error in integration or CSV Files" << std::endl;
- }
-
-
- /*if(i < 10)
- {
- applyChPressureSaturation(1);
- }
- else if(i<40)
- {
- param.reaction.reactionRate = 0;
- applyChPressureSaturation(1);
- applyNsPressureSaturation(1);
- }
- else
- {
- param.reaction.reactionRate = 1;
- applyChReaction(1);
- }*/
+ writeCSV("pressuresaturation.csv", 0, i);
}
}
- void writeVTK(RF time)
- {
- vtkwriter->write(time,Dune::VTK::ascii);
- }
-/*
+//------------------------------------------------------------------------------------
+
+
void grow(RF length) // How much the the solid-fluid interface should grow in normal direction (in m)
{
length = length/poreReferenceLength;
- if(length <= 0)
+ if(length > 0)
+ {
+ param.reaction.reactionRate = 1;
+ }
+ else
{
- std::cout << "Need positive length (dissolution not yet implemented)" << std::endl;
- return;
+ param.reaction.reactionRate = -1;
+ length = -length;
}
+
RF dtmax = param.pTree.get("Time.dtMax",0.01);
while(length > dtmax)
{
//writeVTK(1-length);
length -= dtmax;
param.time.dt = dtmax;
- applyCh(0);
+ applyChReaction(0);
}
param.time.dt = length;
- applyCh(0);
- applyNs(0);
- updateKf();
- updatePhiSolid();
+ applyChReaction(0);
}
void growArea(RF area) // How much Area the solid phase should grow (in m^2)
{
RF phiSolid_beforeGrowth = phiSolid;
area = area/(poreReferenceLength*poreReferenceLength*4); //4 because of Symmetry
- if(area <= 0)
+ int area_sign = 1;
+ if(area > 0)
{
- std::cout << "Need positive area (dissolution not yet implemented)" << std::endl;
- return;
+ param.reaction.reactionRate = 1;
}
+ else
+ {
+ param.reaction.reactionRate = -1;
+ area = -area;
+ area_sign = -1;
+ }
+
RF dtmax = param.pTree.get("Time.dtMax",0.01);
- while(area - (phiSolid - phiSolid_beforeGrowth) > dtmax*calculateGrowthSpeed())
+ while(area - area_sign*(phiSolid - phiSolid_beforeGrowth) > dtmax*calculateGrowthSpeed())
{
//writeVTK(1-length);
param.time.dt = dtmax;
- applyCh(0);
- updatePhiSolid();
+ applyChReaction(0);
}
- param.time.dt = (area - (phiSolid - phiSolid_beforeGrowth)) / calculateGrowthSpeed();
- applyCh(0);
- applyNs(0);
- updateKf();
- updatePhiSolid();
+ param.time.dt = (area - area_sign*(phiSolid - phiSolid_beforeGrowth)) / calculateGrowthSpeed();
+ applyChReaction(0);
+ }
+
+ void growAreaFraction(RF areaFraction)
+ {
+ updateSaturation();
+ growArea(areaFraction * phiFluid1);
+ }
+
+//------------------------------------------------------------------------------------
+
+void calculateKf(RF ViscosityRatio = 0)
+{
+ if(ViscosityRatio == 0)
+ {
+ ViscosityRatio = param.pTree.get("Physics.mu_ratio", (RF)1);
}
+ param.mu_ratio = ViscosityRatio;
+ applyNsTransmissibility();
+ updateKf();
+}
+
+void findSaturation(RF Curvature)
+{
+ //TODO!!!!
+}
- RF getKf()
+//------------------------------------------------------------------------------------
+
+ RF getKf_11()
{
- return kf;
+ return kf_11;
+ }
+
+ RF getKf_12()
+ {
+ return kf_12;
+ }
+
+ RF getKf_21()
+ {
+ return kf_21;
+ }
+
+ RF getKf_22()
+ {
+ return kf_22;
+ }
+
+ RF getPhiFluid1()
+ {
+ return phiFluid1;
+ }
+
+ RF getPhiFluid2()
+ {
+ return phiFluid2;
}
RF getPhiSolid()
@@ -560,6 +690,11 @@ public:
return (1-phiSolid)/(1-initialPhiSolid);
}
+ RF getSaturation()
+ {
+ return saturation;
+ }
+
void setPoreBaseRadius(RF PoreBaseRadius_)
{
poreReferenceLength = PoreBaseRadius_ / (1-param.initial.ptree.get("Radius",0.2) );
@@ -570,5 +705,4 @@ public:
return poreReferenceLength;
}
-*/
};
diff --git a/src/pntest/pn2p.ini b/src/pntest/pn2p.ini
index b955c96f961bb4b367f99d61e7e29b231489afac..ae77650ff893d6138084944998a11e082014b7e3 100644
--- a/src/pntest/pn2p.ini
+++ b/src/pntest/pn2p.ini
@@ -11,17 +11,20 @@ SurfaceTensionScaling=1
VelDissipationRate = 4000
CurvatureAlpha=0#0.001
-[test]
-InitialGrowthTime = 10
-FillTime = 150
-EmptyTime = 400
+[FullPCS_Curve]
+InitialGrowthIter = 10
+FillIter = 150
+EmptyIter = 400
+[PCS_Search]
+VelThreshold = 0.0005 #5e-4
+MaxCHNSIter = 100
[domain]
filename = grids/square_periodic.msh
[output]
-filename = outp3
+filename = outp_kf_testing
[Phasefield]
eps=0.02