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 // begin 2014nonlinear experiment
58 void BDF::nonlinear_algorithm(){
59 // get time from parameter file
60 double t_0 =
61 Dune::Fem::Parameter::getValue<double>("heat.starttime",0.0);
62 double dT =
63 Dune::Fem::Parameter::getValue<double>("heat.timestep",0.1);
64 double t_end =
65 Dune::Fem::Parameter::getValue<double>("heat.endtime",0.6);
66 const int time_step_no_max = (t_end - t_0)/dT + .1;
67
68 // setting up BDF coefficients
69 const int bdf_no =
70 Dune::Fem::Parameter::getValue<double>("heat.bdf", 1);
71 const std::vector<double> bdf_alpha_coeff { NUMERIK::bdf_alpha_coeff(bdf_no) };
72 const std::vector<double> bdf_gamma_coeff { NUMERIK::bdf_gamma_coeff(bdf_no) };
73 // bdf_*_coeff.back() is the lead coefficient of the polynomial
74
75 // debugging: print BDF coeffs
76 // for(double d : bdf_gamma_coeff)
77 // std::cout << d << ' ';
78 // std::cout << std::endl;
79 // for(double d : bdf_alpha_coeff)
80 // std::cout << d << ' ';
81 // std::cout << '\n' << bdf_alpha_coeff.back() << std::endl;
82
83 // prepare grid from DGF file
84 const std::string gridkey =
85 Dune::Fem::IOInterface::defaultGridKey( GridType::dimension );
86 const std::string gridfile =
87 Dune::Fem::Parameter::getValue< std::string >( gridkey );
88 if( Dune::Fem::MPIManager::rank() == 0 )
89 std::cout << "Loading macro grid: " << gridfile << std::endl;
90 Dune::GridPtr< HostGridType > hostGrid {gridfile};
91 hostGrid ->loadBalance();
92
93 // create grid
94 DeformationCoordFunction deformation {};
95 GridType grid(*hostGrid, deformation);
96 GridPartType gridPart(grid);
97 DiscreteFunctionSpaceType dfSpace(gridPart);
98
99 Dune::Fem::GridTimeProvider<GridType> timeProvider(t_0, grid);
100 deformation.set_time_provider(timeProvider);
101
102 // create FE-functions
103 DiscreteFunctionType U_np1("U_np1", dfSpace); // Numerical solution
104 // In a math book U_np1 == U_n+1.
105 DiscreteFunctionType rhs("rhs", dfSpace); // Right hand side for the LES
106 DiscreteFunctionType load_vector("load_vector", dfSpace);
107 DiscreteFunctionType xi("xi", dfSpace); // For the lineary implicit BDF method
108 DiscreteFunctionType exact_solution("exact_solution", dfSpace);
109 DiscreteFunctionType err_vec("dof_U_np1_minus_exact_solution", dfSpace);
110 // Holds U_np1 - exact_solution
111 std::vector<DiscreteFunctionType> prev_steps_U_nmk; // All previous U_{n-k}
112 for(int i = 0; i < bdf_no; ++i)
113 prev_steps_U_nmk.push_back( DiscreteFunctionType
114 {"U_nm" + std::to_string(bdf_no - (i+1)), dfSpace} );
115 std::vector<DiscreteFunctionType> prev_steps_M_U_nmk; // All previous M_U_{n-k}
116 for(int i = 0; i < bdf_no; ++i)
117 prev_steps_M_U_nmk.push_back( DiscreteFunctionType
118 {"M_U_nm" + std::to_string(bdf_no - (i+1)), dfSpace});
119
120 // stupid check
121 if(prev_steps_U_nmk.size() != prev_steps_M_U_nmk.size())
122 throw std::runtime_error("ERROR in bdf_cycle().");
123
124 // File output
125 L2NormType l2norm(gridPart);
126 H1NormType h1norm(gridPart);
127 std::ofstream l2h1error_ofs
128 {Dune::Fem::Parameter::getValue<std::string>("fem.io.errorFile",
129 "../output/l2h1error"),
130 std::ios_base::app};
131 // Paraview output
132 IOTupleType ioTuple(&err_vec);
133 const int step = 0; // is there any resonable choice, whitout adaptivity?
134 DataOutputType
135 dataOutput(grid, ioTuple,
136 DataOutputParameters(Dune::Fem::Parameter::getValue<std::string>
137 ("fem.io.outputName",
138 "../output/ALE_LiteratureExample-"),
139 step) );
140 // helper function for 'err_vec'
141 auto calc_err_vec = [&U_np1, &exact_solution, &err_vec]{
142 err_vec.assign(U_np1);
143 err_vec.axpy((-1.), exact_solution);
144 };
145
146 InitialDataType initialData {timeProvider};
147 RHSFunctionType f {timeProvider};
148 NonlinearModel model {timeProvider};
149 NonlinearOperator ellipticOp {xi, model, bdf_alpha_coeff.back()};
150 const double solverEps =
151 Dune::Fem::Parameter::getValue<double>("heat.solvereps", 1e-8);
152 LinearInverseOperatorType solver(ellipticOp, solverEps, solverEps); // CG-solver
153
154 auto write_error = [&](std::ostream& os)
155 // helper function for 'l2h1error_ofs'
156 {
157 calc_err_vec();
158 os << std::defaultfloat << timeProvider.deltaT() << ' '
159 << std::scientific
160 << linfty_norm(err_vec) << ' '
161 << w1infty_norm(err_vec)
162 << std::endl;
163
164 // this works!
165 // calc_err_vec();
166 // os << std::defaultfloat << timeProvider.deltaT() << ' '
167 // << std::scientific
168 // << l2norm.norm(err_vec) << ' '
169 // << h1norm.norm(err_vec) << std::endl;
170
171 // 2014linfty version
172 // os << std::defaultfloat << timeProvider.deltaT() << ' '
173 // << std::scientific
174 // << l2norm.distance(exact_solution, U_np1) << ' '
175 // << h1norm.distance(exact_solution, U_np1) << std::endl;
176 };
177
178 // initial steps
179 timeProvider.init(dT); // Do your first action before you enter the for loop.
180 int time_step_no = 0;
181
182 // debugging
183 // std::cout << "\n==========" << std::endl;
184
185 // initialize 'prev_steps_U_nmk' and 'prev_steps_M_U_mmk'
186 for(;
187 time_step_no < std::min(prev_steps_U_nmk.size(), size_t(time_step_no_max));
188 ++time_step_no, timeProvider.next(dT)){
189
190 // debugging
191 // std::cout << std::defaultfloat << "time = " << timeProvider.time()
192 // << std::scientific << std::endl;
193
194 InterpolationType::interpolateFunction(initialData, exact_solution);
195 prev_steps_U_nmk.at(time_step_no).assign(exact_solution);
196 ellipticOp.mass_matrix(exact_solution, prev_steps_M_U_nmk.at(time_step_no));
197
198 U_np1.assign(exact_solution);
199
200 // output
201 // calc_err_vec();
202 // dataOutput.write(timeProvider);
203 write_error(l2h1error_ofs);
204
205 // debugging
206 // std::cout << "exact_solution = " << *exact_solution.dbegin() << '\n'
207 // << "U_np1 = " << *(U_np1.dbegin()) << '\n'
208 // << "prev_steps_M_U_nmk.at(time_step_no) = "
209 // << *prev_steps_M_U_nmk.at(time_step_no).dbegin()
210 // << std::endl;
211 }
212
213 // debugging
214 // std::cout << "\n==========" << std::endl;
215
216 for(;
217 time_step_no <= time_step_no_max;
218 timeProvider.next(dT), ++time_step_no){
219
220 // debugging
221 // std::cout << std::defaultfloat << "time = " << timeProvider.time()
222 // << std::scientific << std::endl;
223
224 // 'rhs', i.e.
225 // calculating: rhs = (-1) · (∑ᵢ₌₁ᵏ αᵢ₋₁ (Mu)ⁿ⁻ᵏ⁺ⁱ)
226 // xi = ∑ᵢ₌₁ᵏ γᵢ₋₁ uⁿ⁻ᵏ⁺ⁱ
227 // NOTE for implicit euler: rhs = Mⁿ uⁿ and xi = uⁿ
228 rhs.clear(); // set dof to zero
229 // xi.clear();
230 for(size_t i=0; i < prev_steps_U_nmk.size(); ++i){
231 rhs.axpy(bdf_alpha_coeff.at(i), prev_steps_M_U_nmk.at(i));
232 // xi.axpy(bdf_gamma_coeff.at(i), prev_steps_U_nmk.at(i));
233 // std::cout << "rhs.axpy(bdf_alpha_coeff.at(i), "
234 // "prev_steps_M_U_nmk.at(i)) = " << *rhs.dbegin() << std::endl;
235 }
236 // std::cout << "xi = " << *xi.dbegin() << std::endl;
237 rhs *= (-1);
238 // std::cout << "rhs *= (-1) : " << *rhs.dbegin() << std::endl;
239 assembleRHS(f, load_vector);
240 // std::cout << "assembleRHS(f, load_vector) = "
241 // << *load_vector.dbegin() << std::endl;
242 rhs.axpy(timeProvider.deltaT(), load_vector);
243 // std::cout << "rhs.axpy(timeProvider.deltaT(), load_vector) = "
244 // << *rhs.dbegin() << std::endl;
245
246 // 'solver'
247 solver(rhs, U_np1);
248 // std::cout << "solver(rhs, U_np1) = " << *U_np1.dbegin() << std::endl;
249
250 // cycle bdf values
251 for(size_t i=0; i < prev_steps_U_nmk.size() - 1; ++i){
252 prev_steps_U_nmk.at(i).assign(prev_steps_U_nmk.at(i+1));
253 prev_steps_M_U_nmk.at(i).assign(prev_steps_M_U_nmk.at(i+1));
254 }
255 prev_steps_U_nmk.back().assign(U_np1);
256 ellipticOp.mass_matrix(U_np1, prev_steps_M_U_nmk.back());
257
258 // output
259 InterpolationType::interpolateFunction(initialData, exact_solution);
260 write_error(l2h1error_ofs);
261 // calc_err_vec();
262 // dataOutput.write(timeProvider);
263 // if(time_step_no == time_step_no_max){
264 // std::cout << std::defaultfloat
265 // << "Time = " << timeProvider.time()
266 // << std::scientific << std::endl;
267 // InterpolationType::interpolateFunction(initialData, exact_solution);
268 // write_error(l2h1error_ofs);
269 // write_error(std::cout);
270 // }
271
272 // debugging
273 // std::cout << "U_np1 = " << *U_np1.dbegin()
274 // << "\n----------" << std::endl;
275 }
276 l2h1error_ofs.close();
277 }
278 // end 2014nonlinear experiment
279
280
281 // // begin paperALE experiment
282 //
283 // // This has to be modified for the ODE solver
284 // struct EigenFullPivLuSolv{
285 // Eigen::Vector3d linear_solve(const Eigen::Matrix3d& A,
286 // const Eigen::Vector3d& b) const {
287 // return A.fullPivLu().solve(b);
288 // }
289 // };
290 // struct EigenNorm{
291 // double norm(const Eigen::Vector3d& v) const {
292 // return v.norm();
293 // }
294 // };
295 // struct SurfaceVectorfield{
296 // void set_time(const double time) { t = time; }
297 // double get_time() { return t; }
298 // Eigen::Vector3d evaluate(const Eigen::Vector3d vec) const {
299 // double x = vec(0), y = vec(1), z = vec(2);
300 // return Eigen::Vector3d {
301 // 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) ),
302 // 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)),
303 // 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) ) };
304 // }
305 // Eigen::Matrix3d jacobian(const Eigen::Vector3d vec) const {
306 // double x = vec(0), y = vec(1), z = vec(2);
307 // Eigen::Matrix3d jac_f_t_vec;
308 // jac_f_t_vec <<
309 // // first row
310 // 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)),
311 // -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)),
312 // -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)),
313 // // second row
314 // 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)),
315 // -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)),
316 // -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)),
317 // // third row
318 // 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)),
319 // -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)),
320 // -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));
321 // return jac_f_t_vec;
322 // }
323 // private:
324 // double t;
325 // };
326 // struct EigenIdMatrix3d{
327 // Eigen::Matrix3d identity_matrix() const { return Eigen::Matrix3d::Identity(); }
328 // };
329 //
330 //
331 // typedef NUMERIK::ImplicitEuler<EigenFullPivLuSolv, EigenNorm, SurfaceVectorfield,
332 // EigenIdMatrix3d, Eigen::Matrix3d, Eigen::Vector3d>
333 // EvoMapIntegrator;
334 //
335 // void BDF::algorithm()
336 // {
337 // /**********************************************************************
338 // * Pseudocode:
339 // * u⁰ = interpol(solution) >> plot
340 // * container = M⁰u⁰ >> save
341 // * for-loop n = 0 ↦ N with N·Δt + t_begin = t_end
342 // * t += ∆τ (set time)
343 // * tmp = container == Mⁿuⁿ
344 // * rhs = fⁿ⁺¹ == load Vector
345 // * tmp += Δτ · rhs (solution == Mⁿuⁿ + Δτ·rhs)
346 // * uⁿ⁺¹ == solution = (Mⁿ⁺¹ + Δτ Aⁿ⁺¹)⁻¹ tmp >> save
347 // * container = Mⁿ⁺¹ uⁿ⁺¹
348 // * end-for-loop
349 // *
350 // **********************************************************************/
351 //
352 // // create grid from DGF file
353 // const std::string gridkey =
354 // Dune::Fem::IOInterface::defaultGridKey( GridType::dimension );
355 // const std::string gridfile =
356 // Dune::Fem::Parameter::getValue< std::string >( gridkey );
357 // if( Dune::Fem::MPIManager::rank() == 0 )
358 // std::cout << "Loading macro grid: " << gridfile << std::endl;
359 //
360 // Dune::GridPtr< HostGridType > hostGrid {gridfile};
361 // hostGrid ->loadBalance();
362 //
363 // BDF::EvoMapType evoMap {gridfile};
364 // DeformationCoordFunction deformation {evoMap};
365 // // constructs an unsorted map with domain and range equal the vertexes from the
366 // // initial grid
367 //
368 // GridType grid(*hostGrid, deformation);
369 // GridPartType gridPart(grid);
370 // DiscreteFunctionSpaceType dfSpace(gridPart);
371 //
372 // L2NormType l2norm(gridPart);
373 // H1NormType h1norm(gridPart);
374 // std::ofstream l2h1error_ofs
375 // {Dune::Fem::Parameter::getValue<std::string>("fem.io.errorFile",
376 // "../output/l2h1error")};
377 // //std::ios_base::app};
378 // //l2h1error_ofs << std::scientific;
379 // //l2h1error_ofs << "L2Error\tH1Error" << std::scientific << std::endl;
380 // //std::ostringstream l2h1error_helper_oss;
381 // //l2h1error_helper_oss << "time\ttimestep" << std::endl;
382 //
383 // double t_0 =
384 // Dune::Fem::Parameter::getValue<double>("heat.starttime",0.0);
385 // double dT =
386 // Dune::Fem::Parameter::getValue<double>("heat.timestep",0.1);
387 // double t_end =
388 // Dune::Fem::Parameter::getValue<double>("heat.endtime",0.6);
389 // std::cout << (t_end - t_0)/dT + .1 << std::endl;
390 // const int itno = (t_end - t_0)/dT + .1;
391 //
392 // Dune::Fem::GridTimeProvider< GridType > timeProvider(t_0, grid);
393 // timeProvider.init(dT); // Do your first action before you enter the for loop.
394 // // std::cout << timeProvider.time() << '\t'
395 // // << timeProvider.deltaT() << "\tstarting action."
396 // // << std::endl;
397 // // l2h1error_helper_oss << timeProvider.time() << '\t'
398 // // << timeProvider.deltaT() << "\tstarting action."
399 // // << std::endl;
400 //
401 // DiscreteFunctionType solutionContainer( "U_nContainer", dfSpace ),
402 // U_n("U_n",dfSpace), rhs( "rhs", dfSpace ), tmp( "tmp", dfSpace ),
403 // exactSolution("exactSolution",dfSpace);
404 // InitialDataType initialData;
405 // InterpolationType::interpolateFunction( initialData, U_n );
406 // InterpolationType::interpolateFunction( initialData, exactSolution );
407 // // solutionContainer is used for higher order BDF methods,
408 // // U_n == uⁿ, rhs == tmp1, tmp == tmp2
409 // // (create discrete functions for intermediate functionals)
410 //
411 // RHSFunctionType f;
412 // // This expression is very long; it's literally the same f as in L(u) = f;
413 // // together with the function/method in rhs.hh,
414 //
415 // IOTupleType ioTuple(&U_n);
416 // const int step = 0; // is there any resonable choise, whitout adaptivity?
417 // DataOutputType
418 // dataOutput(grid, ioTuple,
419 // DataOutputParameters(Dune::Fem::Parameter::getValue<std::string>
420 // ("fem.io.outputName",
421 // "../output/ALE_LiteratureExample-"),
422 // step) );
423 //
424 // auto write_error = [&](std::ostream& os){
425 // os << std::defaultfloat << timeProvider.time() << ' '
426 // << std::scientific
427 // << l2norm.distance(exactSolution, U_n) << ' '
428 // << h1norm.distance(exactSolution, U_n) << std::endl;
429 // };
430 // write_error(l2h1error_ofs);
431 //
432 // EllipticOperatorType
433 // ellipticOp( HeatModelType {timeProvider} );
434 //
435 // const double solverEps =
436 // Dune::Fem::Parameter::getValue< double >( "heat.solvereps", 1e-8 );
437 // LinearInverseOperatorType solver( ellipticOp, solverEps, solverEps );
438 // // Initializing CG-Solver
439 // // CAPG: it's very slow. Can't this be speed up somehow?
440 //
441 // dataOutput.write( timeProvider );
442 //
443 // ellipticOp.mass_matrix(U_n, rhs);
444 // // calculating: rhs = Mⁿ uⁿ (uⁿ == solutionContainer)
445 //
446 // // std::vector<double> timecontainer, l2errorcontainer, h1errorcontainer;
447 // int iter = 0;
448 // for(timeProvider.next(dT); iter < itno; timeProvider.next(dT), ++iter)
449 // // put your loop-action here, but not the last action
450 // {
451 // // for debugging reasons
452 // // std::cout << timeProvider.time() << '\t' << timeProvider.deltaT()
453 // // << "\tfor loop action." << std::endl;
454 // // l2h1error_helper_oss << timeProvider.time() << '\t' << timeProvider.deltaT()
455 // // << "\tfor loop action." << std::endl;
456 //
457 // f.setTime(timeProvider.time());
458 // // chance the time for the (long) RHS function f to tⁿ⁺¹
459 // // ('f' from -∆u + \div(v) u + \matdot{u} = f)
460 //
461 // initialData.setTime(timeProvider.time());
462 //
463 // std::array<double,2> time_array {
464 // {timeProvider.time()- timeProvider.deltaT(), timeProvider.time()}
465 // };
466 // EvoMapIntegrator impl_euler { NUMERIK::NewtonParameters<> {1e-10} };
467 // evoMap.evolve(time_array.begin(), time_array.end(), impl_euler);
468 //
469 // assembleRHS( f, tmp );
470 // // assemly stiffness/load vector; the vector is called 'tmp';
471 // // in the sense of above it is fⁿ⁺¹
472 //
473 // rhs.axpy( timeProvider.deltaT(), tmp );
474 // // rhs += Δt * tmp (i.e. rhs += \Delta t * tmp)
475 // // Just calculated: ( uⁿ + ∆t EllipticOperator(uⁿ)) + ∆t fⁿ
476 //
477 // //InterpolationType::interpolateFunction( initialData, U_n );
478 // solver( rhs, U_n );
479 // // Solve: (id + \Delta t EllipticOperator) u_{n+1} = rhs,
480 // // where rhs = ( u_n + \Delta t EllipticOperator(u_n)) + \Delta t f_n
481 // // CAPG remark: this is very very slow.
482 //
483 // dataOutput.write( timeProvider );
484 //
485 // InterpolationType::interpolateFunction( initialData, exactSolution );
486 // write_error(l2h1error_ofs);
487 //
488 // ellipticOp.mass_matrix(U_n, rhs);
489 // // calculating: rhs = Mⁿ uⁿ (uⁿ == solutionContainer)
490 //
491 // // l2errorcontainer.push_back(l2norm.distance(exactSolution, U_n));
492 // // h1errorcontainer.push_back(h1norm.distance(exactSolution, U_n));
493 // // timecontainer.push_back(timeProvider.time());
494 // }
495 // // l2h1error_ofs << "#\n" << l2h1error_helper_oss.str() << '#' << std::endl;
496 //
497 // // l2h1error_ofs << timeProvider.deltaT() << ' '
498 // // << timecontainer.back() << ' '
499 // // << l2errorcontainer.back() << ' '
500 // // << h1errorcontainer.back() << std::endl;
501 // l2h1error_ofs.close();
502 // }
503 // // end paperALE experiment
504 #endif // #ifndef DUNE_HEAT_ALGORITHM_HPP