Example/cahn hilliard

examples/CMakeLists.txt

@@ -18,10 +18,12 @@ add_amdis_executable(NAME stokes1.2d SOURCES stokes1.cc DIM 2 DOW 2)
 add_amdis_executable(NAME stokes3.2d SOURCES stokes3.cc DIM 2 DOW 2)
 add_amdis_executable(NAME navier_stokes.2d SOURCES navier_stokes.cc DIM 2 DOW 2)
 add_amdis_executable(NAME convection_diffusion.2d SOURCES convection_diffusion.cc DIM 2 DOW 2)
+add_amdis_executable(NAME cahn_hilliard.2d SOURCES cahn_hilliard.cc DIM 2 DOW 2 ALBERTA_GRID)
 add_dependencies(examples
   vecellipt.2d
   stokes0.2d
   stokes1.2d
   stokes3.2d
   navier_stokes.2d
-  convection_diffusion.2d)
+  convection_diffusion.2d
+  cahn_hilliard.2d)

examples/cahn_hilliard.cc

+#include <amdis/AMDiS.hpp>
+#include <amdis/AdaptInstationary.hpp>
+#include <amdis/LocalOperators.hpp>
+#include <amdis/ProblemInstat.hpp>
+#include <amdis/ProblemStat.hpp>
+#include <amdis/GridFunctions.hpp>
+#include <amdis/Marker.hpp>
+
+using namespace AMDiS;
+
+using Grid = Dune::AlbertaGrid<GRIDDIM, WORLDDIM>;
+using Param = LagrangeBasis<typename Grid::LeafGridView, 1, 1>;
+
+int main(int argc, char** argv)
+{
+  AMDiS::init(argc, argv);
+
+  ProblemStat<Param> prob("ch");
+  prob.initialize(INIT_ALL);
+
+  ProblemInstat<Param> probInstat("ch", prob);
+  probInstat.initialize(INIT_UH_OLD);
+
+  AdaptInfo adaptInfo("adapt");
+
+  auto invTau = std::ref(probInstat.invTau());
+  auto phi = prob.solution(0);
+  auto phiOld = probInstat.oldSolution(0);
+
+  prob.addMatrixOperator(zot(invTau), 0, 0);
+  prob.addVectorOperator(zot(phiOld * invTau), 0);
+
+  double M = Parameters::get<double>("parameters->mobility").value_or(1.0);
+  prob.addMatrixOperator(sot(M), 0, 1);
+  prob.addMatrixOperator(zot(1.0), 1, 1);
+
+  double a = Parameters::get<double>("parameters->a").value_or(1.0);
+  double b = Parameters::get<double>("parameters->b").value_or(1.0/4.0);
+
+  double eps = Parameters::get<double>("parameters->epsilon").value_or(0.02);
+  prob.addMatrixOperator(sot(-a*eps), 1, 0);
+
+  auto opFimpl = zot(-b/eps * (2 + 12*phi*(phi - 1)));
+  prob.addMatrixOperator(opFimpl, 1, 0);
+
+  auto opFexpl = zot(b/eps * pow<2>(phi)*(6 - 8*phi));
+  prob.addVectorOperator(opFexpl, 1);
+
+  int ref_int  = Parameters::get<int>("refinement->interface").value_or(10);
+  int ref_bulk = Parameters::get<int>("refinement->bulk").value_or(2);
+  auto marker = makeGridFunctionMarker("interface", prob.grid(),
+    invokeAtQP([ref_int, ref_bulk](double phi) {
+      return phi > 0.05 && phi < 0.95 ? ref_int : ref_bulk;
+    }, phi));
+  prob.addMarker(marker);
+
+  double radius1 = Parameters::get<double>("parameters->radius1").value_or(0.15);
+  double radius2 = Parameters::get<double>("parameters->radius2").value_or(0.25);
+  for (int i = 0; i < 6; ++i) {
+    phi << [eps,radius1,radius2](auto const& x) {
+      using Math::sqr;
+      return 0.5*(1 - std::tanh((radius1+radius2)*(std::sqrt(sqr(x[0]/radius1) + sqr(x[1]/radius2)) - 1.0)/(4*std::sqrt(2.0)*eps)));
+    };
+    prob.markElements(adaptInfo);
+    prob.adaptGrid(adaptInfo);
+  }
+
+  AdaptInstationary adapt("adapt", prob, adaptInfo, probInstat, adaptInfo);
+  adapt.adapt();
+
+  AMDiS::finalize();
+  return 0;
+}

examples/init/cahn_hilliard.dat.2d

+dimension of world:   2
+
+chMesh->macro file name:    ./macro/macro.square.2d
+chMesh->global refinements: 8
+
+parameters->radius:  0.5
+parameters->radius1: 0.4
+parameters->radius2: 0.6
+parameters->epsilon: 0.01
+
+refinement->interface: 12
+refinement->bulk: 2
+
+ch->mesh:                       chMesh
+ch->solver->name:               direct
+ch->solver->max iteration:      1000
+ch->solver->absolute tolerance: 1e-6
+ch->solver->break if tolerance not reached: 1
+ch->solver->info: 2
+
+ch->output[0]->filename:           phi.2d
+ch->output[0]->output directory:   output
+ch->output[0]->name:               phi
+ch->output[0]->ParaView format:    1
+ch->output[0]->ParaView mode:      1
+ch->output[0]->ParaView animation: 1
+
+ch->output[1]->filename:           mu.2d
+ch->output[1]->output directory:   output
+ch->output[1]->name:               mu
+ch->output[1]->ParaView format:    1
+ch->output[1]->ParaView mode:      1
+ch->output[1]->ParaView animation: 1
+
+adapt->timestep:    0.001
+adapt->start time:  0.0
+adapt->end time:    1.0

examples/macro/macro.square.2d

@@ -5,15 +5,15 @@ number of elements: 8
 number of vertices: 9
 
 vertex coordinates:
--50 -50
-0 -50
-50 -50
--50 0
+-1.0 -1.0
+0 -1.0
+1.0 -1.0
+-1.0 0
 0 0
-50 0
--50 50
-0 50
-50 50
+1.0 0
+-1.0 1.0
+0 1.0
+1.0 1.0
 
 element vertices:
 0 4 3