Commit c2bd24db authored by Lasse Hinrichsen's avatar Lasse Hinrichsen
Browse files

Example: Allen-Cahn with log potential

parent e09fb0ed
Pipeline #34655 passed with stage
in 5 minutes and 1 second
if (dune-pdelab_FOUND)
set(programs contact
set(programs allencahnlog
contact
obstacle
poisson
trescafriction)
......
#include <config.h>
#include <memory>
#include <dune/common/exceptions.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
#include <dune/grid/io/file/vtk.hh>
#include <dune/grid/utility/structuredgridfactory.hh>
#include <dune/grid/yaspgrid.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
#include <dune/solvers/common/defaultbitvector.hh>
#include <dune/solvers/iterationsteps/multigridstep.hh>
#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/solvers/loopsolver.hh>
#include <dune/solvers/test/common.hh>
#include <dune/solvers/transferoperators/truncateddensemgtransfer.hh>
#include <dune/tnnmg/functionals/bcqfconstrainedlinearization.hh>
#include <dune/tnnmg/functionals/boxconstrainedquadraticfunctional.hh>
#include <dune/tnnmg/functionals/quadraticfunctional.hh>
#include <dune/tnnmg/functionals/quadraticfunctionalconstrainedlinearization.hh>
#include <dune/tnnmg/functionals/semilinearfunctional.hh>
#include <dune/tnnmg/functionals/semilinearfunctionallinearization.hh>
#include <dune/tnnmg/functionals/sumfunctional.hh>
#include <dune/tnnmg/functionals/sumfunctionalconstrainedlinearization.hh>
#include <dune/tnnmg/iterationsteps/nonlineargsstep.hh>
#include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
#include <dune/tnnmg/localsolvers/functionproperties.hh>
#include <dune/tnnmg/localsolvers/scalarbisectionsolver.hh>
#include <dune/tnnmg/localsolvers/scalarobstaclesolver.hh>
#include <dune/tnnmg/projections/obstacledefectprojection.hh>
#include <dune/tnnmg/projections/semilineardefectprojection.hh>
template<typename GV>
auto
stiffness(GV const& gridview)
{
using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
auto matrix = Matrix();
constructPQ1Pattern(gridview, matrix);
matrix = 0.0;
assemblePQ1Stiffness(gridview, matrix);
return matrix;
}
template<typename GV>
auto
mass(GV const& gridview)
{
using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
auto matrix = Matrix();
constructPQ1Pattern(gridview, matrix);
matrix = 0.0;
assemblePQ1Mass(gridview, matrix);
return matrix;
}
template<typename Grid,
typename MatrixType,
typename VectorType,
typename BitVector>
auto
multigridSolver(Grid const& grid)
{
// copied code from obstacle.cc:
// no base solver for now, lets see how far this gets us
using namespace Dune;
auto smoother =
std::make_shared<TruncatedBlockGSStep<MatrixType, VectorType, BitVector>>();
auto linearMultigridStep =
std::make_shared<Solvers::MultigridStep<MatrixType, VectorType>>();
linearMultigridStep->setMGType(1, 3, 3);
linearMultigridStep->setSmoother(smoother);
linearMultigridStep->mgTransfer_.resize(grid.maxLevel());
for (size_t i = 0; i < linearMultigridStep->mgTransfer_.size(); i++) {
// create transfer operator from level i to i+1
auto transferOperator =
std::make_shared<TruncatedDenseMGTransfer<VectorType>>();
transferOperator->setup(grid, i, i + 1);
linearMultigridStep->mgTransfer_[i] = transferOperator;
}
return linearMultigridStep;
}
// supply the potential. This could also be in a more prominent place to
// highlight its interface BEGIN POTENTIAL
// global variables, yay!
const auto temp_theta = 0.1;
struct LogPotentialDerivative
{
double operator()(double x)
{
if (-1 < x and x < 1.)
return temp_theta / 2 * (std::log(1 + x) - std::log(1. - x));
return 0.; // will be truncated
}
};
struct LogPotentialDerivative2
{
double operator()(double x)
{
if (-1 < x and x < 1.)
return temp_theta / 2 * (-2. / (x * x - 1));
return 0.; // will be truncated
}
};
class LogPotential;
auto
derivative(const LogPotential&)
{
return LogPotentialDerivative();
}
auto
derivative(const LogPotentialDerivative&)
{
return LogPotentialDerivative2{};
}
class LogPotential
{
double eps = 1e-12;
double R = 1e5; // sec. derivatives should not pass this bound
public:
using Interval = Dune::Solvers::Interval<double>;
Interval domain_ = { -1 , 1. };
auto operator()(double x) const
{
return temp_theta / 2 *
((1. + x) * std::log(1. + x) + (1. - x) * std::log(1. - x));
}
auto subDifferential(double x) const
{
if (-1 < x and x < 1.) {
return Interval(temp_theta / 2 * (std::log(1 + x) - std::log(1. - x)));
}
// TODO is this correct?
else if (x <= -1.)
return Interval(std::numeric_limits<double>::lowest());
else
return Interval(std::numeric_limits<double>::max());
}
auto domain() const { return domain_; }
bool should_truncate(double x) const
{
if (domain_[0] + eps >= x or x >= domain_[1] - eps) // safety guards on
return true;
if (std::abs(LogPotentialDerivative2()(x))>R)
return true;
return false;
}
};
//// END POTENTIAL
template<typename V, typename M, typename Grid>
auto
solveSystem(V const& u_old, M const& matrix, Grid const& grid)
{
using namespace Dune;
auto u = u_old;
// u = 0.2; // for the lulz DEBUG
// set up ignore vector. We have no Dirichlet nodes, hence all false.
auto ignore = Solvers::DefaultBitVector_t<V>(u_old.size(), false);
using BV = decltype(ignore);
// set up weights (integrate the constant 1)
auto weights = u_old;
{
weights = 0;
assemblePQ1RHS(grid.leafGridView(),
weights,
[](auto) -> FieldVector<double, 1> { return 1.0; });
}
// set up functionals
///* DEBUG: Use obstacle potential
auto phi = LogPotential();
auto quadratic = TNNMG::QuadraticFunctional<M, V, double>(matrix, u_old);
auto nonlinear = TNNMG::Semilinearfunctional<V, decltype(phi)>(weights, phi);
auto functional =
TNNMG::SumFunctional<decltype(quadratic), decltype(nonlinear)>(quadratic,
nonlinear);
// now, set up linearizations
auto quad_lin =
TNNMG::QuadraticFunctionalConstrainedLinearization<decltype(quadratic), BV>(
quadratic, ignore);
auto semi_lin =
TNNMG::SemilinearFunctionalConstrainedLinearization<decltype(nonlinear),
LogPotential,
M,
V,
BV>(nonlinear, ignore);
semi_lin.setMatrixDummy(matrix);
using Linearization =
TNNMG::SumFunctionalConstrainedLinearization<decltype(functional),
BV,
decltype(quad_lin),
decltype(semi_lin)>;
auto functional_linearization = std::make_shared<Linearization>(
std::make_tuple(quad_lin, semi_lin), ignore);
// */ // End Debug obstacle
// DEBUG OBSTACLE
// V lower(u.size());
// V upper(u.size());
// lower = -1;
// upper = 1;
// using Functional =
// TNNMG::BoxConstrainedQuadraticFunctional<M, V, V, V, double>;
// auto functional = Functional(matrix, u_old, lower, upper);
// using Linearization =
// TNNMG::BoxConstrainedQuadraticFunctionalConstrainedLinearization<Functional,
// BV>;
// using DefectProjection = TNNMG::ObstacleDefectProjection;
// using LineSearchSolver = TNNMG::ScalarObstacleSolver;
// END DEBUG OBSTACLE
// set up TNNMGStep
using LocalSolver = TNNMG::ScalarBisectionSolver;
auto ls = LocalSolver();
// auto ls = TNNMG::ScalarObstacleSolver();
auto gs = gaussSeidelLocalSolver(ls);
using GS = decltype(gs);
using Smoother = TNNMG::NonlinearGSStep<decltype(functional), GS>;
auto proj = TNNMG::SemilinearFunctionalDefectProjection();
// auto linesearch = ls; // probably dont even need to copy
auto linesearch = LocalSolver(1.0, 0.8);
auto smoother = std::make_shared<Smoother>(functional, u, gs);
smoother->setIgnore(ignore);
auto mg = multigridSolver<Grid, M, V, BV>(grid);
using Step = TNNMG::TNNMGStep<decltype(functional),
BV,
Linearization,
decltype(proj),
decltype(linesearch)>;
auto step = Step(functional, u, smoother, mg, 1, proj, linesearch);
// using Step = TNNMG::TNNMGStep<Functional,
// BV,
// Linearization,
// DefectProjection,
// LineSearchSolver>;
// auto step = Step(
// functional, u, smoother, mg, 1, DefectProjection(), LineSearchSolver());
step.setIgnore(ignore);
step.setLinearization(functional_linearization);
// solve
using Solver = Dune::Solvers::LoopSolver<V, BV>;
// TODO: make these available
int iter = 100;
double tol = 1e-12;
using Norm = Solvers::EnergyNorm<M, V>;
// TODO: dont assemble each time
auto stiff = stiffness(grid.leafGridView());
auto norm = std::make_shared<Norm>(stiff);
auto solver = Solver(step, iter, tol, norm, Solver::FULL);
// auto solver = Solver(smoother, iter, tol, norm, Solver::FULL);
///* Only relevant for actual TNNMG
solver.addCriterion(
[&]() {
return Dune::formatString(" % 12.5e", step.lastDampingFactor());
},
" damping ");
solver.addCriterion(
[&]() {
return Dune::formatString(" % 12d",
step.linearization().truncated().count());
},
" truncated ");
// */
solver.preprocess();
solver.solve();
return u;
}
template<typename GridView, typename V>
void
plot(const GridView& gv, const V& sol, int timestep)
{
std::string filename = "allen_cahn_solution";
std::ostringstream oss;
oss << std::setfill('0') << std::setw(4) << timestep;
filename += oss.str();
auto basis = Dune::Functions::LagrangeBasis<GridView, 1>(gv);
auto xFunction =
Dune::Functions::makeDiscreteGlobalBasisFunction<double>(basis, sol);
Dune::VTKWriter<GridView> vtkWriter(basis.gridView());
vtkWriter.addVertexData(
xFunction,
Dune::VTK::FieldInfo("x", Dune::VTK::FieldInfo::Type::scalar, 1));
vtkWriter.write(filename);
}
int
main(int argc, char** argv)
{
using namespace Dune;
MPIHelper::instance(argc, argv);
// Model parameters:
auto eps = 0.1; // width of the smeared interface
auto tau = 0.005; // time step size
auto num_timesteps = 5;
// Discretization parameters:
auto gridlevels = 6;
// set up a grid
constexpr const int dim = 2;
using Grid = YaspGrid<dim>;
auto gridptr = StructuredGridFactory<Grid>::createCubeGrid(
{ { 0.0, 0.0 } }, { { 1.0, 1.0 } }, { { 2, 2 } });
gridptr->globalRefine(gridlevels - 1);
auto gridview = gridptr->leafGridView();
// prepare system matrix
auto matrix = stiffness(gridview);
auto massMatrix = mass(gridview);
{
matrix *= eps * eps;
auto m = massMatrix;
m *= (eps * eps / tau - 1.);
matrix += m;
}
using Vector = BlockVector<FieldVector<double, 1>>;
// Prepare RHS. Set up some initial function u_0 and integrate it.
auto u0 = Vector(matrix.M());
{
u0 = 0.0;
auto u0_lambda = [](auto const& pos) -> FieldVector<double, 1> {
return pos * pos > 0.5 ? -1.0 : 1.0;
};
assemblePQ1RHS(gridview, u0, u0_lambda);
}
auto solutions = std::vector<Vector>(num_timesteps);
solutions[0] = u0; // TODO: This is not the function but the functional.
for (int timestep = 1; timestep < num_timesteps; ++timestep) {
auto u_old = Vector();
if (timestep == 1) {
u_old = u0;
} else {
u_old = solutions[timestep - 1];
massMatrix.mv(solutions[timestep - 1], u_old);
}
u_old *= eps * eps / tau;
solutions[timestep] = solveSystem(u_old, matrix, *gridptr);
}
// print vtk files
for (int t = 1; t < solutions.size(); ++t) {
plot(gridview, solutions[t], t);
}
}
\ No newline at end of file
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