vpfc.cc 7.87 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
#include "AMDiS.h"
#include "Helpers.h"
#include "base_problems/PhaseFieldCrystal.h"

#if (defined HAVE_SEQ_PETSC) || (defined HAVE_PETSC)
#include "preconditioner/PetscPreconPfc.h"
#elif !defined(HAVE_PARALLEL_DOMAIN_AMDIS)
#include "preconditioner/MTLPreconPfc.h"
#else
#include "preconditioner/PetscSolverPfc.h"
#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
98
99
100
101
102
103
104
  void finalizeData() override
  {
    super::finalizeData();

#if (defined HAVE_SEQ_PETSC) || (defined HAVE_PETSC)
    PetscPreconPfc* precon = dynamic_cast<PetscPreconPfc*>(prob->getSolver()->getRightPrecon());
    if (precon)
      precon->setData(getTau(), M0);
#elif !defined(HAVE_PARALLEL_DOMAIN_AMDIS)
    MTLPreconPfc* precon = dynamic_cast<MTLPreconPfc*>(prob->getSolver()->getRightPrecon());
    if (precon)
105
      precon->setData(getTau(), M0);
106
107
108
109
110
111
#else
    Parallel::PetscSolverPfc* solver = dynamic_cast<Parallel::PetscSolverPfc*>(prob->getSolver());
    if (solver)
      solver->setData(getTau(), M0);
#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

236
237
238
239
240
241
242
243
244
#if (defined HAVE_SEQ_PETSC) || (defined HAVE_PETSC)
  CreatorMap<PetscPreconditionerNested>::addCreator("pfc", new PetscPreconPfc::Creator);
#elif !defined(HAVE_PARALLEL_DOMAIN_AMDIS)
  CreatorMap<typename MTLPreconPfc::precon_base>::addCreator("pfc", new MTLPreconPfc::Creator);
#else
  CreatorMap<LinearSolverInterface>::addCreator("p_petsc_pfc", new Parallel::PetscSolverPfc::Creator);
#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;
}