vpfc.cc 7.81 KB
Newer Older
1 2 3 4
#include "AMDiS.h"
#include "Helpers.h"
#include "base_problems/PhaseFieldCrystal.h"

Praetorius, Simon's avatar
Praetorius, Simon committed
5
#if defined(USE_SEQ_PETSC)
6
#include "preconditioner/PetscPreconPfc.h"
Praetorius, Simon's avatar
Praetorius, Simon committed
7
#elif defined(USE_PARALLEL_PETSC)
8
#include "preconditioner/PetscSolverPfc.h"
Praetorius, Simon's avatar
Praetorius, Simon committed
9 10
#elif defined(USE_MTL)
#include "preconditioner/MTLPreconPfc.h"
11 12 13 14 15 16 17 18 19 20 21 22 23 24
#endif

// #include "OneModeApproximation.h"
#include "Expressions.h"

#include "positions.hpp"
#include "maxima.hpp"

using namespace AMDiS;


struct Sinus : AbstractFunction<double, WorldVector<double> >
{
  Sinus(double density, double amplitude)  : AbstractFunction<double, WorldVector<double> >(4), density(density), amplitude(amplitude) {}
25

26 27
  double operator()(const WorldVector<double>& x) const
  {
28
    return density + 0.5*amplitude*(sin(162.0*sqr(x[0])+6123*x[1])*cos(243*sin(x[0])+152*x[0]*x[1]));
29 30 31 32 33 34 35 36
  }
private:
  double density, amplitude;
};


struct VPfcPeak : AbstractFunction<double, WorldVector<double> >
{
37
  VPfcPeak(WorldVector<double> pos_, double A_, double d_=4.0*m_pi/sqrt(3.0))
38 39 40 41 42 43
    : AbstractFunction<double, WorldVector<double> >(6),pos(pos_),A(A_),d(d_)
  {
    q = (4.0*m_pi/sqrt(3.0))/d;
  }

  double operator()(const WorldVector<double> &x) const
44
  {
45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
    WorldVector<double> x1; x1 = x-pos;
    double nrmX = sqrt(x1*x1);
    return (nrmX < d/2.0 ? A*(cos(q*sqrt(3.0)/2.0 * nrmX)+1.0) : 0.0);
  }
protected:
  WorldVector<double> pos;
  double A;
  double d;
  double q;
};

class VPfcPC : public base_problems::detail::PhaseFieldCrystal<ExtendedProblemStat>
{
  typedef base_problems::detail::PhaseFieldCrystal<ExtendedProblemStat> super;

public:
61 62 63 64 65
  VPfcPC(std::string name)
    : super(name),
      overall_solution_time(0.0),
      overall_build_time(0.0),
      overall_solver_iter(0),
66 67 68 69 70 71 72 73 74 75 76 77 78
      n_timesteps(0),
      H(1500.0),
      B0(10.0),
      B1(1.0),
      rho1(1.0),
      A(1.0),
      d(4.0*m_pi/sqrt(3.0)),
      max_psi(1.0)
  {
    Parameters::get(name + "->H", H);
    Parameters::get(name + "->B0",B0);
    Parameters::get(name + "->rho1",rho1);
  }
79

80 81 82
  //////////////////////////////////////////////////////////////////////////////////////////
  void initData()
  { FUNCNAME("VPfcPC::solveInitialProblem()");
83

84 85 86 87 88 89 90 91 92
    super::initData();
    initPositions(name, pos);
    MSG("Anzahl Partikel: %d\n", pos.size());

    double Bp =  M_PI*sqr(2.0*M_PI/std::sqrt(3.0));
    B1 = pos.size() * Bp;
    A = rho1;
    density = B1/B0*rho1;
  }
93

94 95 96 97
  void finalizeData() override
  {
    super::finalizeData();

Praetorius, Simon's avatar
Praetorius, Simon committed
98
#if defined(USE_SEQ_PETSC)
99 100 101
    PetscPreconPfc* precon = dynamic_cast<PetscPreconPfc*>(prob->getSolver()->getRightPrecon());
    if (precon)
      precon->setData(getTau(), M0);
Praetorius, Simon's avatar
Praetorius, Simon committed
102
#elif defined(USE_PARALLEL_PETSC)
103 104 105
    Parallel::PetscSolverPfc* solver = dynamic_cast<Parallel::PetscSolverPfc*>(prob->getSolver());
    if (solver)
      solver->setData(getTau(), M0);
Praetorius, Simon's avatar
Praetorius, Simon committed
106 107 108 109
#elif defined(USE_MTL)
    MTLPreconPfc* precon = dynamic_cast<MTLPreconPfc*>(prob->getSolver()->getRightPrecon());
    if (precon)
      precon->setData(getTau(), M0);
110 111
#endif
  }
112

113 114 115 116 117 118 119
  // generate initial solution for evolution equation
  void solveInitialProblem(AdaptInfo *adaptInfo) override
  { FUNCNAME("VPfcPC::solveInitialProblem()");

    Flag initFlag = initDataFromFile(adaptInfo);
    if (initFlag.isSet(DATA_ADOPTED))
      return;
120

121
    DOFVector<double>* psi = getDensity();
122

123 124 125
    psi->set(0.0);
    for (size_t i = 0; i < pos.size(); i++)
      *psi << valueOf(psi) + eval(new VPfcPeak(pos[i], A));
126

127 128 129
    double mean_density = psi->Int() / B0;
    MSG("rho0: %e, mean density: %e, rho1: %e\n", density, mean_density, rho1);
    *psi << valueOf(psi) * (density/mean_density);
130

131
    writeMaxima(name, adaptInfo, psi, pos);
132

133
    MSG("Number of degrees of freedom: %d\n", getDensity()->getUsedSize());
134

135 136 137 138
    double minH, maxH;
    int minLevel, maxLevel;
    MSG("Mesh size: %d\n", Helpers::calcMeshSizes(prob->getMesh(), minH, maxH, minLevel, maxLevel));
  }
139

140 141 142 143
  void initTimestep(AdaptInfo *adaptInfo) override
  {
    using namespace AMDiS::io;
    super::initTimestep(adaptInfo);
144

145 146 147 148 149
    bool write_matrix = false;
    Parameters::get("write matrix",write_matrix);
    if (write_matrix) {
      std::string xFilename="x.dat";
      Parameters::get("x filename", xFilename);
150

151 152 153 154 155
      std::vector<DOFVector<double>*> vecs;
      for (int i = 0; i < prob->getNumComponents(); i++)
	vecs.push_back(prob->getSolution(i));
      DofWriter::writeFile(vecs, xFilename);
    }
156

157 158
    max_psi = max(*getDensity());
  }
159

160 161 162 163 164 165 166 167 168
  void closeTimestep(AdaptInfo *adaptInfo) override
  {
    using namespace AMDiS::io;
    super::closeTimestep(adaptInfo);

    n_timesteps++;
    overall_solution_time += prob->getSolutionTime();
    overall_build_time += prob->getBuildTime();
    overall_solver_iter += adaptInfo->getSolverIterations();
169

170 171 172 173 174 175 176
    bool write_matrix = false;
    Parameters::get("write matrix",write_matrix);
    if (write_matrix) {
      std::string matrixFilename="matrix.mtx", rhsFilename="rhs.dat";
      Parameters::get("matrix filename", matrixFilename);
      Parameters::get("rhs filename", rhsFilename);
      prob->writeMatrix(matrixFilename);
177

178 179 180 181 182 183 184 185 186 187
      std::vector<DOFVector<double>*> vecs;
      for (int i = 0; i < prob->getNumComponents(); i++)
	vecs.push_back(prob->getRhsVector(i));
      DofWriter::writeFile(vecs, rhsFilename);
    }
  }

  void fillOperators() override
  {
    super::fillOperators();
188

189 190
    DOFVector<double>* psi = getDensity();
    static const int n = 3;
191

192 193
    Operator *opPenalty = new Operator(prob->getFeSpace(0));
    addZOT(opPenalty, (H*n*(2.0 - n)) * (signum(valueOf(psi)) - 1.0) * pow<n-1>(valueOf(psi)) );
194 195
    prob->addVectorOperator(opPenalty, 0);

196 197
    Operator *opPenaltyDerivLhs = new Operator(prob->getFeSpace(0), prob->getFeSpace(1));
    addZOT(opPenaltyDerivLhs, (-H*n*(n-1)) * (signum(valueOf(psi)) - 1.0) * pow<n-2>(valueOf(psi)) );
198 199 200
    prob->addMatrixOperator(opPenaltyDerivLhs, 0, 1);

    if (mobility_type == 3) {
201 202 203 204 205
      Operator *opLM = new Operator(prob->getFeSpace(1), prob->getFeSpace(0));
      addSOT(opLM, max(valueOf(psi)*(M0 / ref_(&max_psi)), 1.e-2));
      prob->addMatrixOperator(opLM, 1, 0, getTau(), getTau()); // -laplace(mu)
    }
  }
206 207


208 209 210 211 212 213 214 215 216 217
  double getSolutionTime() { return overall_solution_time; }
  double getBuildTime() { return overall_build_time; }
  int getSolverIterations() { return overall_solver_iter; }
  int getNumTimesteps() { return n_timesteps; }

private:
  double overall_solution_time;
  double overall_build_time;
  int overall_solver_iter;
  int n_timesteps;
218

219 220 221 222 223 224 225
  std::vector<WorldVector<double> > pos;

  double H;
  double B0, B1;
  double rho1;
  double A;
  double d;
226

227 228 229 230 231 232 233 234
  double max_psi;
};


int main(int argc, char** argv)
{ FUNCNAME("main");

  AMDiS::init(argc, argv);
235

Praetorius, Simon's avatar
Praetorius, Simon committed
236
#if defined(USE_SEQ_PETSC)
237
  CreatorMap<PetscPreconditionerNested>::addCreator("pfc", new PetscPreconPfc::Creator);
Praetorius, Simon's avatar
Praetorius, Simon committed
238
#elif defined(USE_PARALLEL_PETSC)
239
  CreatorMap<LinearSolverInterface>::addCreator("p_petsc_pfc", new Parallel::PetscSolverPfc::Creator);
Praetorius, Simon's avatar
Praetorius, Simon committed
240 241
#elif defined(USE_MTL)
  CreatorMap<typename MTLPreconPfc::precon_base>::addCreator("pfc", new MTLPreconPfc::Creator);
242 243 244
#endif

  Timer t;
245

246 247 248 249 250 251 252 253 254 255 256 257 258 259 260
  VPfcPC pfcProb("vpfc");
  pfcProb.initialize(INIT_ALL);

  // Adapt-Infos
  AdaptInfo adaptInfo("adapt", pfcProb.getNumComponents());
  AdaptInstationary adaptInstat("adapt", pfcProb, adaptInfo, pfcProb, adaptInfo);

  // Scale Mesh
  bool scaleMesh = false;
  Initfile::get("mesh->scale mesh",scaleMesh);
  if (scaleMesh) {
    WorldVector<double> scale; scale.set(1.0);
    Initfile::get("mesh->dimension",scale);
    Helpers::scaleMesh(pfcProb.getMesh(), scale);
  }
261

262
  pfcProb.initTimeInterface(); // fill operators and BC
263
  int error_code = adaptInstat.adapt();
264 265 266 267 268 269

  MSG("elapsed time= %f sec\n", t.elapsed());

  MSG("solution time= %f sec\n", pfcProb.getSolutionTime() / pfcProb.getNumTimesteps() );
  MSG("build time= %f sec\n", pfcProb.getBuildTime() / pfcProb.getNumTimesteps() );
  MSG("solver iterations= %f\n", static_cast<double>(pfcProb.getSolverIterations()) / pfcProb.getNumTimesteps() );
270

271 272 273
  AMDiS::finalize();
  return error_code;
}