Discontinuous Galerkin Library
#include "dg/algorithm.h"
Loading...
Searching...
No Matches
runge_kutta.h
Go to the documentation of this file.
1#ifndef _DG_RK_
2#define _DG_RK_
3
4#include <cassert>
5#include <array>
6#include <tuple>
7
8#include "ode.h"
10#include "tableau.h"
11#include "blas1.h"
12#include "implicit.h"
13
18namespace dg{
20namespace detail{
21// this is a dense matrix-matrix multiplication
22// a generalization of the corresponding blas2::gemv DenseMatrix algorithm
23template<class ContainerType>
24void gemm(
25 const std::array<double,2>& alpha,
26 const DenseMatrix<ContainerType>& m,
27 const std::array<const std::vector<double>*,2> x,
28 const std::array<double,2>& beta,
29 std::array<ContainerType*,2>&& y)
30{
31 const unsigned size = m.num_cols();
32 unsigned i=0;
33 if( size >= 8)
34 for( i=0; i<size/8; i++)
36 *y[0], *y[1], i==0 ? beta[0] : 1., i==0 ? beta[1] : 1.,
37 alpha[0]*(*x[0])[i*8+0], alpha[1]*(*x[1])[i*8+0], *m.get()[i*8+0],
38 alpha[0]*(*x[0])[i*8+1], alpha[1]*(*x[1])[i*8+1], *m.get()[i*8+1],
39 alpha[0]*(*x[0])[i*8+2], alpha[1]*(*x[1])[i*8+2], *m.get()[i*8+2],
40 alpha[0]*(*x[0])[i*8+3], alpha[1]*(*x[1])[i*8+3], *m.get()[i*8+3],
41 alpha[0]*(*x[0])[i*8+4], alpha[1]*(*x[1])[i*8+4], *m.get()[i*8+4],
42 alpha[0]*(*x[0])[i*8+5], alpha[1]*(*x[1])[i*8+5], *m.get()[i*8+5],
43 alpha[0]*(*x[0])[i*8+6], alpha[1]*(*x[1])[i*8+6], *m.get()[i*8+6],
44 alpha[0]*(*x[0])[i*8+7], alpha[1]*(*x[1])[i*8+7], *m.get()[i*8+7]);
45 unsigned l=0;
46 if( size%8 >= 4)
47 for( l=0; l<(size%8)/4; l++)
49 *y[0], *y[1], size < 8 ? beta[0] : 1., size < 8 ? beta[1] : 1.,
50 alpha[0]*(*x[0])[i*8+l*4+0], alpha[1]*(*x[1])[i*8+l*4+0], *m.get()[i*8+l*4+0],
51 alpha[0]*(*x[0])[i*8+l*4+1], alpha[1]*(*x[1])[i*8+l*4+1], *m.get()[i*8+l*4+1],
52 alpha[0]*(*x[0])[i*8+l*4+2], alpha[1]*(*x[1])[i*8+l*4+2], *m.get()[i*8+l*4+2],
53 alpha[0]*(*x[0])[i*8+l*4+3], alpha[1]*(*x[1])[i*8+l*4+3], *m.get()[i*8+l*4+3]);
54 unsigned k=0;
55 if( (size%8)%4 >= 2)
56 for( k=0; k<((size%8)%4)/2; k++)
58 *y[0], *y[1], size < 4 ? beta[0] : 1., size < 4 ? beta[1] : 1.,
59 alpha[0]*(*x[0])[i*8+l*4+k*2+0], alpha[1]*(*x[1])[i*8+l*4+k*2+0], *m.get()[i*8+l*4+k*2+0],
60 alpha[0]*(*x[0])[i*8+l*4+k*2+1], alpha[1]*(*x[1])[i*8+l*4+k*2+1], *m.get()[i*8+l*4+k*2+1]);
61 if( ((size%8)%4)%2 == 1)
63 *y[0], *y[1], size < 2 ? beta[0] : 1., size < 2 ? beta[1] : 1.,
64 alpha[0]*(*x[0])[i*8+l*4+k*2], alpha[1]*(*x[1])[i*8+l*4+k*2], *m.get()[i*8+l*4+k*2]);
65
66}
67} // namespace detail
69
119{
127 template<class ContainerType1>
128 void operator()( ContainerType1&) const{ }
129};
131template<class ContainerType>
132struct FilteredERKStep;
134//Should we try if filters inside RHS are equally usable?
135
138
139
162template< class ContainerType>
164{
166 using container_type = ContainerType;
168 ERKStep() = default;
177 ERKStep( ConvertsToButcherTableau<value_type> tableau, const ContainerType&
178 copyable): m_ferk(tableau, copyable)
179 {
180 }
181
183 template<class ...Params>
184 void construct( Params&& ...ps)
185 {
186 //construct and swap
187 *this = ERKStep( std::forward<Params>( ps)...);
188 }
190 const ContainerType& copyable()const{ return m_ferk.copyable();}
191
193 void ignore_fsal(){ m_ferk.ignore_fsal();}
195 void enable_fsal(){ m_ferk.enable_fsal();}
196
200 template<class ExplicitRHS>
201 void step( ExplicitRHS& rhs, value_type t0, const ContainerType& u0, value_type& t1, ContainerType& u1, value_type dt, ContainerType& delta)
202 {
204 m_ferk.step( std::tie( rhs, id), t0, u0, t1, u1, dt, delta);
205 }
232 template<class ExplicitRHS>
233 void step( ExplicitRHS& rhs, value_type t0, const ContainerType& u0, value_type&
234 t1, ContainerType& u1, value_type dt)
235 {
237 m_ferk.step( std::tie( rhs, id), t0, u0, t1, u1, dt);
238 }
240 unsigned order() const {
241 return m_ferk.order();
242 }
244 unsigned embedded_order() const {
245 return m_ferk.embedded_order();
246 }
248 unsigned num_stages() const{
249 return m_ferk.num_stages();
250 }
251 private:
253};
254
276template< class ContainerType>
278{
280 using container_type = ContainerType;
282 FilteredERKStep() { m_k.resize(1); }
285 ContainerType& copyable): m_rk(tableau), m_k(m_rk.num_stages(),
286 copyable)
287 {
288 m_rkb.resize(m_k.size()), m_rkd.resize( m_k.size());
289 for( unsigned i=0; i<m_k.size(); i++)
290 {
291 m_rkb[i] = m_rk.b(i);
292 m_rkd[i] = m_rk.d(i);
293 }
294 }
295
297 template<class ...Params>
298 void construct( Params&& ...ps)
299 {
300 //construct and swap
301 *this = FilteredERKStep( std::forward<Params>( ps)...);
302 }
304 const ContainerType& copyable()const{ return m_k[0];}
305
307 void ignore_fsal(){ m_ignore_fsal = true;}
309 void enable_fsal(){ m_ignore_fsal = false;}
310
315 template<class ExplicitRHS, class Limiter>
316 void step( const std::tuple<ExplicitRHS,Limiter>& rhs, value_type t0, const ContainerType& u0, value_type& t1, ContainerType& u1, value_type dt, ContainerType& delta)
317 {
318 step ( rhs, t0, u0, t1, u1, dt, delta, true);
319 }
322 template<class ExplicitRHS, class Limiter>
323 void step( const std::tuple<ExplicitRHS, Limiter>& rhs, value_type t0, const ContainerType& u0, value_type&
324 t1, ContainerType& u1, value_type dt)
325 {
326 if( !m_tmp_allocated)
327 {
328 dg::assign( m_k[0], m_tmp);
329 m_tmp_allocated = true;
330 }
331 step ( rhs, t0, u0, t1, u1, dt, m_tmp, false);
332 }
334 unsigned order() const {
335 return m_rk.order();
336 }
338 unsigned embedded_order() const {
339 return m_rk.embedded_order();
340 }
342 unsigned num_stages() const{
343 return m_rk.num_stages();
344 }
345 private:
346 template<class ExplicitRHS, class Limiter>
347 void step( const std::tuple<ExplicitRHS, Limiter>& rhs, value_type t0, const ContainerType& u0, value_type& t1, ContainerType& u1, value_type dt, ContainerType& delta, bool);
349 std::vector<value_type> m_rkb, m_rkd;
350 std::vector<ContainerType> m_k;
351 value_type m_t1 = 1e300;//remember the last timestep at which ERK is called
352 bool m_ignore_fsal = false;
353 ContainerType m_tmp; // only conditionally allocated
354 bool m_tmp_allocated = false;
355};
356
358
359template< class ContainerType>
360template<class ExplicitRHS, class Limiter>
361void FilteredERKStep<ContainerType>::step( const std::tuple<ExplicitRHS, Limiter>& ode, value_type t0, const ContainerType& u0, value_type& t1, ContainerType& u1, value_type dt, ContainerType& delta, bool compute_delta)
362{
363 unsigned s = m_rk.num_stages();
364 std::vector<const ContainerType*> k_ptrs = dg::asPointers( m_k);
365 //0 stage: probe
366 value_type tu = t0;
367 if( t0 != m_t1 || m_ignore_fsal)
368 std::get<0>(ode)(t0, u0, m_k[0]); //freshly compute k_0
369 //else take from last call
370 for ( unsigned i=1; i<s; i++)
371 {
372 std::vector<value_type> rka( i);
373 for( unsigned l=0; l<i; l++)
374 rka[l] = m_rk.a(i,l);
375
376 tu = DG_FMA( dt,m_rk.c(i),t0); //l=0
377 dg::blas1::copy( u0, delta); // can't use u1 here cause u0 can alias
378 dg::blas2::gemv( dt, dg::asDenseMatrix( k_ptrs, i), rka, 1., delta);
379
380 std::get<1>(ode)( delta);
381 std::get<0>(ode)( tu, delta, m_k[i]);
382 }
383 //Now add everything up to get solution and error estimate
384 dg::blas1::copy( u0, u1);
385 if( compute_delta)
386 detail::gemm( {dt,dt}, dg::asDenseMatrix(k_ptrs), {&m_rkb, &m_rkd},
387 {1.,0.}, {&u1, &delta});
388 else
389 blas2::gemv( dt, dg::asDenseMatrix(k_ptrs), m_rkb, 1., u1);
390 std::get<1>(ode)( u1);
391 //make sure (t1,u1) is the last call to f
392 m_t1 = t1 = t0 + dt;
393 if(!m_rk.isFsal() )
394 std::get<0>(ode)( t1, u1, m_k[0]);
395 else
396 {
397 using std::swap;
398 swap( m_k[0], m_k[s-1]); //enable free swap functions
399 }
400}
402
403
428template<class ContainerType>
430{
432 using container_type = ContainerType;
434 ARKStep(){ m_kI.resize(1); }
441 ARKStep( std::string name, const ContainerType& copyable)
442 {
443 std::string exp_name = name+" (explicit)";
444 std::string imp_name = name+" (implicit)";
445 *this = ARKStep( exp_name, imp_name, copyable);
446 }
463 const ContainerType& copyable
464 ):
465 m_rkE(ex_tableau),
466 m_rkI(im_tableau),
467 m_kE(m_rkE.num_stages(), copyable),
468 m_kI(m_rkI.num_stages(), copyable)
469 {
470 assert( m_rkE.num_stages() == m_rkI.num_stages());
471 // check fsal
472 assert( m_rkI.a(0,0) == 0);
473 assert( m_rkI.c(m_rkI.num_stages()-1) == 1);
474 check_implicit_fsal();
475 assign_coeffs();
476 }
478 template<class ...Params>
479 void construct( Params&& ...ps)
480 {
481 //construct and swap
482 *this = ARKStep( std::forward<Params>( ps)...);
483 }
485 const ContainerType& copyable()const{ return m_kI[0];}
486
509 template< class ExplicitRHS, class ImplicitRHS, class Solver>
510 void step( const std::tuple<ExplicitRHS, ImplicitRHS, Solver>& ode, value_type t0, const ContainerType& u0, value_type& t1, ContainerType& u1, value_type dt, ContainerType& delta);
512 unsigned order() const {
513 return m_rkE.order();
514 }
516 unsigned embedded_order() const {
517 return m_rkE.order();
518 }
520 unsigned num_stages() const{
521 return m_rkE.num_stages();
522 }
523 private:
524 ButcherTableau<value_type> m_rkE, m_rkI;
525 std::vector<ContainerType> m_kE, m_kI;
526 std::vector<value_type> m_rkb, m_rkd;
527 value_type m_t1 = 1e300;
528 bool m_implicit_fsal = false;
529 void check_implicit_fsal(){
530 m_implicit_fsal = true;
531 for( unsigned i=0; i<m_rkI.num_stages(); i++)
532 if( m_rkI.a(i,0) != 0)
533 m_implicit_fsal = false;
534 }
535 void assign_coeffs()
536 {
537 m_rkb.resize( 2*m_rkI.num_stages());
538 m_rkd.resize( 2*m_rkI.num_stages());
539 for( unsigned i=0; i<m_rkI.num_stages(); i++)
540 {
541 m_rkb[2*i] = m_rkE.b(i);
542 m_rkb[2*i+1] = m_rkI.b(i);
543 m_rkd[2*i] = m_rkE.d(i);
544 m_rkd[2*i+1] = m_rkI.d(i);
545 }
546 }
547};
548
550template<class ContainerType>
551template< class Explicit, class Implicit, class Solver>
552void ARKStep<ContainerType>::step( const std::tuple<Explicit,Implicit,Solver>& ode, value_type t0, const ContainerType& u0, value_type& t1, ContainerType& u1, value_type dt, ContainerType& delta)
553{
554 unsigned s = m_rkE.num_stages();
555 value_type tu = t0;
556 //0 stage
558 if( t0 != m_t1)
559 std::get<0>(ode)(t0, u0, m_kE[0]); //freshly compute k_0
560 if( !m_implicit_fsal) // all a(i,0) == 0
561 std::get<1>(ode)(t0, u0, m_kI[0]);
562 // DO NOT HOLD POINTERS AS PRVIATE MEMBERS
563 std::vector<const ContainerType*> k_ptrs( 2*m_rkI.num_stages());
564 for( unsigned i=0; i<m_rkI.num_stages(); i++)
565 {
566 k_ptrs[2*i ] = &m_kE[i];
567 k_ptrs[2*i+1] = &m_kI[i];
568 }
569
570 for( unsigned i=1; i<s; i++)
571 {
572 std::vector<value_type> rka( 2*i);
573 for(unsigned l=0; l<i; l++)
574 {
575 rka[2*l] = m_rkE.a(i,l);
576 rka[2*l+1] = m_rkI.a(i,l);
577 }
578 tu = DG_FMA( m_rkI.c(i),dt, t0);
579 dg::blas1::copy( u0, m_kI[i]);
580 dg::blas2::gemv( dt, dg::asDenseMatrix( k_ptrs, 2*i), rka, 1., m_kI[i]);
581 value_type alpha = dt*m_rkI.a(i,i);
582 std::get<2>(ode)( alpha, tu, delta, m_kI[i]);
583 dg::blas1::axpby( 1./alpha, delta, -1./alpha, m_kI[i]);
584 std::get<0>(ode)(tu, delta, m_kE[i]);
585 }
586 m_t1 = t1 = tu;
587 //Now compute result and error estimate
588
589 dg::blas1::copy( u0, u1);
590 detail::gemm( {dt,dt}, dg::asDenseMatrix(k_ptrs), {&m_rkb, &m_rkd},
591 {1.,0.}, {&u1, &delta});
592 //make sure (t1,u1) is the last call to ex
593 std::get<0>(ode)(t1,u1,m_kE[0]);
594}
596
614template<class ContainerType>
616{
617 //MW Dirk methods cannot have stage order greater than 1
619 using container_type = ContainerType;
621 DIRKStep(){ m_kI.resize(1); }
622
634 const ContainerType& copyable):
635 m_rkI(im_tableau),
636 m_kI(m_rkI.num_stages(), copyable)
637 {
638 m_rkIb.resize(m_kI.size()), m_rkId.resize(m_kI.size());
639 for( unsigned i=0; i<m_kI.size(); i++)
640 {
641 m_rkIb[i] = m_rkI.b(i);
642 m_rkId[i] = m_rkI.d(i);
643 }
644 }
645
647 template<class ...Params>
648 void construct( Params&& ...ps)
649 {
650 //construct and swap
651 *this = DIRKStep( std::forward<Params>( ps)...);
652 }
654 const ContainerType& copyable()const{ return m_kI[0];}
655
662 template< class ImplicitRHS, class Solver>
663 void step( const std::tuple<ImplicitRHS,Solver>& ode, value_type t0, const
664 ContainerType& u0, value_type& t1, ContainerType& u1, value_type dt,
665 ContainerType& delta)
666 {
667 step( ode, t0, u0, t1, u1, dt, delta, true);
668 }
685 template< class ImplicitRHS, class Solver>
686 void step( const std::tuple<ImplicitRHS, Solver>& ode, value_type t0, const
687 ContainerType& u0, value_type& t1, ContainerType& u1, value_type dt)
688 {
689 if( !m_allocated)
690 {
691 dg::assign( m_kI[0], m_tmp);
692 m_allocated = true;
693 }
694 step( ode, t0, u0, t1, u1, dt, m_tmp, false);
695 }
697 unsigned order() const {
698 return m_rkI.order();
699 }
701 unsigned embedded_order() const {
702 return m_rkI.order();
703 }
705 unsigned num_stages() const{
706 return m_rkI.num_stages();
707 }
708
709 private:
710 template< class ImplicitRHS, class Solver>
711 void step( const std::tuple<ImplicitRHS, Solver>& ode, value_type t0, const
712 ContainerType& u0, value_type& t1, ContainerType& u1, value_type
713 dt, ContainerType& delta, bool compute_delta);
715 std::vector<ContainerType> m_kI;
716 ContainerType m_tmp;
717 std::vector<value_type> m_rkIb, m_rkId;
718 bool m_allocated = false;
719};
720
722template<class ContainerType>
723template< class ImplicitRHS, class Solver>
724void DIRKStep<ContainerType>::step( const std::tuple<ImplicitRHS,Solver>& ode, value_type t0, const ContainerType& u0, value_type& t1, ContainerType& u1, value_type dt, ContainerType& delta, bool compute_delta)
725{
726 unsigned s = m_rkI.num_stages();
727 value_type tu = t0;
728 //0 stage
729 //rhs = u0
730 tu = DG_FMA( m_rkI.c(0),dt, t0);
732 if( m_rkI.a(0,0) !=0 )
733 {
734 alpha = dt*m_rkI.a(0,0);
735 std::get<1>(ode)( alpha, tu, delta, u0);
736 dg::blas1::axpby( 1./alpha, delta, -1./alpha, u0, m_kI[0]);
737 }
738 else
739 std::get<0>(ode)(tu, u0, m_kI[0]);
740 std::vector<const ContainerType*> kIptr = dg::asPointers( m_kI);
741
742 for( unsigned i=1; i<s; i++)
743 {
744 tu = DG_FMA( m_rkI.c(i),dt, t0);
745 dg::blas1::copy( u0, m_kI[i]);
746 std::vector<value_type> rkIa( i);
747 for( unsigned l=0; l<i; l++)
748 rkIa[l] = m_rkI.a(i,l);
749 dg::blas2::gemv( dt, dg::asDenseMatrix(kIptr,i), rkIa, 1., m_kI[i]);
750 if( m_rkI.a(i,i) !=0 )
751 {
752 alpha = dt*m_rkI.a(i,i);
753 std::get<1>(ode)( alpha, tu, delta, m_kI[i]);
754 dg::blas1::axpby( 1./alpha, delta, -1./alpha, m_kI[i]);
755 }
756 else
757 std::get<0>(ode)(tu, delta, m_kI[i]);
758 }
759 t1 = t0 + dt;
760 //Now compute result and error estimate
761 dg::blas1::copy( u0, u1);
762 if( compute_delta)
763 detail::gemm( {dt,dt}, dg::asDenseMatrix(kIptr), {&m_rkIb, &m_rkId},
764 {1.,0.}, {&u1, &delta});
765 else
766 blas2::gemv( dt, dg::asDenseMatrix(kIptr), m_rkIb, 1., u1);
767}
769
794template<class ContainerType>
796
807template<class ContainerType>
809
842template<class ContainerType>
844{
846 using container_type = ContainerType;
848 ShuOsher(){m_u.resize(1);}
857 ShuOsher( dg::ConvertsToShuOsherTableau<value_type> tableau, const ContainerType& copyable): m_t( tableau), m_u( m_t.num_stages(), copyable), m_k(m_u)
858 { }
860 template<class ...Params>
861 void construct( Params&& ...ps)
862 {
863 //construct and swap
864 *this = ShuOsher( std::forward<Params>( ps)...);
865 }
867 const ContainerType& copyable()const{ return m_u[0];}
868
883 template<class ExplicitRHS, class Limiter>
884 void step( const std::tuple<ExplicitRHS, Limiter>& ode, value_type t0, const ContainerType& u0, value_type& t1, ContainerType& u1, value_type dt){
885 unsigned s = m_t.num_stages();
886 std::vector<value_type> ts( m_t.num_stages()+1);
887 ts[0] = t0;
888 dg::blas1::copy(u0, m_u[0]);
889 if( t0 != m_t1 ) //this is the first time we call step
890 std::get<0>(ode)(ts[0], m_u[0], m_k[0]); //freshly compute k_0
891 for( unsigned i=1; i<=s; i++)
892 {
893 dg::blas1::axpbypgz( m_t.alpha(i-1,0), m_u[0], dt*m_t.beta(i-1,0), m_k[0], 0., i==s ? u1 : m_u[i]);
894 ts[i] = m_t.alpha(i-1,0)*ts[0] + dt*m_t.beta(i-1,0);
895 for( unsigned j=1; j<i; j++)
896 {
897 //about the i-1: it is unclear to me how the ShuOsher tableau makes implicit schemes
898 dg::blas1::axpbypgz( m_t.alpha(i-1,j), m_u[j], dt*m_t.beta(i-1,j), m_k[j], 1., i==s ? u1 : m_u[i]);
899 ts[i] += m_t.alpha(i-1,j)*ts[j] + dt*m_t.beta(i-1,j);
900
901 }
902 std::get<1>(ode)( i==s ? u1 : m_u[i]);
903 if(i!=s)
904 std::get<0>(ode)(ts[i], m_u[i], m_k[i]);
905 else
906 //make sure (t1,u1) is the last call to f
907 std::get<0>(ode)(ts[i], u1, m_k[0]);
908 }
909 m_t1 = t1 = ts[s];
910 }
912 unsigned order() const {
913 return m_t.order();
914 }
916 unsigned num_stages() const{
917 return m_t.num_stages();
918 }
919 private:
921 std::vector<ContainerType> m_u, m_k;
922 value_type m_t1 = 1e300;
923};
943template<class ContainerType>
945
946
949inline bool is_same( double x, double y, double eps = 1e-15)
950{
951 return fabs(x - y) < eps * std::max(1.0, std::max( fabs(x), fabs(y)));
952}
955inline bool is_same( float x, float y, float eps = 1e-6)
956{
957 return fabsf(x - y) < eps * std::max(1.0f, std::max( fabsf(x), fabsf(y)));
958}
959
962template<class T>
963inline bool is_same( T x, T y)
964{
965 return x == y;
966}
967
971inline bool is_divisable( double a, double b, double eps = 1e-15)
972{
973 return is_same( round(a/b)*b, a, eps);
974}
978inline bool is_divisable( float a, float b, float eps = 1e-6)
979{
980 return is_same( (float)round(a/b)*b, (float)a, eps);
981}
982
998template<class ContainerType>
999struct SinglestepTimeloop : public aTimeloop<ContainerType>
1000{
1001 using container_type = ContainerType;
1005
1012 SinglestepTimeloop( std::function<void ( value_type, const ContainerType&,
1013 value_type&, ContainerType&, value_type)> step, value_type dt = 0 )
1014 : m_step(step), m_dt(dt){}
1015
1030 template<class Stepper, class ODE>
1031 SinglestepTimeloop( Stepper&& stepper, ODE&& ode, value_type dt = 0)
1032 {
1033 // Open/Close Principle (OCP) Software entities should be open for extension but closed for modification
1034 m_step = [=, cap = std::tuple<Stepper, ODE>(std::forward<Stepper>(stepper),
1035 std::forward<ODE>(ode)) ]( auto t0, const auto& y0, auto& t1, auto& y1,
1036 auto dtt) mutable
1037 {
1038 std::get<0>(cap).step( std::get<1>(cap), t0, y0, t1, y1, dtt);
1039 };
1040 m_dt = dt;
1041 }
1042
1044 template<class ...Params>
1045 void construct( Params&& ...ps)
1046 {
1047 //construct and swap
1048 *this = SinglestepTimeloop( std::forward<Params>( ps)...);
1049 }
1050
1056 void set_dt( value_type dt){ m_dt = dt;}
1057
1073 container_type& u1, unsigned steps)
1074 {
1075 set_dt( (t1-t0)/(value_type)steps );
1076 this->integrate( t0, u0, t1, u1);
1077 }
1078
1079 virtual SinglestepTimeloop* clone() const{ return new
1080 SinglestepTimeloop(*this);}
1081 private:
1082 virtual void do_integrate(value_type& t0, const container_type& u0, value_type t1, container_type& u1, enum to mode) const;
1083 std::function<void ( value_type, const ContainerType&, value_type&,
1084 ContainerType&, value_type)> m_step;
1085 virtual value_type do_dt( ) const { return m_dt;}
1086 value_type m_dt;
1087};
1088
1090template< class ContainerType>
1091void SinglestepTimeloop<ContainerType>::do_integrate(
1092 value_type& t_begin, const ContainerType&
1093 begin, value_type t_end, ContainerType& end,
1094 enum to mode ) const
1095{
1096 bool forward = (t_end - t_begin > 0);
1097 if( (m_dt < 0 && forward) || ( m_dt > 0 && !forward) )
1098 throw dg::Error( dg::Message(_ping_)<<"Timestep has wrong sign! dt "<<m_dt);
1099 dg::blas1::copy( begin, end);
1100 if( m_dt == 0)
1101 throw dg::Error( dg::Message(_ping_)<<"Timestep may not be zero in SinglestepTimeloop!");
1102 if( is_divisable( t_end-t_begin, m_dt))
1103 {
1104 unsigned N = (unsigned)round((t_end - t_begin)/m_dt);
1105 for( unsigned i=0; i<N; i++)
1106 m_step( t_begin, end, t_begin, end, m_dt);
1107 }
1108 else
1109 {
1110 unsigned N = (unsigned)floor( (t_end-t_begin)/m_dt);
1111 for( unsigned i=0; i<N; i++)
1112 m_step( t_begin, end, t_begin, end, m_dt);
1113 if( dg::to::exact == mode)
1114 {
1115 value_type dt_final = t_end - t_begin;
1116 m_step( t_begin, end, t_begin, end, dt_final);
1117 }
1118 else
1119 m_step( t_begin, end, t_begin, end, m_dt);
1120 }
1121 return;
1122}
1124
1126
1127} //namespace dg
1128
1129#endif //_DG_RK_
class intended for the use in throw statements
Definition exceptions.h:83
small class holding a stringstream
Definition exceptions.h:29
Error classes or the dg library.
#define _ping_
Definition exceptions.h:12
void rhs(double t, double, double &yp)
bool is_divisable(double a, double b, double eps=1e-15)
Definition runge_kutta.h:971
bool is_same(double x, double y, double eps=1e-15)
Definition runge_kutta.h:949
void copy(const ContainerTypeIn &source, ContainerTypeOut &target)
Definition blas1.h:243
void axpbypgz(value_type alpha, const ContainerType1 &x, value_type1 beta, const ContainerType2 &y, value_type2 gamma, ContainerType &z)
Definition blas1.h:337
void axpby(value_type alpha, const ContainerType1 &x, value_type1 beta, ContainerType &y)
Definition blas1.h:306
void subroutine(Subroutine f, ContainerType &&x, ContainerTypes &&... xs)
; Customizable and generic blas1 function
Definition blas1.h:677
void assign(const from_ContainerType &from, ContainerType &to, Params &&... ps)
Generic way to assign the contents of a from_ContainerType object to a ContainerType object optionall...
Definition blas1.h:767
void gemv(get_value_type< ContainerType1 > alpha, MatrixType &&M, const ContainerType1 &x, get_value_type< ContainerType1 > beta, ContainerType2 &y)
Alias for blas2::symv ;.
Definition blas2.h:339
@ forward
forward derivative (cell to the right and current cell)
Definition enums.h:98
@ y
y direction
@ x
x direction
std::vector< const ContainerType * > asPointers(const std::vector< ContainerType > &in)
Convert a vector of vectors to a vector of pointers.
Definition densematrix.h:110
auto asDenseMatrix(const std::vector< const ContainerType * > &in)
Lightweight DenseMatrix for dg::blas2::gemv.
Definition densematrix.h:75
typename TensorTraits< std::decay_t< Vector > >::value_type get_value_type
Definition tensor_traits.h:45
to
Switch for the Timeloop integrate function.
Definition ode.h:17
@ exact
match the ending exactly
const double m
const double alpha
This is the namespace for all functions and classes defined and used by the discontinuous Galerkin li...
An abstract ODE integrator.
Additive Runge Kutta (semi-implicit) time-step with error estimate following The ARKode library
Definition runge_kutta.h:430
ContainerType container_type
Definition runge_kutta.h:432
void step(const std::tuple< ExplicitRHS, ImplicitRHS, Solver > &ode, value_type t0, const ContainerType &u0, value_type &t1, ContainerType &u1, value_type dt, ContainerType &delta)
Advance one step.
get_value_type< ContainerType > value_type
the value type of the time variable (float or double)
Definition runge_kutta.h:431
ARKStep()
No memory allocation.
Definition runge_kutta.h:434
unsigned order() const
global order of the method given by the current Butcher Tableau
Definition runge_kutta.h:512
ARKStep(ConvertsToButcherTableau< value_type > ex_tableau, ConvertsToButcherTableau< value_type > im_tableau, const ContainerType &copyable)
Construct with two Butcher Tableaus.
Definition runge_kutta.h:461
void construct(Params &&...ps)
Perfect forward parameters to one of the constructors.
Definition runge_kutta.h:479
const ContainerType & copyable() const
Return an object of same size as the object used for construction.
Definition runge_kutta.h:485
ARKStep(std::string name, const ContainerType &copyable)
Construct with given name.
Definition runge_kutta.h:441
unsigned embedded_order() const
global order of the embedding given by the current Butcher Tableau
Definition runge_kutta.h:516
unsigned num_stages() const
number of stages of the method given by the current Butcher Tableau
Definition runge_kutta.h:520
Manage coefficients of a (extended) Butcher tableau.
Definition tableau.h:33
unsigned order() const
global order of accuracy for the method represented by b
Definition tableau.h:132
real_type c(unsigned i) const
Read the c_i coefficients.
Definition tableau.h:99
unsigned embedded_order() const
global order of accuracy for the embedded method represented by bt
Definition tableau.h:136
real_type a(unsigned i, unsigned j) const
Read the a_ij coefficients.
Definition tableau.h:91
real_type d(unsigned j) const
Return the coefficients for the error estimate Equivalent to b(j)-bt(j)
Definition tableau.h:124
real_type b(unsigned j) const
Read the b_j coefficients.
Definition tableau.h:107
unsigned num_stages() const
The number of stages s.
Definition tableau.h:128
Convert identifiers to their corresponding dg::ButcherTableau.
Definition tableau.h:1695
Convert identifiers to their corresponding dg::ShuOsherTableau.
Definition tableau.h:1743
Embedded diagonally implicit Runge Kutta time-step with error estimate .
Definition runge_kutta.h:616
DIRKStep(ConvertsToButcherTableau< value_type > im_tableau, const ContainerType &copyable)
Construct with a diagonally implicit Butcher Tableau.
Definition runge_kutta.h:633
unsigned order() const
global order of the method given by the current Butcher Tableau
Definition runge_kutta.h:697
void step(const std::tuple< ImplicitRHS, Solver > &ode, value_type t0, const ContainerType &u0, value_type &t1, ContainerType &u1, value_type dt, ContainerType &delta)
Advance one step with error estimate.
Definition runge_kutta.h:663
get_value_type< ContainerType > value_type
the value type of the time variable (float or double)
Definition runge_kutta.h:618
unsigned embedded_order() const
global order of the embedding given by the current Butcher Tableau
Definition runge_kutta.h:701
ContainerType container_type
Definition runge_kutta.h:619
DIRKStep()
No memory allocation.
Definition runge_kutta.h:621
void step(const std::tuple< ImplicitRHS, Solver > &ode, value_type t0, const ContainerType &u0, value_type &t1, ContainerType &u1, value_type dt)
Advance one step ignoring error estimate and embedded method.
Definition runge_kutta.h:686
void construct(Params &&...ps)
Perfect forward parameters to one of the constructors.
Definition runge_kutta.h:648
unsigned num_stages() const
number of stages of the method given by the current Butcher Tableau
Definition runge_kutta.h:705
const ContainerType & copyable() const
Return an object of same size as the object used for construction.
Definition runge_kutta.h:654
Embedded Runge Kutta explicit time-step with error estimate .
Definition runge_kutta.h:164
unsigned num_stages() const
number of stages of the method given by the current Butcher Tableau
Definition runge_kutta.h:248
void step(ExplicitRHS &rhs, value_type t0, const ContainerType &u0, value_type &t1, ContainerType &u1, value_type dt, ContainerType &delta)
Advance one step with error estimate.
Definition runge_kutta.h:201
void enable_fsal()
All subsequent calls to step method will enable the check for the first same as last property.
Definition runge_kutta.h:195
const ContainerType & copyable() const
Return an object of same size as the object used for construction.
Definition runge_kutta.h:190
ContainerType container_type
Definition runge_kutta.h:166
void construct(Params &&...ps)
Perfect forward parameters to one of the constructors.
Definition runge_kutta.h:184
ERKStep()=default
No memory allocation.
void ignore_fsal()
All subsequent calls to step method will ignore the first same as last property (useful if you want t...
Definition runge_kutta.h:193
ERKStep(ConvertsToButcherTableau< value_type > tableau, const ContainerType &copyable)
Reserve internal workspace for the integration.
Definition runge_kutta.h:177
get_value_type< ContainerType > value_type
the value type of the time variable (float or double)
Definition runge_kutta.h:165
void step(ExplicitRHS &rhs, value_type t0, const ContainerType &u0, value_type &t1, ContainerType &u1, value_type dt)
Advance one step ignoring error estimate and embedded method.
Definition runge_kutta.h:233
unsigned embedded_order() const
global order of the embedding given by the current Butcher Tableau
Definition runge_kutta.h:244
unsigned order() const
global order of the method given by the current Butcher Tableau
Definition runge_kutta.h:240
Definition subroutines.h:174
EXPERIMENTAL: Filtered Embedded Runge Kutta explicit time-step with error estimate .
Definition runge_kutta.h:278
get_value_type< ContainerType > value_type
the value type of the time variable (float or double)
Definition runge_kutta.h:279
unsigned order() const
global order of the method given by the current Butcher Tableau
Definition runge_kutta.h:334
FilteredERKStep()
No memory allocation.
Definition runge_kutta.h:282
void step(const std::tuple< ExplicitRHS, Limiter > &rhs, value_type t0, const ContainerType &u0, value_type &t1, ContainerType &u1, value_type dt)
Advance one step ignoring error estimate and embedded method.
Definition runge_kutta.h:323
void step(const std::tuple< ExplicitRHS, Limiter > &rhs, value_type t0, const ContainerType &u0, value_type &t1, ContainerType &u1, value_type dt, ContainerType &delta)
Advance one step with error estimate.
Definition runge_kutta.h:316
void ignore_fsal()
All subsequent calls to step method will ignore the first same as last property (useful if you want t...
Definition runge_kutta.h:307
void enable_fsal()
All subsequent calls to step method will enable the check for the first same as last property.
Definition runge_kutta.h:309
void construct(Params &&...ps)
Perfect forward parameters to one of the constructors.
Definition runge_kutta.h:298
unsigned embedded_order() const
global order of the embedding given by the current Butcher Tableau
Definition runge_kutta.h:338
ContainerType container_type
Definition runge_kutta.h:280
unsigned num_stages() const
number of stages of the method given by the current Butcher Tableau
Definition runge_kutta.h:342
const ContainerType & copyable() const
Return an object of same size as the object used for construction.
Definition runge_kutta.h:304
FilteredERKStep(ConvertsToButcherTableau< value_type > tableau, const ContainerType &copyable)
Reserve internal workspace for the integration.
Definition runge_kutta.h:284
A filter that does nothing.
Definition runge_kutta.h:119
void operator()(ContainerType1 &) const
Do nothing.
Definition runge_kutta.h:128
Shu-Osher fixed-step explicit ODE integrator with Slope Limiter / Filter .
Definition runge_kutta.h:844
ContainerType container_type
Definition runge_kutta.h:846
void step(const std::tuple< ExplicitRHS, Limiter > &ode, value_type t0, const ContainerType &u0, value_type &t1, ContainerType &u1, value_type dt)
Advance one step.
Definition runge_kutta.h:884
const ContainerType & copyable() const
Return an object of same size as the object used for construction.
Definition runge_kutta.h:867
ShuOsher()
No memory allocation.
Definition runge_kutta.h:848
ShuOsher(dg::ConvertsToShuOsherTableau< value_type > tableau, const ContainerType &copyable)
Reserve internal workspace for the integration.
Definition runge_kutta.h:857
unsigned order() const
global order of the method given by the current Butcher Tableau
Definition runge_kutta.h:912
get_value_type< ContainerType > value_type
the value type of the time variable (float or double)
Definition runge_kutta.h:845
void construct(Params &&...ps)
Perfect forward parameters to one of the constructors.
Definition runge_kutta.h:861
unsigned num_stages() const
number of stages of the method given by the current Butcher Tableau
Definition runge_kutta.h:916
Manage coefficients in Shu-Osher form.
Definition tableau.h:181
real_type beta(unsigned i, unsigned j)
Read the beta_ij coefficients.
Definition tableau.h:257
unsigned num_stages() const
The number of stages s.
Definition tableau.h:259
unsigned order() const
global order of accuracy for the method
Definition tableau.h:263
real_type alpha(unsigned i, unsigned j)
Read the alpha_ij coefficients.
Definition tableau.h:250
Integrate using a for loop and a fixed time-step.
Definition runge_kutta.h:1000
SinglestepTimeloop(Stepper &&stepper, ODE &&ode, value_type dt=0)
Bind the step function of a single step stepper.
Definition runge_kutta.h:1031
void construct(Params &&...ps)
Perfect forward parameters to one of the constructors.
Definition runge_kutta.h:1045
SinglestepTimeloop(std::function< void(value_type, const ContainerType &, value_type &, ContainerType &, value_type)> step, value_type dt=0)
Construct using a std::function.
Definition runge_kutta.h:1012
virtual SinglestepTimeloop * clone() const
Abstract copy method that returns a copy of *this on the heap.
Definition runge_kutta.h:1079
void integrate_steps(value_type t0, const container_type &u0, value_type t1, container_type &u1, unsigned steps)
Integrate differential equation with a fixed number of steps.
Definition runge_kutta.h:1072
ContainerType container_type
Definition runge_kutta.h:1001
void set_dt(value_type dt)
Set the constant timestep to be used in the integrate functions.
Definition runge_kutta.h:1056
dg::get_value_type< ContainerType > value_type
Definition runge_kutta.h:1002
SinglestepTimeloop()=default
no allocation
Abstract timeloop independent of stepper and ODE.
Definition ode.h:60
void integrate(value_type t0, const ContainerType &u0, value_type t1, ContainerType &u1)
Integrate a differential equation between given bounds.
Definition ode.h:83
double value_type