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 // local includes
20 #include "dune_typedef_heat.hpp"
21 #include "/Users/christianpower/cpp/ODE_Solver/implicit_euler.h"
22 #include "/Users/christianpower/cpp/ODE_Solver/bdf.h"
23
24 // Remember: For every experiment modify
25 // in 'heat.hh'
26 // - RHSfunction
27 // - InitialData
28 // - HeatModel
29 // in 'elliptic.hh'
30 // - EllitpicOperator
31 // in 'deformation.hh'
32 // - DeformationCoordFunction
33
34 // DataOutputParameters
35 // --------------------
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
58 void BDF::ie_heat_algorithm(){
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 itno = (t_end - t_0)/dT + .1;
66
67 // prepare grid from DGF file
68 const std::string gridkey =
69 Dune::Fem::IOInterface::defaultGridKey( GridType::dimension );
70 const std::string gridfile =
71 Dune::Fem::Parameter::getValue< std::string >( gridkey );
72 if( Dune::Fem::MPIManager::rank() == 0 )
73 std::cout << "Loading macro grid: " << gridfile << std::endl;
74 Dune::GridPtr< HostGridType > hostGrid {gridfile};
75 hostGrid ->loadBalance();
76
77 // create grid
78 DeformationCoordFunction deformation {};
79 GridType grid(*hostGrid, deformation);
80 GridPartType gridPart(grid);
81 DiscreteFunctionSpaceType dfSpace(gridPart);
82
83 L2NormType l2norm(gridPart);
84 H1NormType h1norm(gridPart);
85 std::ofstream l2h1error_ofs
86 {Dune::Fem::Parameter::getValue<std::string>("fem.io.errorFile",
87 "../output/l2h1error")};
88
89 Dune::Fem::GridTimeProvider<GridType> timeProvider(t_0, grid);
90 timeProvider.init(dT); // Do your first action before you enter the for loop.
91
92 deformation.set_time_provider(timeProvider);
93 // deformation is aware of timeProvider
94
95 DiscreteFunctionType U_n("U_n",dfSpace), rhs("rhs", dfSpace),
96 load_vector("load_vector", dfSpace), exactSolution("exactSolution",dfSpace);
97
98 InitialDataType initialData {timeProvider};
99 InterpolationType::interpolateFunction( initialData, U_n );
100 InterpolationType::interpolateFunction( initialData, exactSolution );
101 // solutionContainer is used for higher order BDF methods,
102 // U_n == uⁿ, rhs == tmp1, tmp == tmp2
103 // (create discrete functions for intermediate functionals)
104
105 // to get the degree of freedom vector
106 // for(auto it = U_n.dbegin(); it != U_n.dend(); ++it)
107 // std::cout << *it << std::endl;
108
109 RHSFunctionType f {timeProvider};
110 // This expression is very long; it's literally the same f as in L(u) = f;
111 // together with the function/method in rhs.hh,
112
113 IOTupleType ioTuple(&U_n);
114 const int step = 0; // is there any resonable choise, whitout adaptivity?
115 DataOutputType
116 dataOutput(grid, ioTuple,
117 DataOutputParameters(Dune::Fem::Parameter::getValue<std::string>
118 ("fem.io.outputName",
119 "../output/ALE_LiteratureExample-"),
120 step) );
121 // dataOutput.write(timeProvider);
122
123 auto write_error = [&](std::ostream& os){
124 os << std::defaultfloat << timeProvider.time() << ' '
125 << std::scientific
126 << l2norm.distance(exactSolution, U_n) << ' '
127 << h1norm.distance(exactSolution, U_n) << std::endl;
128 };
129
130 // write_error(l2h1error_ofs);
131
132 NonlinearModel model {timeProvider};
133 NonlinearOperator ellipticOp {U_n, model};
134 const double solverEps =
135 Dune::Fem::Parameter::getValue<double>( "heat.solvereps", 1e-8 );
136 LinearInverseOperatorType solver(ellipticOp, solverEps, solverEps);
137 // Initializing CG-Solver
138 // CAPG: it's very slow. Can't this be speed up somehow?
139
140 // dataOutput.write(timeProvider);
141
142 ellipticOp.get_xi(U_n);
143 ellipticOp.mass_matrix(U_n, rhs);
144 // calculating: rhs = Mⁿ uⁿ
145
146 // debugging
147 std::cout << std::scientific << "\n==========\n"
148 << "U_n.dbegin() = " << *(U_n.dbegin()) << std::endl;
149 std::cout << "rhs.dbegin() = " << *(rhs.dbegin()) << std::endl;
150
151 int iter = 0;
152 for(timeProvider.next(dT); iter < itno; timeProvider.next(dT), ++iter)
153 // put your loop-action here, but not the last action
154 {
155 std::cout << std::defaultfloat << "==========\n"
156 << "In for loop: iter = "<< iter << ", "
157 << "time = " << timeProvider.time() << std::endl;
158
159 // f.setTime(timeProvider.time()); // obselete
160 // chance the time for the (long) RHS function f to tⁿ⁺¹
161 // ('f' from -∆u + \div(v) u + \matdot{u} = f)
162
163 // initialData.setTime(timeProvider.time()); // obselete
164
165 assembleRHS(f, load_vector);
166 std::cout << "assembleRHS(f, load_vector) = " << *(load_vector.dbegin()) << std::endl;
167
168 // assemly stiffness/load vector; the vector is called 'rhs';
169 // in the sense of above it is fⁿ⁺¹
170
171 rhs.axpy(timeProvider.deltaT(), load_vector);
172 // rhs += Δt * load_vector
173 // Just calculated: rhs == MⁿUⁿ + Δt fⁿ⁺¹
174
175 std::cout << "rhs.axpy = " << *(rhs.dbegin()) << std::endl;
176
177 //InterpolationType::interpolateFunction( initialData, U_n );
178 solver(rhs, U_n);
179 // Solve: (Mⁿ⁺¹ + Δt A(Uⁿ)) uⁿ⁺¹ = rhs,
180 // CAPG remark: this is very very slow.
181
182 std::cout << "solver(rhs, U_n) = " << *(U_n.dbegin()) << std::endl;
183
184 // dataOutput.write(timeProvider);
185
186 InterpolationType::interpolateFunction(initialData, exactSolution);
187 // write_error(l2h1error_ofs);
188
189 ellipticOp.mass_matrix(U_n, rhs);
190 std::cout << "mass_matrix(U_n, rhs) = " << *(rhs.dbegin()) << std::endl;
191
192 // calculating: rhs = Mⁿ uⁿ
193 ellipticOp.get_xi(U_n);
194
195 // debugging
196 std::cout << std::scientific
197 << "U_n.dbegin() = " << *(U_n.dbegin()) << std::endl;
198 std::cout << "rhs.dbegin() = " << *(rhs.dbegin()) << std::endl;
199 std::cout << "load_vector.dbegin() = " << *(load_vector.dbegin())
200 << std::endl;
201 }
202 l2h1error_ofs.close();
203 }
204
205 void BDF::bdf_algorithm(const int bdf_no){
206 if( bdf_no < 1 || bdf_no > 6)
207 throw std::runtime_error("ERROR in BDF::bdf_algorithm():"
208 " Bad bdf_no = " + std::to_string(bdf_no) + "\n"
209 "Enter an integer between 1 and 6");
210
211 double t_0 =
212 Dune::Fem::Parameter::getValue<double>("heat.starttime",0.0);
213 double dT =
214 Dune::Fem::Parameter::getValue<double>("heat.timestep",0.1);
215 double t_end =
216 Dune::Fem::Parameter::getValue<double>("heat.endtime",0.6);
217 const int time_step_no_max = (t_end - t_0)/dT + .1;
218
219 // prepare grid from DGF file
220 const std::string gridkey =
221 Dune::Fem::IOInterface::defaultGridKey(GridType::dimension);
222 const std::string gridfile =
223 Dune::Fem::Parameter::getValue<std::string>(gridkey);
224 if(Dune::Fem::MPIManager::rank() == 0)
225 std::cout << "Loading macro grid: " << gridfile << std::endl;
226 Dune::GridPtr<HostGridType> hostGrid {gridfile};
227 hostGrid->loadBalance();
228
229 // create grid
230 DeformationCoordFunction deformation {};
231 GridType grid(*hostGrid, deformation);
232 GridPartType gridPart(grid);
233 DiscreteFunctionSpaceType dfSpace(gridPart);
234
235 L2NormType l2norm(gridPart);
236 H1NormType h1norm(gridPart);
237 std::ofstream l2h1error_ofs
238 {Dune::Fem::Parameter::getValue<std::string>("fem.io.errorFile",
239 "../output/l2h1error")};
240
241 Dune::Fem::GridTimeProvider<GridType> timeProvider(t_0, grid);
242 timeProvider.init(dT); // Do your first action before you enter the for loop.
243
244 deformation.set_time_provider(timeProvider);
245 // deformation is aware of timeProvider
246
247 DiscreteFunctionType U_n("U_n",dfSpace), rhs("rhs", dfSpace),
248 load_vector("load_vector", dfSpace), exactSolution("exactSolution",dfSpace);
249 DiscreteFunctionType ie_U_n("ie_U_n", dfSpace);
250 DiscreteFunctionType ie_rhs("ie_rhs", dfSpace);
251
252 std::vector<DiscreteFunctionType> prev_M_n_U_n;
253 for(size_t i = 0; i < bdf_no; ++i)
254 prev_M_n_U_n.push_back(
255 DiscreteFunctionType {"U_n-" + std::to_string(i+1), dfSpace});
256
257
258 IOTupleType ioTuple(&U_n);
259 const int step = 0; // is there any resonable choise, whitout adaptivity?
260 DataOutputType
261 dataOutput(grid, ioTuple,
262 DataOutputParameters(Dune::Fem::Parameter::getValue<std::string>
263 ("fem.io.outputName",
264 "../output/ALE_LiteratureExample-"),
265 step) );
266 // dataOutput.write(timeProvider);
267
268 auto write_error = [&](std::ostream& os){
269 os << std::defaultfloat << timeProvider.time() << ' '
270 << std::scientific
271 << l2norm.distance(exactSolution, U_n) << ' '
272 << h1norm.distance(exactSolution, U_n) << std::endl;
273 };
274
275 InitialDataType initialData {timeProvider};
276 RHSFunctionType f {timeProvider};
277 NonlinearModel model {timeProvider};
278 NonlinearOperator ellipticOp {U_n, model};
279 const double solverEps =
280 Dune::Fem::Parameter::getValue<double>( "heat.solvereps", 1e-8 );
281 LinearInverseOperatorType solver(ellipticOp, solverEps, solverEps);
282
283 // first time step
284 int time_step_no = 0;
285 // timeProvider.time() == t_0;
286
287 InterpolationType::interpolateFunction(initialData, U_n);
288 InterpolationType::interpolateFunction(initialData, exactSolution);
289 InterpolationType::interpolateFunction(initialData, ie_U_n);
290
291 // write_error(l2h1error_ofs);
292 // dataOutput.write(timeProvider);
293
294 ellipticOp.get_xi(U_n);
295 ellipticOp.mass_matrix(ie_U_n, ie_rhs);
296 ellipticOp.mass_matrix(ie_U_n, prev_M_n_U_n.at(0) );
297 // calculating: rhs = Mⁿ uⁿ
298
299 // debugging
300 std::cout << std::scientific << "\n==========\n"
301 << "ie_U_n = " << *(ie_U_n.dbegin()) << '\t'
302 << "ie_rhs = " << *(ie_rhs.dbegin()) << '\n'
303 << "U_n = " << *(U_n.dbegin()) << '\t'
304 << "prev_M_n_U_n.at(0) = " << *(prev_M_n_U_n.at(0).dbegin()) << std::endl;
305
306
307 for(timeProvider.next(dT); time_step_no < time_step_no_max; timeProvider.next(dT), ++time_step_no)
308 // put your loop-action here, but not the last action
309 {
310 std::cout << std::defaultfloat << "==========\n"
311 << "In for loop: time_step_no = "<< time_step_no << ", "
312 << "time = " << timeProvider.time() << std::endl;
313
314 // f.setTime(timeProvider.time()); // obselete
315 // chance the time for the (long) RHS function f to tⁿ⁺¹
316 // ('f' from -∆u + \div(v) u + \matdot{u} = f)
317
318 // initialData.setTime(timeProvider.time()); // obselete
319
320 assembleRHS(f, load_vector);
321 std::cout << "assembleRHS(f, load_vector) = " << *(load_vector.dbegin()) << std::endl;
322
323 // assemly stiffness/load vector; the vector is called 'rhs';
324 // in the sense of above it is fⁿ⁺¹
325
326 rhs.axpy(timeProvider.deltaT(), load_vector);
327 // rhs += Δt * load_vector
328 // Just calculated: rhs == MⁿUⁿ + Δt fⁿ⁺¹
329
330 std::cout << "rhs.axpy = " << *(rhs.dbegin()) << std::endl;
331
332 //InterpolationType::interpolateFunction( initialData, U_n );
333 solver(rhs, U_n);
334 // Solve: (Mⁿ⁺¹ + Δt A(Uⁿ)) uⁿ⁺¹ = rhs,
335 // CAPG remark: this is very very slow.
336
337 std::cout << "solver(rhs, U_n) = " << *(U_n.dbegin()) << std::endl;
338
339 // dataOutput.write(timeProvider);
340
341 InterpolationType::interpolateFunction(initialData, exactSolution);
342 // write_error(l2h1error_ofs);
343
344 ellipticOp.mass_matrix(U_n, rhs);
345 std::cout << "mass_matrix(U_n, rhs) = " << *(rhs.dbegin()) << std::endl;
346
347 // calculating: rhs = Mⁿ uⁿ
348 ellipticOp.get_xi(U_n);
349
350 // debugging
351 std::cout << std::scientific
352 << "U_n.dbegin() = " << *(U_n.dbegin()) << std::endl;
353 std::cout << "rhs.dbegin() = " << *(rhs.dbegin()) << std::endl;
354 std::cout << "load_vector.dbegin() = " << *(load_vector.dbegin())
355 << std::endl;
356 }
357 l2h1error_ofs.close();
358 }
359
360
361 // begin 2014nonlinear experiment
362 void BDF::nonlinear_algorithm(){
363 double t_0 =
364 Dune::Fem::Parameter::getValue<double>("heat.starttime",0.0);
365 double dT =
366 Dune::Fem::Parameter::getValue<double>("heat.timestep",0.1);
367 double t_end =
368 Dune::Fem::Parameter::getValue<double>("heat.endtime",0.6);
369 const int itno = (t_end - t_0)/dT + .1;
370
371 // setting up BDF coefficients
372 const int bdf_no =
373 Dune::Fem::Parameter::getValue<double>("heat.bdf", 1);
374 const std::vector<double> bdf_alpha_coeff { NUMERIK::bdf_alpha_coeff(bdf_no) };
375 const std::vector<double> bdf_gamma_coeff { NUMERIK::bdf_gamma_coeff(bdf_no) };
376
377 // prepare grid from DGF file
378 const std::string gridkey =
379 Dune::Fem::IOInterface::defaultGridKey( GridType::dimension );
380 const std::string gridfile =
381 Dune::Fem::Parameter::getValue< std::string >( gridkey );
382 if( Dune::Fem::MPIManager::rank() == 0 )
383 std::cout << "Loading macro grid: " << gridfile << std::endl;
384 Dune::GridPtr< HostGridType > hostGrid {gridfile};
385 hostGrid ->loadBalance();
386
387 // create grid
388 DeformationCoordFunction deformation {};
389 GridType grid(*hostGrid, deformation);
390 GridPartType gridPart(grid);
391 DiscreteFunctionSpaceType dfSpace(gridPart);
392
393 L2NormType l2norm(gridPart);
394 H1NormType h1norm(gridPart);
395 std::ofstream l2h1error_ofs
396 {Dune::Fem::Parameter::getValue<std::string>("fem.io.errorFile",
397 "../output/l2h1error")};
398
399 Dune::Fem::GridTimeProvider<GridType> timeProvider(t_0, grid);
400 timeProvider.init(dT); // Do your first action before you enter the for loop.
401
402 deformation.set_time_provider(timeProvider);
403 // deformation is aware of timeProvider
404
405 // Setting up dof vectors
406 DiscreteFunctionType U_n("U_n", dfSpace); // Numerical solution
407 DiscreteFunctionType rhs("rhs", dfSpace); // Right hand side for the LES
408 DiscreteFunctionType bdf_tmp("bdf_tmp", dfSpace);
409 // Purpose of bdf_tmp:
410 // For the first k-time steps it is a tmp variable.
411 // Later it holds the sum from the previous BDF steps
412 DiscreteFunctionType load_vector("load_vector", dfSpace);
413 DiscreteFunctionType exact_solution("exact_solution", dfSpace);
414 DiscreteFunctionType err_vec("U_n_minus_exact_solution", dfSpace);
415 // Holds U_n - exact_solution
416
417 auto calc_err_vec = [&]
418 // helper function: caclulates err_vec - exact_solution
419 {
420 err_vec.assign(U_n);
421 err_vec.axpy((-1.), exact_solution);
422 };
423
424 std::vector<DiscreteFunctionType> prev_steps_vec;
425 for(int i = 0; i < bdf_no; ++i)
426 prev_steps_vec.push_back(
427 DiscreteFunctionType {"U_n" + std::to_string(i+1), dfSpace} );
428
429 // debugging
430 DiscreteFunctionType ie_U_n("ie_U_n", dfSpace);
431 DiscreteFunctionType ie_rhs("ie_rhs", dfSpace);
432 DiscreteFunctionType ie_load_vector("ie_load_vector", dfSpace);
433
434 // Setting up PDE/ discrete operator
435 InitialDataType initialData {timeProvider};
436 InterpolationType::interpolateFunction( initialData, ie_U_n );
437
438 RHSFunctionType f {timeProvider};
439 NonlinearModel model {timeProvider};
440 NonlinearOperator ellipticOp {ie_U_n, model};
441 const double solverEps =
442 Dune::Fem::Parameter::getValue<double>( "heat.solvereps", 1e-8 );
443 LinearInverseOperatorType solver(ellipticOp, solverEps, solverEps);
444
445 // Setting up paraview output
446 IOTupleType ioTuple(&U_n);
447 const int step = 0; // is there any resonable choise, whitout adaptivity?
448 DataOutputType
449 dataOutput(grid, ioTuple,
450 DataOutputParameters(Dune::Fem::Parameter::getValue<std::string>
451 ("fem.io.outputName",
452 "../output/ALE_LiteratureExample-"),
453 step) );
454 // helper function
455 auto write_error = [&](std::ostream& os){
456 os << std::defaultfloat << timeProvider.time() << ' '
457 << std::scientific
458 << l2norm.distance(exact_solution, U_n) << ' '
459 << h1norm.distance(exact_solution, U_n) << std::endl;
460 };
461
462
463
464 std::cout << '\n' << std::endl;
465 // initialize prev_step_vec[i] dof
466 int iter = (-1) * bdf_no;
467 for(int i =0;
468 (i < prev_steps_vec.size()) && (iter < itno);
469 ++i, ++iter, timeProvider.next(dT)){
470 std::cout << std::defaultfloat << "==========\n"
471 << "In initial steps: iter = "<< iter << ", "
472 << "time = " << timeProvider.time() << std::endl;
473
474 InterpolationType::interpolateFunction(initialData, bdf_tmp);
475 ellipticOp.mass_matrix(bdf_tmp, prev_steps_vec.at(i));
476 InterpolationType::interpolateFunction(initialData, exact_solution);
477 U_n.assign(bdf_tmp);
478
479 // debugging
480 // if(i == 0){
481 // InterpolationType::interpolateFunction(initialData, ie_U_n);
482 // std::cout << "interpolate = " << *(ie_U_n.dbegin()) << std::endl;
483 // ellipticOp.mass_matrix(ie_U_n, ie_rhs);
484 // ellipticOp.get_xi(ie_U_n);
485 // }
486 // else{
487 // assembleRHS(f, ie_load_vector);
488 // ie_rhs.axpy(timeProvider.deltaT(), ie_load_vector);
489 // solver(ie_rhs, ie_U_n);
490 // ellipticOp.get_xi(ie_U_n);
491 // }
492
493 // calc_err_vec();
494 // dataOutput.write(timeProvider);
495 // write_error(l2h1error_ofs);
496
497 std::cout << std::scientific
498 << "U_n.dbegin() = " << *(U_n.dbegin()) << '\t'
499 << "ie_U_n.dbegin() = " << *(ie_U_n.dbegin()) << '\n'
500 << "ie_rhs.dbegin() = " << *(ie_rhs.dbegin()) << std::endl;
501 }
502
503 ellipticOp.get_xi(U_n);
504 ellipticOp.mass_matrix(U_n, rhs);
505
506 // two helper functions
507
508 auto bdf_rhs = [&]
509 // calculating: rhs = (-1) · (∑ᵢ₌₁ᵏ αᵢ₋₁ (Mu)ⁿ⁻ᵏ⁺ⁱ)
510 // bdf_tmp = ∑ᵢ₌₁ᵏ γᵢ₋₁ (Mu)ⁿ⁻ᵏ⁺ⁱ
511 // NOTE for implicit euler: rhs = bdf_tmp = Mⁿ uⁿ
512 {
513
514 rhs.clear(); // set dof to zero
515 bdf_tmp.clear();
516
517 for(size_t i=0; i < prev_steps_vec.size(); ++i){
518 rhs.axpy(bdf_alpha_coeff.at(i), prev_steps_vec.at(i));
519 bdf_tmp.axpy(bdf_gamma_coeff.at(i), prev_steps_vec.at(i));
520 }
521 rhs *= (-1);
522 ellipticOp.get_xi(bdf_tmp);
523 };
524
525 auto bdf_cycle = [&]
526 // use this after calculating the numerical solution
527 {
528 // std::cout << "bdf_cycle at timeProvider.time() = "
529 // << timeProvider.time() << std::endl;
530 for(size_t i=0; i < prev_steps_vec.size() - 1; ++i)
531 prev_steps_vec.at(i).assign(prev_steps_vec.at(i+1));
532 ellipticOp.mass_matrix(U_n, prev_steps_vec.back());
533 };
534
535
536 for(;
537 iter < itno;
538 timeProvider.next(dT), ++iter)
539 // put your loop-action here, but not the last action
540 {
541 std::cout << std::defaultfloat << "==========\n"
542 << "In for loop: iter = "<< iter << ", "
543 << "time = " << timeProvider.time() << std::endl;
544
545 bdf_rhs();
546 assembleRHS(f, load_vector);
547 // assemly stiffness/load vector; the vector is called 'load_vector';
548 // in the sense of above it is fⁿ⁺¹
549
550 rhs.axpy(timeProvider.deltaT(), load_vector);
551 // rhs += Δt * load_vector
552 // Just calculated: rhs == MⁿUⁿ + Δt fⁿ⁺¹
553
554 rhs /= bdf_alpha_coeff.back();
555
556 //InterpolationType::interpolateFunction( initialData, U_n );
557 solver(rhs, U_n);
558 // Solve: (Mⁿ⁺¹ + Δt A(Uⁿ)) uⁿ⁺¹ = rhs,
559 // CAPG remark: this is very very slow.
560
561 bdf_cycle();
562 InterpolationType::interpolateFunction(initialData, exact_solution);
563 calc_err_vec();
564
565 // dataOutput.write(timeProvider);
566 // if(timeProvider.time() > t_end - 3*dT)
567 // write_error(l2h1error_ofs); // write_error(std::cout);
568
569 assembleRHS(f, ie_load_vector);
570 std::cout << "assembleRHS(f, ie_load_vector) = " << *(ie_load_vector.dbegin()) << std::endl;
571 ie_rhs.axpy(timeProvider.deltaT(), ie_load_vector);
572 std::cout << "ie_rhs.axpy = " << *(ie_rhs.dbegin()) << std::endl;
573 solver(ie_rhs, ie_U_n);
574 std::cout << "solver(ie_rhs, ie_U_n) = " << *(ie_U_n.dbegin()) << std::endl;
575 ellipticOp.mass_matrix(ie_U_n, ie_rhs);
576 std::cout << "mass_matrix(ie_U_n, ie_rhs) = " << *(ie_rhs.dbegin()) << std::endl;
577 ellipticOp.get_xi(ie_U_n);
578
579 // debugging
580 std::cout << std::scientific
581 << "U_n.dbegin() = " << *(U_n.dbegin()) << "\t\t"
582 << "ie_U_n.dbegin() = " << *(ie_U_n.dbegin()) << std::endl;
583 std::cout << "rhs.dbegin() = " << *(rhs.dbegin()) << "\t\t"
584 << "ie_rhs.dbegin() = " << *(ie_rhs.dbegin()) << std::endl;
585 std::cout << "load_vector.dbegin() = " << *(load_vector.dbegin()) << '\t'
586 << "ie_load_vector.dbegin() = " << *(ie_load_vector.dbegin())
587 << std::endl;
588 std::cout << "bdf_tmp.dbegin() = " << *(bdf_tmp.dbegin()) << std::endl;
589
590 for(size_t i=0; i < prev_steps_vec.size(); ++i){
591 std::cout << "prev_steps_vec.at(" << i << ").dbegin() = "
592 << *(prev_steps_vec.at(i).dbegin()) << std::endl;
593 }
594
595 }
596 l2h1error_ofs.close();
597 }
598 // end 2014nonlinear experiment
599
600
601 // // begin paperALE experiment
602 //
603 // // This has to be modified for the ODE solver
604 // struct EigenFullPivLuSolv{
605 // Eigen::Vector3d linear_solve(const Eigen::Matrix3d& A,
606 // const Eigen::Vector3d& b) const {
607 // return A.fullPivLu().solve(b);
608 // }
609 // };
610 // struct EigenNorm{
611 // double norm(const Eigen::Vector3d& v) const {
612 // return v.norm();
613 // }
614 // };
615 // struct SurfaceVectorfield{
616 // void set_time(const double time) { t = time; }
617 // double get_time() { return t; }
618 // Eigen::Vector3d evaluate(const Eigen::Vector3d vec) const {
619 // double x = vec(0), y = vec(1), z = vec(2);
620 // return Eigen::Vector3d {
621 // 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) ),
622 // 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)),
623 // 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) ) };
624 // }
625 // Eigen::Matrix3d jacobian(const Eigen::Vector3d vec) const {
626 // double x = vec(0), y = vec(1), z = vec(2);
627 // Eigen::Matrix3d jac_f_t_vec;
628 // jac_f_t_vec <<
629 // // first row
630 // 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)),
631 // -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)),
632 // -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)),
633 // // second row
634 // 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)),
635 // -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)),
636 // -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)),
637 // // third row
638 // 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)),
639 // -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)),
640 // -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));
641 // return jac_f_t_vec;
642 // }
643 // private:
644 // double t;
645 // };
646 // struct EigenIdMatrix3d{
647 // Eigen::Matrix3d identity_matrix() const { return Eigen::Matrix3d::Identity(); }
648 // };
649 //
650 //
651 // typedef NUMERIK::ImplicitEuler<EigenFullPivLuSolv, EigenNorm, SurfaceVectorfield,
652 // EigenIdMatrix3d, Eigen::Matrix3d, Eigen::Vector3d>
653 // EvoMapIntegrator;
654 //
655 // void BDF::algorithm()
656 // {
657 // /**********************************************************************
658 // * Pseudocode:
659 // * u⁰ = interpol(solution) >> plot
660 // * container = M⁰u⁰ >> save
661 // * for-loop n = 0 ↦ N with N·Δt + t_begin = t_end
662 // * t += ∆τ (set time)
663 // * tmp = container == Mⁿuⁿ
664 // * rhs = fⁿ⁺¹ == load Vector
665 // * tmp += Δτ · rhs (solution == Mⁿuⁿ + Δτ·rhs)
666 // * uⁿ⁺¹ == solution = (Mⁿ⁺¹ + Δτ Aⁿ⁺¹)⁻¹ tmp >> save
667 // * container = Mⁿ⁺¹ uⁿ⁺¹
668 // * end-for-loop
669 // *
670 // **********************************************************************/
671 //
672 // // create grid from DGF file
673 // const std::string gridkey =
674 // Dune::Fem::IOInterface::defaultGridKey( GridType::dimension );
675 // const std::string gridfile =
676 // Dune::Fem::Parameter::getValue< std::string >( gridkey );
677 // if( Dune::Fem::MPIManager::rank() == 0 )
678 // std::cout << "Loading macro grid: " << gridfile << std::endl;
679 //
680 // Dune::GridPtr< HostGridType > hostGrid {gridfile};
681 // hostGrid ->loadBalance();
682 //
683 // BDF::EvoMapType evoMap {gridfile};
684 // DeformationCoordFunction deformation {evoMap};
685 // // constructs an unsorted map with domain and range equal the vertexes from the
686 // // initial grid
687 //
688 // GridType grid(*hostGrid, deformation);
689 // GridPartType gridPart(grid);
690 // DiscreteFunctionSpaceType dfSpace(gridPart);
691 //
692 // L2NormType l2norm(gridPart);
693 // H1NormType h1norm(gridPart);
694 // std::ofstream l2h1error_ofs
695 // {Dune::Fem::Parameter::getValue<std::string>("fem.io.errorFile",
696 // "../output/l2h1error")};
697 // //std::ios_base::app};
698 // //l2h1error_ofs << std::scientific;
699 // //l2h1error_ofs << "L2Error\tH1Error" << std::scientific << std::endl;
700 // //std::ostringstream l2h1error_helper_oss;
701 // //l2h1error_helper_oss << "time\ttimestep" << std::endl;
702 //
703 // double t_0 =
704 // Dune::Fem::Parameter::getValue<double>("heat.starttime",0.0);
705 // double dT =
706 // Dune::Fem::Parameter::getValue<double>("heat.timestep",0.1);
707 // double t_end =
708 // Dune::Fem::Parameter::getValue<double>("heat.endtime",0.6);
709 // std::cout << (t_end - t_0)/dT + .1 << std::endl;
710 // const int itno = (t_end - t_0)/dT + .1;
711 //
712 // Dune::Fem::GridTimeProvider< GridType > timeProvider(t_0, grid);
713 // timeProvider.init(dT); // Do your first action before you enter the for loop.
714 // // std::cout << timeProvider.time() << '\t'
715 // // << timeProvider.deltaT() << "\tstarting action."
716 // // << std::endl;
717 // // l2h1error_helper_oss << timeProvider.time() << '\t'
718 // // << timeProvider.deltaT() << "\tstarting action."
719 // // << std::endl;
720 //
721 // DiscreteFunctionType solutionContainer( "U_nContainer", dfSpace ),
722 // U_n("U_n",dfSpace), rhs( "rhs", dfSpace ), tmp( "tmp", dfSpace ),
723 // exactSolution("exactSolution",dfSpace);
724 // InitialDataType initialData;
725 // InterpolationType::interpolateFunction( initialData, U_n );
726 // InterpolationType::interpolateFunction( initialData, exactSolution );
727 // // solutionContainer is used for higher order BDF methods,
728 // // U_n == uⁿ, rhs == tmp1, tmp == tmp2
729 // // (create discrete functions for intermediate functionals)
730 //
731 // RHSFunctionType f;
732 // // This expression is very long; it's literally the same f as in L(u) = f;
733 // // together with the function/method in rhs.hh,
734 //
735 // IOTupleType ioTuple(&U_n);
736 // const int step = 0; // is there any resonable choise, whitout adaptivity?
737 // DataOutputType
738 // dataOutput(grid, ioTuple,
739 // DataOutputParameters(Dune::Fem::Parameter::getValue<std::string>
740 // ("fem.io.outputName",
741 // "../output/ALE_LiteratureExample-"),
742 // step) );
743 //
744 // auto write_error = [&](std::ostream& os){
745 // os << std::defaultfloat << timeProvider.time() << ' '
746 // << std::scientific
747 // << l2norm.distance(exactSolution, U_n) << ' '
748 // << h1norm.distance(exactSolution, U_n) << std::endl;
749 // };
750 // write_error(l2h1error_ofs);
751 //
752 // EllipticOperatorType
753 // ellipticOp( HeatModelType {timeProvider} );
754 //
755 // const double solverEps =
756 // Dune::Fem::Parameter::getValue< double >( "heat.solvereps", 1e-8 );
757 // LinearInverseOperatorType solver( ellipticOp, solverEps, solverEps );
758 // // Initializing CG-Solver
759 // // CAPG: it's very slow. Can't this be speed up somehow?
760 //
761 // dataOutput.write( timeProvider );
762 //
763 // ellipticOp.mass_matrix(U_n, rhs);
764 // // calculating: rhs = Mⁿ uⁿ (uⁿ == solutionContainer)
765 //
766 // // std::vector<double> timecontainer, l2errorcontainer, h1errorcontainer;
767 // int iter = 0;
768 // for(timeProvider.next(dT); iter < itno; timeProvider.next(dT), ++iter)
769 // // put your loop-action here, but not the last action
770 // {
771 // // for debugging reasons
772 // // std::cout << timeProvider.time() << '\t' << timeProvider.deltaT()
773 // // << "\tfor loop action." << std::endl;
774 // // l2h1error_helper_oss << timeProvider.time() << '\t' << timeProvider.deltaT()
775 // // << "\tfor loop action." << std::endl;
776 //
777 // f.setTime(timeProvider.time());
778 // // chance the time for the (long) RHS function f to tⁿ⁺¹
779 // // ('f' from -∆u + \div(v) u + \matdot{u} = f)
780 //
781 // initialData.setTime(timeProvider.time());
782 //
783 // std::array<double,2> time_array {
784 // {timeProvider.time()- timeProvider.deltaT(), timeProvider.time()}
785 // };
786 // EvoMapIntegrator impl_euler { NUMERIK::NewtonParameters<> {1e-10} };
787 // evoMap.evolve(time_array.begin(), time_array.end(), impl_euler);
788 //
789 // assembleRHS( f, tmp );
790 // // assemly stiffness/load vector; the vector is called 'tmp';
791 // // in the sense of above it is fⁿ⁺¹
792 //
793 // rhs.axpy( timeProvider.deltaT(), tmp );
794 // // rhs += Δt * tmp (i.e. rhs += \Delta t * tmp)
795 // // Just calculated: ( uⁿ + ∆t EllipticOperator(uⁿ)) + ∆t fⁿ
796 //
797 // //InterpolationType::interpolateFunction( initialData, U_n );
798 // solver( rhs, U_n );
799 // // Solve: (id + \Delta t EllipticOperator) u_{n+1} = rhs,
800 // // where rhs = ( u_n + \Delta t EllipticOperator(u_n)) + \Delta t f_n
801 // // CAPG remark: this is very very slow.
802 //
803 // dataOutput.write( timeProvider );
804 //
805 // InterpolationType::interpolateFunction( initialData, exactSolution );
806 // write_error(l2h1error_ofs);
807 //
808 // ellipticOp.mass_matrix(U_n, rhs);
809 // // calculating: rhs = Mⁿ uⁿ (uⁿ == solutionContainer)
810 //
811 // // l2errorcontainer.push_back(l2norm.distance(exactSolution, U_n));
812 // // h1errorcontainer.push_back(h1norm.distance(exactSolution, U_n));
813 // // timecontainer.push_back(timeProvider.time());
814 // }
815 // // l2h1error_ofs << "#\n" << l2h1error_helper_oss.str() << '#' << std::endl;
816 //
817 // // l2h1error_ofs << timeProvider.deltaT() << ' '
818 // // << timecontainer.back() << ' '
819 // // << l2errorcontainer.back() << ' '
820 // // << h1errorcontainer.back() << std::endl;
821 // l2h1error_ofs.close();
822 // }
823 // // end paperALE experiment
824
825 #endif // #ifndef DUNE_HEAT_ALGORITHM_HPP