Different behaviour in usage of (1 - solution) and (1.0 - solution) while defining operators
I have a strange problem in the operators. For the below code, when I define opMob in line 39 with (1 - phi), there is no evolution observed for the cahn-hilliard model and if I change it to (1.0 - phi), the model behaves as expected, where phi is one of the solution variables. Attaching the macro and init files along. macro.square.2dSurface_Diffusion.dat
#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>
#include <dune/grid/albertagrid.hh>
using namespace AMDiS;
using Grid = Dune::AlbertaGrid<GRIDDIM, WORLDDIM>;
using Param = LagrangeBasis<Grid, 1, 1>;
int main(int argc, char** argv)
{
Environment env(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);
auto mu = prob.solution(1);
double M = Parameters::get<double>("parameters->mobility").value_or(1.0);
double eps = Parameters::get<double>("parameters->epsilon").value_or(0.02);
prob.addMatrixOperator(makeOperator(tag::test_trial{}, invTau), 0, 0);
prob.addVectorOperator(makeOperator(tag::test{}, phiOld * invTau), 0);
auto opMob = (M/eps) * 36 * pow<2>(phi) * pow<2>(1 - phi);
prob.addMatrixOperator(makeOperator(tag::gradtest_gradtrial{}, opMob), 0, 1);
prob.addMatrixOperator(makeOperator(tag::test_trial{}, 1.0), 1, 1);
prob.addMatrixOperator(makeOperator(tag::gradtest_gradtrial{}, -eps), 1, 0);
auto opFimpl = makeOperator(tag::test_trial{}, -1/eps * (36 - 216 * phi + 216 * pow<2>(phi)));
prob.addMatrixOperator(opFimpl, 1, 0);
auto opFexpl1 = makeOperator(tag::test{}, 1/eps * (36 * phi - 108 * pow<2>(phi) + 72 * pow<3>(phi)));
prob.addVectorOperator(opFexpl1, 1);
auto opFexpl2 = makeOperator(tag::test{}, -phi/eps * (36 - 216 * phi + 216 * pow<2>(phi)));
prob.addVectorOperator(opFexpl2, 1);
int ref_int = Parameters::get<int>("refinement->interface").value_or(10);
int ref_bulk = Parameters::get<int>("refinement->bulk").value_or(2);
GridFunctionMarker marker("interface", prob.grid(),
invokeAtQP([ref_int, ref_bulk](double xx) {
return xx > 0.05 && xx < 0.95 ? ref_int : ref_bulk;
}, phi));
prob.addMarker(marker);
for (int i = 0; i < 6; ++i) {
phi << [eps](auto const& x) {
using Math::sqr;
double theta = std::atan2(x[1], x[0]);
double r = (0.6 * 0.4) / std::sqrt(std::pow(0.4 * std::cos(theta), 2) + std::pow(0.6 * std::sin(theta),2));
double dist = std::sqrt(std::pow(x[0],2) + std::pow(x[1],2));
return 0.5 * (1 - std::tanh(3.0 * ((std::abs(dist) - r) / eps)));
};
prob.markElements(adaptInfo);
prob.adaptGrid(adaptInfo);
}
AdaptInstationary adapt("adapt", prob, adaptInfo, probInstat, adaptInfo);
adapt.adapt();
return 0;
}