vpfc.cc 7.98 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
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
#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) {}
  
  double operator()(const WorldVector<double>& x) const
  {
    return density + 0.5*amplitude*(sin(162.0*sqr(x[0])+6123*x[1])*cos(243*sin(x[0])+152*x[0]*x[1]));    
  }
private:
  double density, amplitude;
};


struct VPfcPeak : AbstractFunction<double, WorldVector<double> >
{
  VPfcPeak(WorldVector<double> pos_, double A_, double d_=4.0*m_pi/sqrt(3.0)) 
    : 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
  {    
    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:
  VPfcPC(std::string name) 
    : super(name), 
      overall_solution_time(0.0), 
      overall_build_time(0.0), 
      overall_solver_iter(0), 
      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);
  }
  
  //////////////////////////////////////////////////////////////////////////////////////////
  void initData()
  { FUNCNAME("VPfcPC::solveInitialProblem()");
  
    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;
  }
  
  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)
      precon->setData(getTau(), M0);	
#else
    Parallel::PetscSolverPfc* solver = dynamic_cast<Parallel::PetscSolverPfc*>(prob->getSolver());
    if (solver)
      solver->setData(getTau(), M0);
#endif
  }
  
  // generate initial solution for evolution equation
  void solveInitialProblem(AdaptInfo *adaptInfo) override
  { FUNCNAME("VPfcPC::solveInitialProblem()");

    Flag initFlag = initDataFromFile(adaptInfo);
    if (initFlag.isSet(DATA_ADOPTED))
      return;
      
    DOFVector<double>* psi = getDensity();
    
    psi->set(0.0);
    for (size_t i = 0; i < pos.size(); i++)
      *psi << valueOf(psi) + eval(new VPfcPeak(pos[i], A));
    
    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);
    
    writeMaxima(name, adaptInfo, psi, pos);
    
    MSG("Number of degrees of freedom: %d\n", getDensity()->getUsedSize());
	      
    double minH, maxH;
    int minLevel, maxLevel;
    MSG("Mesh size: %d\n", Helpers::calcMeshSizes(prob->getMesh(), minH, maxH, minLevel, maxLevel));
  }
  
  void initTimestep(AdaptInfo *adaptInfo) override
  {
    using namespace AMDiS::io;
    super::initTimestep(adaptInfo);
    
    bool write_matrix = false;
    Parameters::get("write matrix",write_matrix);
    if (write_matrix) {
      std::string xFilename="x.dat";
      Parameters::get("x filename", xFilename);
      
      std::vector<DOFVector<double>*> vecs;
      for (int i = 0; i < prob->getNumComponents(); i++)
	vecs.push_back(prob->getSolution(i));
      DofWriter::writeFile(vecs, xFilename);
    }
    
    max_psi = max(*getDensity());
  }
  
  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();
    
    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);
      
      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();
    
    DOFVector<double>* psi = getDensity();
    static const int n = 3;
      
    Operator *opPenalty = new Operator(prob->getFeSpace(0));
    addZOT(opPenalty, (H*n*(2.0 - n)) * (signum(valueOf(psi)) - 1.0) * pow<n-1>(valueOf(psi)) );
    prob->addVectorOperator(opPenalty, 0); 
    
    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)) );
    prob->addMatrixOperator(opPenaltyDerivLhs, 0, 1); 
    
    if (mobility_type == 3) { 
      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)
    }
  }
  
  
  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;
  
  std::vector<WorldVector<double> > pos;

  double H;
  double B0, B1;
  double rho1;
  double A;
  double d;
  
  double max_psi;
};


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

  AMDiS::init(argc, argv);
  
#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;
  
  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);
  }
  
  pfcProb.initTimeInterface(); // fill operators and BC
  int error_code = adaptInstat.adapt(); 

  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() );
  
  AMDiS::finalize();
  return error_code;
}