1 /*********************************************************************
2 * dune_heat_algorithm.hpp *
3 * *
4 * This header constains the algorithm and an auxiliar class which *
5 * where original in the heat.cc file *
6 * *
7 * Revision history: *
8 * none *
9 * *
10 * *
11 * Created by Christian Power on 19.06.14. *
12 * Copyright (c) 2014 Christian Power. All rights reserved. *
13 * *
14 *********************************************************************/
15
16 #ifndef DUNE_HEAT_ALGORITHM_HPP
17 #define DUNE_HEAT_ALGORITHM_HPP
18
19 #include <algorithm>
20
21 // local includes
22 #include "dune_typedef_heat.hpp"
23 // #include "/Users/christianpower/cpp/ODE_Solver/implicit_euler.h"
24 #include "/Users/christianpower/cpp/ODE_Solver/bdf.h"
25 #include "w1i_norm.hh"
26
27 // Remember: For every experiment modify
28 // in 'heat.hh'
29 // - RHSfunction
30 // - InitialData
31 // - HeatModel
32 // in 'elliptic.hh'
33 // - EllipticOperator
34 // in 'deformation.hh'
35 // - DeformationCoordFunction
36
37 struct DataOutputParameters:
38 public Dune::Fem::
39 LocalParameter<Dune::Fem::DataOutputParameters, DataOutputParameters> {
40 DataOutputParameters(const std::string name, const int step)
41 : name_(name), step_( step )
42 {}
43 DataOutputParameters(const DataOutputParameters& other)
44 : step_( other.step_ )
45 {}
46 std::string prefix() const {
47 std::stringstream s;
48 // s << "ALE_LiteratureExample-" << step_ << "-";
49 s << name_ << step_ << "-";
50 return s.str();
51 }
52 private:
53 std::string name_;
54 int step_;
55 };
56
57 void BDF::nonlinear_algorithm(){
58 // get time from parameter file
59 double t_0 =
60 Dune::Fem::Parameter::getValue<double>("heat.starttime",0.0);
61 double dT =
62 Dune::Fem::Parameter::getValue<double>("heat.timestep",0.1);
63 double t_end =
64 Dune::Fem::Parameter::getValue<double>("heat.endtime",0.6);
65 const int time_step_no_max = (t_end - t_0)/dT + .1;
66
67 // setting up BDF coefficients
68 const int bdf_no =
69 Dune::Fem::Parameter::getValue<double>("heat.bdf", 1);
70 const std::vector<double> bdf_alpha_coeff { NUMERIK::bdf_alpha_coeff(bdf_no) };
71 const std::vector<double> bdf_gamma_coeff { NUMERIK::bdf_gamma_coeff(bdf_no) };
72 // bdf_*_coeff.back() is the lead coefficient of the polynomial
73
74 // debugging: print BDF coeffs
75 // for(double d : bdf_gamma_coeff)
76 // std::cout << d << ' ';
77 // std::cout << std::endl;
78 // for(double d : bdf_alpha_coeff)
79 // std::cout << d << ' ';
80 // std::cout << '\n' << bdf_alpha_coeff.back() << std::endl;
81
82 // prepare grid from DGF file
83 const std::string gridkey =
84 Dune::Fem::IOInterface::defaultGridKey( GridType::dimension );
85 const std::string gridfile =
86 Dune::Fem::Parameter::getValue< std::string >( gridkey );
87 if( Dune::Fem::MPIManager::rank() == 0 )
88 std::cout << "Loading macro grid: " << gridfile << std::endl;
89 Dune::GridPtr< HostGridType > hostGrid {gridfile};
90 hostGrid ->loadBalance();
91
92 // create grid
93 DeformationCoordFunction deformation {};
94 GridType grid(*hostGrid, deformation);
95 GridPartType gridPart(grid);
96 DiscreteFunctionSpaceType dfSpace(gridPart);
97
98 Dune::Fem::GridTimeProvider<GridType> timeProvider(t_0, grid);
99 deformation.set_time_provider(timeProvider);
100
101 // create FE-functions
102 DiscreteFunctionType U_np1("U_np1", dfSpace); // Numerical solution
103 // In a math book U_np1 == U_n+1.
104 DiscreteFunctionType rhs("rhs", dfSpace); // Right hand side for the LES
105 DiscreteFunctionType load_vector("load_vector", dfSpace);
106 DiscreteFunctionType xi("xi", dfSpace); // For the lineary implicit BDF method
107 DiscreteFunctionType exact_solution("exact_solution", dfSpace);
108 DiscreteFunctionType err_vec("dof_U_np1_minus_exact_solution", dfSpace);
109 // Holds U_np1 - exact_solution
110 std::vector<DiscreteFunctionType> prev_steps_U_nmk; // All previous U_{n-k}
111 for(int i = 0; i < bdf_no; ++i)
112 prev_steps_U_nmk.push_back( DiscreteFunctionType
113 {"U_nm" + std::to_string(bdf_no - (i+1)), dfSpace} );
114 std::vector<DiscreteFunctionType> prev_steps_M_U_nmk; // All previous M_U_{n-k}
115 for(int i = 0; i < bdf_no; ++i)
116 prev_steps_M_U_nmk.push_back( DiscreteFunctionType
117 {"M_U_nm" + std::to_string(bdf_no - (i+1)), dfSpace});
118
119 // stupid check
120 if(prev_steps_U_nmk.size() != prev_steps_M_U_nmk.size())
121 throw std::runtime_error("ERROR in bdf_cycle().");
122
123 // File output
124 L2NormType l2norm(gridPart);
125 H1NormType h1norm(gridPart);
126 std::ofstream l2h1error_ofs
127 {Dune::Fem::Parameter::getValue<std::string>("fem.io.errorFile",
128 "../output/l2h1error"),
129 std::ios_base::app};
130 // Paraview output
131 IOTupleType ioTuple(&err_vec);
132 const int step = 0; // is there any resonable choice, whitout adaptivity?
133 DataOutputType
134 dataOutput(grid, ioTuple,
135 DataOutputParameters(Dune::Fem::Parameter::getValue<std::string>
136 ("fem.io.outputName",
137 "../output/ALE_LiteratureExample-"),
138 step) );
139 // helper function for 'err_vec'
140 auto calc_err_vec = [&U_np1, &exact_solution, &err_vec]{
141 err_vec.assign(U_np1);
142 err_vec.axpy((-1.), exact_solution);
143 };
144
145 InitialDataType initialData {timeProvider};
146 RHSFunctionType f {timeProvider};
147 NonlinearModel model {timeProvider};
148 NonlinearOperator ellipticOp {xi, model, bdf_alpha_coeff.back()};
149 const double solverEps =
150 Dune::Fem::Parameter::getValue<double>("heat.solvereps", 1e-8);
151 LinearInverseOperatorType solver(ellipticOp, solverEps, solverEps); // CG-solver
152
153 auto write_error = [&](std::ostream& os)
154 // helper function for 'l2h1error_ofs'
155 {
156 calc_err_vec();
157 os << std::defaultfloat << timeProvider.deltaT() << ' '
158 << std::scientific
159 << linfty_norm(err_vec) << ' '
160 << w1infty_norm(err_vec)
161 << std::endl;
162
163 // this works!
164 // calc_err_vec();
165 // os << std::defaultfloat << timeProvider.deltaT() << ' '
166 // << std::scientific
167 // << l2norm.norm(err_vec) << ' '
168 // << h1norm.norm(err_vec) << std::endl;
169
170 // 2014linfty version
171 // os << std::defaultfloat << timeProvider.deltaT() << ' '
172 // << std::scientific
173 // << l2norm.distance(exact_solution, U_np1) << ' '
174 // << h1norm.distance(exact_solution, U_np1) << std::endl;
175 };
176
177 // initial steps
178 timeProvider.init(dT); // Do your first action before you enter the for loop.
179 int time_step_no = 0;
180
181 // debugging
182 // std::cout << "\n==========" << std::endl;
183
184 // initialize 'prev_steps_U_nmk' and 'prev_steps_M_U_mmk'
185 for(;
186 time_step_no < std::min(prev_steps_U_nmk.size(), size_t(time_step_no_max));
187 ++time_step_no, timeProvider.next(dT)){
188
189 // debugging
190 // std::cout << std::defaultfloat << "time = " << timeProvider.time()
191 // << std::scientific << std::endl;
192
193 InterpolationType::interpolateFunction(initialData, exact_solution);
194 prev_steps_U_nmk.at(time_step_no).assign(exact_solution);
195 ellipticOp.mass_matrix(exact_solution, prev_steps_M_U_nmk.at(time_step_no));
196
197 U_np1.assign(exact_solution);
198
199 // output
200 // calc_err_vec();
201 // dataOutput.write(timeProvider);
202 write_error(l2h1error_ofs);
203
204 // debugging
205 // std::cout << "exact_solution = " << *exact_solution.dbegin() << '\n'
206 // << "U_np1 = " << *(U_np1.dbegin()) << '\n'
207 // << "prev_steps_M_U_nmk.at(time_step_no) = "
208 // << *prev_steps_M_U_nmk.at(time_step_no).dbegin()
209 // << std::endl;
210 }
211
212 // debugging
213 // std::cout << "\n==========" << std::endl;
214
215 for(;
216 time_step_no <= time_step_no_max;
217 timeProvider.next(dT), ++time_step_no){
218
219 // debugging
220 // std::cout << std::defaultfloat << "time = " << timeProvider.time()
221 // << std::scientific << std::endl;
222
223 // 'rhs', i.e.
224 // calculating: rhs = (-1) · (∑ᵢ₌₁ᵏ αᵢ₋₁ (Mu)ⁿ⁻ᵏ⁺ⁱ)
225 // xi = ∑ᵢ₌₁ᵏ γᵢ₋₁ uⁿ⁻ᵏ⁺ⁱ
226 // NOTE for implicit euler: rhs = Mⁿ uⁿ and xi = uⁿ
227 rhs.clear(); // set dof to zero
228 // xi.clear();
229 for(size_t i=0; i < prev_steps_U_nmk.size(); ++i){
230 rhs.axpy(bdf_alpha_coeff.at(i), prev_steps_M_U_nmk.at(i));
231 // xi.axpy(bdf_gamma_coeff.at(i), prev_steps_U_nmk.at(i));
232 // std::cout << "rhs.axpy(bdf_alpha_coeff.at(i), "
233 // "prev_steps_M_U_nmk.at(i)) = " << *rhs.dbegin() << std::endl;
234 }
235 // std::cout << "xi = " << *xi.dbegin() << std::endl;
236 rhs *= (-1);
237 // std::cout << "rhs *= (-1) : " << *rhs.dbegin() << std::endl;
238 assembleRHS(f, load_vector);
239 // std::cout << "assembleRHS(f, load_vector) = "
240 // << *load_vector.dbegin() << std::endl;
241 rhs.axpy(timeProvider.deltaT(), load_vector);
242 // std::cout << "rhs.axpy(timeProvider.deltaT(), load_vector) = "
243 // << *rhs.dbegin() << std::endl;
244
245 // 'solver'
246 solver(rhs, U_np1);
247 // std::cout << "solver(rhs, U_np1) = " << *U_np1.dbegin() << std::endl;
248
249 // cycle bdf values
250 for(size_t i=0; i < prev_steps_U_nmk.size() - 1; ++i){
251 prev_steps_U_nmk.at(i).assign(prev_steps_U_nmk.at(i+1));
252 prev_steps_M_U_nmk.at(i).assign(prev_steps_M_U_nmk.at(i+1));
253 }
254 prev_steps_U_nmk.back().assign(U_np1);
255 ellipticOp.mass_matrix(U_np1, prev_steps_M_U_nmk.back());
256
257 // output
258 InterpolationType::interpolateFunction(initialData, exact_solution);
259 write_error(l2h1error_ofs);
260 // calc_err_vec();
261 // dataOutput.write(timeProvider);
262 // if(time_step_no == time_step_no_max){
263 // std::cout << std::defaultfloat
264 // << "Time = " << timeProvider.time()
265 // << std::scientific << std::endl;
266 // InterpolationType::interpolateFunction(initialData, exact_solution);
267 // write_error(l2h1error_ofs);
268 // write_error(std::cout);
269 // }
270
271 // debugging
272 // std::cout << "U_np1 = " << *U_np1.dbegin()
273 // << "\n----------" << std::endl;
274 }
275 l2h1error_ofs.close();
276 }
277
278 // // solution_driven paper
279 // void nonlinear_evolution(){ //mcf model
280 // // get time from parameter file
281 // double t_0 =
282 // Dune::Fem::Parameter::getValue<double>("heat.starttime",0.0);
283 // double dT =
284 // Dune::Fem::Parameter::getValue<double>("heat.timestep",0.1);
285 // double t_end =
286 // Dune::Fem::Parameter::getValue<double>("heat.endtime",0.6);
287 // const int time_step_no_max = (t_end - t_0)/dT + .1;
288 //
289 // // prepare grid from DGF file
290 // const std::string gridkey =
291 // Dune::Fem::IOInterface::defaultGridKey( GridType::dimension );
292 // const std::string gridfile =
293 // Dune::Fem::Parameter::getValue< std::string >( gridkey );
294 // if( Dune::Fem::MPIManager::rank() == 0 )
295 // std::cout << "Loading macro grid: " << gridfile << std::endl;
296 //
297 // Dune::GridPtr<HostGridType> hostGrid {gridfile};
298 // hostGrid->loadBalance();
299 //
300 // // create grid
301 // BDF::EvoMapType evoMap {gridfile};
302 // evoMap.save_original_vertices();
303 // DeformationCoordFunction deformation {evoMap};
304 // GridType grid(*hostGrid, deformation);
305 // GridPartType gridPart(grid);
306 // DiscreteFunctionSpaceType dfSpace(gridPart);
307 // Vec_FE_Space r3dfSpace(gridPart);
308 //
309 // Dune::Fem::GridTimeProvider<GridType> timeProvider(t_0, grid);
310 // deformation.set_time_provider(timeProvider);
311 //
312 // // create vec_FE-functions
313 // Vec_FE_Fun next_surface {"next_surface", r3dfSpace}; // surface for next step
314 // Vec_FE_Fun surface_cg_rhs {"surface_cg_rhs", r3dfSpace};
315 //
316 // // create FE-functions
317 // DiscreteFunctionType U_np1 {"U_np1", dfSpace}; // Numerical solution
318 //
319 //
320 // // Abstract functions
321 // R3Identity r3identity {timeProvider};
322 //
323 // // Surface FEM operator
324 // const double alpha =
325 // Dune::Fem::Parameter::
326 // getValue<double>("heat.surface.mass_matrix.regularisation_parameter",0.);
327 // Vec_Model surface_model {timeProvider, alpha};
328 // Vec_Operator surface_operator {surface_model};
329 // const double solverEps =
330 // Dune::Fem::Parameter::getValue<double>("heat.solvereps", 1e-8);
331 // Vec_CG_Solver surface_solver {surface_operator, solverEps, solverEps};
332 //
333 //
334 // // File output
335 // // L2NormType l2norm(gridPart);
336 // // H1NormType h1norm(gridPart);
337 // // std::ofstream l2h1error_ofs {Dune::Fem::
338 // // Parameter::getValue<std::string>("fem.io.errorFile", "../output/l2h1error"),
339 // // std::ios_base::app};
340 //
341 // // Paraview output
342 // IOTupleType ioTuple(&U_np1);
343 // const int step = 0; // is there any resonable choice, whitout adaptivity?
344 // DataOutputType
345 // dataOutput(grid, ioTuple,
346 // DataOutputParameters(Dune::Fem::Parameter::getValue<std::string>
347 // ("fem.io.outputName",
348 // "../output/ALE_LiteratureExample-"),
349 // step) );
350 // // helper function for 'err_vec'
351 // // auto calc_err_vec = [&U_np1, &exact_solution, &err_vec]{
352 // // err_vec.assign(U_np1);
353 // // err_vec.axpy((-1.), exact_solution);
354 // // };
355 //
356 // timeProvider.init(dT); // Do your first action before you enter the for loop.
357 // int time_step_no = 0;
358 //
359 // R3Interpolation::interpolateFunction(r3identity, next_surface);
360 // // next_surface = u_0
361 //
362 // surface_operator.reg_mass_matrix(next_surface, surface_cg_rhs);
363 // // surface_cg_rhs = M_0 u_0
364 //
365 // for(;
366 // time_step_no <= time_step_no_max;
367 // timeProvider.next(dT), ++time_step_no){
368 //
369 // surface_solver(surface_cg_rhs, next_surface);
370 //
371 // evoMap.evolve(next_surface);
372 //
373 // surface_operator.reg_mass_matrix(next_surface, surface_cg_rhs);
374 //
375 // dataOutput.write(timeProvider);
376 // }
377 // }
378
379 void tumor_growth_code(){
380 // get time from parameter file
381 double t_0 =
382 Dune::Fem::Parameter::getValue<double>("heat.starttime",0.0);
383 double dT =
384 Dune::Fem::Parameter::getValue<double>("heat.timestep",0.1);
385 double t_end =
386 Dune::Fem::Parameter::getValue<double>("heat.endtime",0.6);
387 const int time_step_no_max = (t_end - t_0)/dT + .1;
388
389 // prepare grid from DGF file
390 const std::string gridkey =
391 Dune::Fem::IOInterface::defaultGridKey( GridType::dimension );
392 const std::string gridfile =
393 Dune::Fem::Parameter::getValue< std::string >( gridkey );
394 if( Dune::Fem::MPIManager::rank() == 0 )
395 std::cout << "Loading macro grid: " << gridfile << std::endl;
396
397 Dune::GridPtr<HostGridType> hostGrid {gridfile};
398 hostGrid->loadBalance();
399
400 // create grid
401 BDF::EvoMapType evoMap {gridfile};
402 evoMap.save_original_vertices();
403 DeformationCoordFunction deformation {evoMap};
404 GridType grid (*hostGrid, deformation);
405 GridPartType gridPart (grid);
406 DiscreteFunctionSpaceType dfSpace (gridPart);
407 Vec_FE_Space r3dfSpace (gridPart);
408
409 Dune::Fem::GridTimeProvider<GridType> timeProvider(t_0, grid);
410 deformation.set_time_provider(timeProvider);
411
412 // create vec_FE-functions
413 Vec_FE_Fun next_surface {"next_surface", r3dfSpace}; // surface for next step
414 Vec_FE_Fun surface_cg_rhs {"surface_cg_rhs", r3dfSpace};
415
416 // create FE-functions
417 DiscreteFunctionType u_n {"u_n", dfSpace}; // Numerical solution
418 DiscreteFunctionType u_approx {"u_approx", dfSpace}; // Approximaion to u_n
419 DiscreteFunctionType u_cg_rhs {"u_cg_rhs", dfSpace};
420 DiscreteFunctionType w_n {"w_n", dfSpace}; // Numerical solution
421 DiscreteFunctionType w_approx {"w_approx", dfSpace}; // Approximaion to u_n
422 DiscreteFunctionType w_cg_rhs {"w_cg_rhs", dfSpace};
423 DiscreteFunctionType tmp_fef {"tmp_fef", dfSpace}; // tmp finite element function
424
425 // Abstract functions
426 const double mean_value_u
427 = Dune::Fem::Parameter::getValue<double>("tumor_growth.heat.init"
428 "_data.mean_value.u", 1.);
429 const double variance_u
430 = Dune::Fem::Parameter::getValue<double>("tumor_growth.heat.init"
431 "_data.variance.u", .5);
432 const double mean_value_w
433 = Dune::Fem::Parameter::getValue<double>("tumor_growth.heat.init"
434 "_data.mean_value.w", .9);
435 const double variance_w
436 = Dune::Fem::Parameter::getValue<double>("tumor_growth.heat.init"
437 "_data.variance.w", .5);
438
439 const Initial_data_u initial_data_u {mean_value_u, variance_u};
440 const Initial_data_w initial_data_w {mean_value_w, variance_w};
441 R3Identity r3identity {timeProvider};
442
443 // Tumor growth model & operator
444 // Surface model
445 const double alpha
446 = Dune::Fem::Parameter::getValue<double>("tumor_growth.surface.alpha", 1e-3);
447 const double epsilon
448 = Dune::Fem::Parameter::getValue<double>("tumor_growth.surface.epsilon", 1e-2);
449 const double delta
450 = Dune::Fem::Parameter::getValue<double>("tumor_growth.surface.delta", .4);
451
452 const Tumor_surface_model surface_model {timeProvider, alpha, epsilon, delta};
453
454 // Surface operator
455 const double solverEps
456 = Dune::Fem::Parameter::getValue<double>("heat.solvereps", 1e-8);
457
458 const Tumor_surface_operator surface_operator {surface_model, u_approx};
459 Vec_CG_Solver surface_solver {surface_operator, solverEps, solverEps};
460
461 // Brusselator model
462 const double a
463 = Dune::Fem::Parameter::getValue<double>("tumor_growth.heat.a", .1);
464 const double b
465 = Dune::Fem::Parameter::getValue<double>("tumor_growth.heat.b", .9);
466 const double dc
467 = Dune::Fem::Parameter::getValue<double>("tumor_growth.heat.Dc", 10.);
468 const double gamma
469 = Dune::Fem::Parameter::getValue<double>("tumor_growth.heat.gamma", 100.);
470 using Growth = Tumor_growth::Growth;
471
472 const Tumor_Brusselator_model u_model {dT, gamma, a, Growth::promoting};
473 const Tumor_Brusselator_model w_model {dT, gamma, b, Growth::inhibiting, dc};
474
475 // Brusselator operators
476 const Tumor_Brusselator_operator u_operator {u_model, u_approx, w_approx};
477 const Tumor_Brusselator_operator w_operator {w_model, u_n, u_n};
478 LinearInverseOperatorType u_solver {u_operator, solverEps, solverEps};
479 LinearInverseOperatorType w_solver {w_operator, solverEps, solverEps};
480
481 // Norms and file output
482 L2NormType l2norm(gridPart);
483 H1NormType h1norm(gridPart);
484 // std::ofstream l2h1error_ofs {Dune::Fem::
485 // Parameter::getValue<std::string>("fem.io.errorFile", "../output/l2h1error"),
486 // std::ios_base::app};
487
488 // Paraview output
489 IOTupleType u_ioTuple(&u_n);
490 IOTupleType w_ioTuple(&w_n);
491 const int step = 0; // is there any resonable choice, whitout adaptivity?
492 DataOutputType
493 u_dataOutput(grid, u_ioTuple,
494 DataOutputParameters(Dune::Fem::Parameter::getValue<std::string>
495 ("tumor_growth.io.u.output_name",
496 "../output/ALE_LiteratureExample-"),
497 step) );
498 DataOutputType
499 w_dataOutput(grid, w_ioTuple,
500 DataOutputParameters(Dune::Fem::Parameter::getValue<std::string>
501 ("tumor_growth.io.w.output_name",
502 "../output/ALE_LiteratureExample2-"),
503 step) );
504 auto data_writer = [](DataOutputType& u_data, DataOutputType& w_data,
505 const Dune::Fem::GridTimeProvider<GridType>& tp){
506 u_data.write(tp);
507 w_data.write(tp);
508 };
509 // helper function for 'err_vec'
510 // auto calc_err_vec = [&u_np1, &exact_solution, &err_vec]{
511 // err_vec.assign(u_np1);
512 // err_vec.axpy((-1.), exact_solution);
513 // };
514
515 // starting calculation
516 timeProvider.init(dT); // Do your first action before you enter the for-loop.
517 int time_step_no = 0;
518
519 // evoMap.save_original_vertices(); already done
520 InterpolationType::interpolateFunction(initial_data_u, u_n);
521 InterpolationType::interpolateFunction(initial_data_w, w_n);
522
523 u_approx.assign(u_n);
524 w_approx.assign(w_n);
525
526 // This space is useful for BDF codes
527
528 // dataOutput.write(timeProvider);
529 data_writer(u_dataOutput, w_dataOutput, timeProvider);
530
531 for(;
532 time_step_no <= time_step_no_max;
533 timeProvider.next(dT), ++time_step_no){
534
535 R3Interpolation::interpolateFunction(r3identity, next_surface);
536
537 surface_operator.generate_rhs(next_surface, surface_cg_rhs);
538
539 u_operator.generate_rhs_old_surface(u_n, u_cg_rhs);
540 w_operator.generate_rhs_old_surface(w_n, w_cg_rhs);
541
542 surface_solver(surface_cg_rhs, next_surface);
543
544 evoMap.evolve(next_surface);
545
546 u_operator.generate_rhs_new_surface(tmp_fef);
547 u_cg_rhs.axpy(1., tmp_fef);
548 w_operator.generate_rhs_new_surface(tmp_fef); // this can be done faster
549 w_cg_rhs.axpy(1., tmp_fef);
550
551 u_solver(u_cg_rhs, u_n);
552 w_solver(w_cg_rhs, w_n);
553
554 u_approx.assign(u_n);
555 w_approx.assign(w_n);
556
557 // dataOutput.write(timeProvider);
558 data_writer(u_dataOutput, w_dataOutput, timeProvider);
559 }
560 }
561
562 // // begin paperALE experiment
563 //
564 // // This has to be modified for the ODE solver
565 // struct EigenFullPivLuSolv{
566 // Eigen::Vector3d linear_solve(const Eigen::Matrix3d& A,
567 // const Eigen::Vector3d& b) const {
568 // return A.fullPivLu().solve(b);
569 // }
570 // };
571 // struct EigenNorm{
572 // double norm(const Eigen::Vector3d& v) const {
573 // return v.norm();
574 // }
575 // };
576 // struct SurfaceVectorfield{
577 // void set_time(const double time) { t = time; }
578 // double get_time() { return t; }
579 // Eigen::Vector3d evaluate(const Eigen::Vector3d vec) const {
580 // double x = vec(0), y = vec(1), z = vec(2);
581 // return Eigen::Vector3d {
582 // 16.*M_PI*pow(x,3)*cos(2*M_PI*t)/( (pow(y,2) + pow(z,2) + 16.*pow(x,2)/pow(sin(2.*M_PI*t) + 4.,2) )*pow(sin(2.*M_PI*t) + 4.,3) ),
583 // 4.*M_PI*pow(x,2)*y*cos(2*M_PI*t)/((pow(y,2) + pow(z,2) + 16.*pow(x,2)/pow(sin(2.*M_PI*t) + 4.,2) )*pow(sin(2.*M_PI*t) + 4.,2)),
584 // 4.*M_PI*pow(x,2)*z*cos(2.*M_PI*t)/((pow(y,2) + pow(z,2) + 16.*pow(x,2)/pow(sin(2*M_PI*t) + 4.,2) )*pow(sin(2.*M_PI*t) + 4,2) ) };
585 // }
586 // Eigen::Matrix3d jacobian(const Eigen::Vector3d vec) const {
587 // double x = vec(0), y = vec(1), z = vec(2);
588 // Eigen::Matrix3d jac_f_t_vec;
589 // jac_f_t_vec <<
590 // // first row
591 // 48*M_PI*pow(x,2)*cos(2*M_PI*t)/((pow(y,2) + pow(z,2) + 16*pow(x,2)/pow(sin(2*M_PI*t) + 4,2))*pow(sin(2*M_PI*t) + 4,3)) - 512*M_PI*pow(x,4)*cos(2*M_PI*t)/(pow(pow(y,2) + pow(z,2) + 16*pow(x,2)/pow(sin(2*M_PI*t) + 4,2),2)*pow(sin(2*M_PI*t) + 4,5)),
592 // -32*M_PI*pow(x,3)*y*cos(2*M_PI*t)/(pow(pow(y,2) + pow(z,2) + 16*pow(x,2)/pow(sin(2*M_PI*t) + 4,2),2)*pow(sin(2*M_PI*t) + 4,3)),
593 // -32*M_PI*pow(x,3)*z*cos(2*M_PI*t)/(pow(pow(y,2) + pow(z,2) + 16*pow(x,2)/pow(sin(2*M_PI*t) + 4,2),2)*pow(sin(2*M_PI*t) + 4,3)),
594 // // second row
595 // 8*M_PI*x*y*cos(2*M_PI*t)/((pow(y,2) + pow(z,2) + 16*pow(x,2)/pow(sin(2*M_PI*t) + 4,2))*pow(sin(2*M_PI*t) + 4,2)) - 128*M_PI*pow(x,3)*y*cos(2*M_PI*t)/(pow(pow(y,2) + pow(z,2) + 16*pow(x,2)/pow(sin(2*M_PI*t) + 4,2),2)*pow(sin(2*M_PI*t) + 4,4)),
596 // -8*M_PI*pow(x,2)*pow(y,2)*cos(2*M_PI*t)/(pow(pow(y,2) + pow(z,2) + 16*pow(x,2)/pow(sin(2*M_PI*t) + 4,2),2)*pow(sin(2*M_PI*t) + 4,2)) + 4*M_PI*pow(x,2)*cos(2*M_PI*t)/((pow(y,2) + pow(z,2) + 16*pow(x,2)/pow(sin(2*M_PI*t) + 4,2))*pow(sin(2*M_PI*t) + 4,2)),
597 // -8*M_PI*pow(x,2)*y*z*cos(2*M_PI*t)/(pow(pow(y,2) + pow(z,2) + 16*pow(x,2)/pow(sin(2*M_PI*t) + 4,2),2)*pow(sin(2*M_PI*t) + 4,2)),
598 // // third row
599 // 8*M_PI*x*z*cos(2*M_PI*t)/((pow(y,2) + pow(z,2) + 16*pow(x,2)/pow(sin(2*M_PI*t) + 4,2))*pow(sin(2*M_PI*t) + 4,2)) - 128*M_PI*pow(x,3)*z*cos(2*M_PI*t)/(pow(pow(y,2) + pow(z,2) + 16*pow(x,2)/pow(sin(2*M_PI*t) + 4,2),2)*pow(sin(2*M_PI*t) + 4,4)),
600 // -8*M_PI*pow(x,2)*y*z*cos(2*M_PI*t)/(pow(pow(y,2) + pow(z,2) + 16*pow(x,2)/pow(sin(2*M_PI*t) + 4,2),2)*pow(sin(2*M_PI*t) + 4,2)),
601 // -8*M_PI*pow(x,2)*pow(z,2)*cos(2*M_PI*t)/(pow(pow(y,2) + pow(z,2) + 16*pow(x,2)/pow(sin(2*M_PI*t) + 4,2),2)*pow(sin(2*M_PI*t) + 4,2)) + 4*M_PI*pow(x,2)*cos(2*M_PI*t)/((pow(y,2) + pow(z,2) + 16*pow(x,2)/pow(sin(2*M_PI*t) + 4,2))*pow(sin(2*M_PI*t) + 4,2));
602 // return jac_f_t_vec;
603 // }
604 // private:
605 // double t;
606 // };
607 // struct EigenIdMatrix3d{
608 // Eigen::Matrix3d identity_matrix() const { return Eigen::Matrix3d::Identity(); }
609 // };
610 //
611 //
612 // typedef NUMERIK::ImplicitEuler<EigenFullPivLuSolv, EigenNorm, SurfaceVectorfield,
613 // EigenIdMatrix3d, Eigen::Matrix3d, Eigen::Vector3d>
614 // EvoMapIntegrator;
615 //
616 // void BDF::algorithm()
617 // {
618 // /**********************************************************************
619 // * Pseudocode:
620 // * u⁰ = interpol(solution) >> plot
621 // * container = M⁰u⁰ >> save
622 // * for-loop n = 0 ↦ N with N·Δt + t_begin = t_end
623 // * t += ∆τ (set time)
624 // * tmp = container == Mⁿuⁿ
625 // * rhs = fⁿ⁺¹ == load Vector
626 // * tmp += Δτ · rhs (solution == Mⁿuⁿ + Δτ·rhs)
627 // * uⁿ⁺¹ == solution = (Mⁿ⁺¹ + Δτ Aⁿ⁺¹)⁻¹ tmp >> save
628 // * container = Mⁿ⁺¹ uⁿ⁺¹
629 // * end-for-loop
630 // *
631 // **********************************************************************/
632 //
633 // // create grid from DGF file
634 // const std::string gridkey =
635 // Dune::Fem::IOInterface::defaultGridKey( GridType::dimension );
636 // const std::string gridfile =
637 // Dune::Fem::Parameter::getValue< std::string >( gridkey );
638 // if( Dune::Fem::MPIManager::rank() == 0 )
639 // std::cout << "Loading macro grid: " << gridfile << std::endl;
640 //
641 // Dune::GridPtr< HostGridType > hostGrid {gridfile};
642 // hostGrid ->loadBalance();
643 //
644 // BDF::EvoMapType evoMap {gridfile};
645 // DeformationCoordFunction deformation {evoMap};
646 // // constructs an unsorted map with domain and range equal the vertexes from the
647 // // initial grid
648 //
649 // GridType grid(*hostGrid, deformation);
650 // GridPartType gridPart(grid);
651 // DiscreteFunctionSpaceType dfSpace(gridPart);
652 //
653 // L2NormType l2norm(gridPart);
654 // H1NormType h1norm(gridPart);
655 // std::ofstream l2h1error_ofs
656 // {Dune::Fem::Parameter::getValue<std::string>("fem.io.errorFile",
657 // "../output/l2h1error")};
658 // //std::ios_base::app};
659 // //l2h1error_ofs << std::scientific;
660 // //l2h1error_ofs << "L2Error\tH1Error" << std::scientific << std::endl;
661 // //std::ostringstream l2h1error_helper_oss;
662 // //l2h1error_helper_oss << "time\ttimestep" << std::endl;
663 //
664 // double t_0 =
665 // Dune::Fem::Parameter::getValue<double>("heat.starttime",0.0);
666 // double dT =
667 // Dune::Fem::Parameter::getValue<double>("heat.timestep",0.1);
668 // double t_end =
669 // Dune::Fem::Parameter::getValue<double>("heat.endtime",0.6);
670 // std::cout << (t_end - t_0)/dT + .1 << std::endl;
671 // const int itno = (t_end - t_0)/dT + .1;
672 //
673 // Dune::Fem::GridTimeProvider< GridType > timeProvider(t_0, grid);
674 // timeProvider.init(dT); // Do your first action before you enter the for loop.
675 // // std::cout << timeProvider.time() << '\t'
676 // // << timeProvider.deltaT() << "\tstarting action."
677 // // << std::endl;
678 // // l2h1error_helper_oss << timeProvider.time() << '\t'
679 // // << timeProvider.deltaT() << "\tstarting action."
680 // // << std::endl;
681 //
682 // DiscreteFunctionType solutionContainer( "U_nContainer", dfSpace ),
683 // U_n("U_n",dfSpace), rhs( "rhs", dfSpace ), tmp( "tmp", dfSpace ),
684 // exactSolution("exactSolution",dfSpace);
685 // InitialDataType initialData;
686 // InterpolationType::interpolateFunction( initialData, U_n );
687 // InterpolationType::interpolateFunction( initialData, exactSolution );
688 // // solutionContainer is used for higher order BDF methods,
689 // // U_n == uⁿ, rhs == tmp1, tmp == tmp2
690 // // (create discrete functions for intermediate functionals)
691 //
692 // RHSFunctionType f;
693 // // This expression is very long; it's literally the same f as in L(u) = f;
694 // // together with the function/method in rhs.hh,
695 //
696 // IOTupleType ioTuple(&U_n);
697 // const int step = 0; // is there any resonable choise, whitout adaptivity?
698 // DataOutputType
699 // dataOutput(grid, ioTuple,
700 // DataOutputParameters(Dune::Fem::Parameter::getValue<std::string>
701 // ("fem.io.outputName",
702 // "../output/ALE_LiteratureExample-"),
703 // step) );
704 //
705 // auto write_error = [&](std::ostream& os){
706 // os << std::defaultfloat << timeProvider.time() << ' '
707 // << std::scientific
708 // << l2norm.distance(exactSolution, U_n) << ' '
709 // << h1norm.distance(exactSolution, U_n) << std::endl;
710 // };
711 // write_error(l2h1error_ofs);
712 //
713 // EllipticOperatorType
714 // ellipticOp( HeatModelType {timeProvider} );
715 //
716 // const double solverEps =
717 // Dune::Fem::Parameter::getValue< double >( "heat.solvereps", 1e-8 );
718 // LinearInverseOperatorType solver( ellipticOp, solverEps, solverEps );
719 // // Initializing CG-Solver
720 // // CAPG: it's very slow. Can't this be speed up somehow?
721 //
722 // dataOutput.write( timeProvider );
723 //
724 // ellipticOp.mass_matrix(U_n, rhs);
725 // // calculating: rhs = Mⁿ uⁿ (uⁿ == solutionContainer)
726 //
727 // // std::vector<double> timecontainer, l2errorcontainer, h1errorcontainer;
728 // int iter = 0;
729 // for(timeProvider.next(dT); iter < itno; timeProvider.next(dT), ++iter)
730 // // put your loop-action here, but not the last action
731 // {
732 // // for debugging reasons
733 // // std::cout << timeProvider.time() << '\t' << timeProvider.deltaT()
734 // // << "\tfor loop action." << std::endl;
735 // // l2h1error_helper_oss << timeProvider.time() << '\t' << timeProvider.deltaT()
736 // // << "\tfor loop action." << std::endl;
737 //
738 // f.setTime(timeProvider.time());
739 // // chance the time for the (long) RHS function f to tⁿ⁺¹
740 // // ('f' from -∆u + \div(v) u + \matdot{u} = f)
741 //
742 // initialData.setTime(timeProvider.time());
743 //
744 // std::array<double,2> time_array {
745 // {timeProvider.time()- timeProvider.deltaT(), timeProvider.time()}
746 // };
747 // EvoMapIntegrator impl_euler { NUMERIK::NewtonParameters<> {1e-10} };
748 // evoMap.evolve(time_array.begin(), time_array.end(), impl_euler);
749 //
750 // assembleRHS( f, tmp );
751 // // assemly stiffness/load vector; the vector is called 'tmp';
752 // // in the sense of above it is fⁿ⁺¹
753 //
754 // rhs.axpy( timeProvider.deltaT(), tmp );
755 // // rhs += Δt * tmp (i.e. rhs += \Delta t * tmp)
756 // // Just calculated: ( uⁿ + ∆t EllipticOperator(uⁿ)) + ∆t fⁿ
757 //
758 // //InterpolationType::interpolateFunction( initialData, U_n );
759 // solver( rhs, U_n );
760 // // Solve: (id + \Delta t EllipticOperator) u_{n+1} = rhs,
761 // // where rhs = ( u_n + \Delta t EllipticOperator(u_n)) + \Delta t f_n
762 // // CAPG remark: this is very very slow.
763 //
764 // dataOutput.write( timeProvider );
765 //
766 // InterpolationType::interpolateFunction( initialData, exactSolution );
767 // write_error(l2h1error_ofs);
768 //
769 // ellipticOp.mass_matrix(U_n, rhs);
770 // // calculating: rhs = Mⁿ uⁿ (uⁿ == solutionContainer)
771 //
772 // // l2errorcontainer.push_back(l2norm.distance(exactSolution, U_n));
773 // // h1errorcontainer.push_back(h1norm.distance(exactSolution, U_n));
774 // // timecontainer.push_back(timeProvider.time());
775 // }
776 // // l2h1error_ofs << "#\n" << l2h1error_helper_oss.str() << '#' << std::endl;
777 //
778 // // l2h1error_ofs << timeProvider.deltaT() << ' '
779 // // << timecontainer.back() << ' '
780 // // << l2errorcontainer.back() << ' '
781 // // << h1errorcontainer.back() << std::endl;
782 // l2h1error_ofs.close();
783 // }
784 // // end paperALE experiment
785 #endif // #ifndef DUNE_HEAT_ALGORITHM_HPP