pntest.cc

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include<iostream>
#include<fstream>
// dune includes
#include<dune/common/parallel/mpihelper.hh>
#include<dune/common/parametertreeparser.hh>
#include<dune/common/timer.hh>
#if HAVE_DUNE_ALUGRID
#include<dune/alugrid/grid.hh>
#endif
// phasefield includes
#include<dune/phasefield/porenetwork/pn_1pflow.hh>

//===============================================================
// Main program
//===============================================================


int main(int argc, char** argv)
{
  try{
    // Maybe initialize Mpi
    Dune::MPIHelper&
      helper = Dune::MPIHelper::instance(argc, argv);
    if(Dune::MPIHelper::isFake)
      std::cout<< "This is a sequential program." << std::endl;
    else
      std::cout << "Parallel code run on "
                << helper.size() << " process(es)" << std::endl;

    // open ini file
    Dune::ParameterTree ptree;
    Dune::ParameterTreeParser ptreeparser;
    ptreeparser.readINITree("pn1p.ini",ptree);
    ptreeparser.readOptions(argc,argv,ptree);

    // Define Template Options
    typedef double RF;
    typedef Params_fs_r< RF,
                          //DoubleWell_limited< RF, DoubleWell_poly<RF>, Limiter_FakeDiverging<RF> >,
                          DoubleWell_limited< RF, DoubleWell_8th_order_poly<RF>, Limiter_FakeDiverging<RF> >,
                          Mobility_Const<RF>,
                          Reaction_Const<RF>,
                          VelDissipation_Quadratic<RF>
                          //VelDissipation_Quadratic_Shifted<RF>
                          > 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;
    //param.initial.shape = INITIAL_SHAPE_ELLIPSE;
    param.initial.shape = INITIAL_SHAPE_TWOCIRCLES;

    //Create Simulation
    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 << "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");
    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;

    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;

  }

  }
  catch (Dune::Exception &e){
    std::cerr << "Dune reported error: " << e << std::endl;
    return 1;
  }
  catch (...){
    std::cerr << "Unknown exception thrown!" << std::endl;
    return 1;
  }
}