Commit 3efd88e3 authored by Praetorius, Simon's avatar Praetorius, Simon

Merge branch 'example/cahn_hilliard' into 'develop'

Example/cahn hilliard

See merge request spraetor/dune-amdis!102
parents d5cac81c f02bb6f6
Pipeline #1834 passed with stage
in 26 minutes and 49 seconds
......@@ -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)
#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;
}
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
......@@ -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
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment