gear_module contains adaptive high-order Gear solver More...

Functions/Subroutines
subroutine, public	init_gear_solver (Y, t_init, h_init)
	Initialize iteration variables. More...

subroutine, public	init_dydt (dYdt)
	Initialize dYdt component of the Nordsieck vector. More...

real(r_kind) function, public	get_time ()
	The current time. More...

real(r_kind) function, public	get_timestep ()
	The current timestep. More...

subroutine, public	set_timestep (timestep)
	Set the current timestep. More...

subroutine, public	revert_timestep (timestep)
	Revert a timestep. More...

real(r_kind) function, public	get_l1 ()
	The current l_1 (in Fortran: l(2)) More...

real(r_kind) function, dimension(net_size), public	get_solution ()
	Gives a copy of the solution at the current time. More...

real(r_kind) function, dimension(net_size), public	get_predictor_y ()
	Determines the predicted y_next. More...

real(r_kind) function, dimension(net_size), public	get_predictor_dydt ()
	Determines the predicted dydt_next. More...

subroutine, public	get_solution_at (time_inter, Y_inter)
	Determines the solution at a given point in time for [ti-h, ti]. More...

subroutine, public	shift_tau
	Shifts tau vector and controls the current time. More...

subroutine, public	nordsieck_product
	Nordsieck vector-matrix product. More...

subroutine, public	set_xi
	Function to calculate xi (needed for getting l in the corrector step) More...

subroutine, public	set_l
	Function to calculate \( \ell \) (needed for the corrector step) More...

subroutine, public	set_new_result (ydiff)
	Updates the solution (corrector step) More...

subroutine, public	prepare_next_step
	Stepsize and order control to pepare the next step. More...

subroutine, private	geterrors
	Calculate errors (helper functions below) More...

subroutine	errorq
	Calculate the errors for the current order. More...

subroutine, private	errorqm1
	Calculate the errors for the order (q-1) More...

subroutine, private	errorqp1
	Calculate the errors for the order (q+1) More...

subroutine, private	qfun (helpprod, qval)
	Calculate Helper function Q. More...

subroutine	cfun (qp1, xi_in, helpprod, cval)
	Calculate helper function C. More...

subroutine, private	shiftorder
	Adjusts the order, if possible and necessary. This is done every qth step, or earlier in case the time step was modified. More...

real(r_kind) function, private	get_cfactor (q_val)
	Function to set the conservative timestep factor depending on the current order q. More...

Variables
integer, parameter, private	histsize = 13
	length of history in tau, e, ... More...

integer, parameter, private	qmax = 5
	maximum order of q More...

real(r_kind), dimension(6, 6), private	amat = reshape( (/ 1.0d0, 1.0d0, 1.0d0, 1.0d0, 1.0d0, 1.0d0, 0.0d0, 1.0d0, 2.0d0, 3.0d0, 4.0d0, 5.0d0, 0.0d0, 0.0d0, 1.0d0, 3.0d0, 6.0d0,10.0d0, 0.0d0, 0.0d0, 0.0d0, 1.0d0, 4.0d0,10.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 1.0d0, 5.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 1.0d0 /), (/6, 6/) )
	Helper Matrix. More...

integer, private	ifac

integer, dimension(qmax+1), private	helperarr = (/ (ifac, ifac=1,qmax+1) /)

integer, dimension(qmax+1), private	factorial

integer, public	gear_nr_count = 0
	counter for NR steps More...

integer, private	q

integer, private	n

integer, private	nq

real(r_kind), private	ti

real(r_kind), private	h

integer, private	q_ini = 1
	start with order 1 More...

integer, private	n_ini = 1
	iteration step More...

integer, private	nq_ini = 1
	iteration step at current order More...

real(r_kind), dimension(histsize), private	tau
	hlist inverted More...

real(r_kind), dimension(:,:), allocatable, private	z
	Nordsieck vector It is defined as. More...

real(r_kind), dimension(:,:), allocatable, private	znew
	Predictor step. More...

real(r_kind), dimension(histsize), private	xi

real(r_kind), dimension(:), allocatable, private	en

real(r_kind), dimension(:,:), allocatable, private	e

real(r_kind), dimension(histsize), private	l

real(r_kind), dimension(histsize), private	el

real(r_kind), dimension(:), allocatable, private	errq

real(r_kind), dimension(:), allocatable, private	errqm1

real(r_kind), dimension(:), allocatable, private	errqp1

integer, private	alloc_stat

real(r_kind), private	hq

real(r_kind), private	hqp1

real(r_kind), private	hqm1

Detailed Description

gear_module contains adaptive high-order Gear solver

The solver is described by, e.g., Byrne & Hindmarsh (1975), ACM TOMS 1(1), p.71. The basic idea is to estimate the most efficient way between a larger timestep, but a higher order (and therefore computationally more expensive) or a smaller timestep and a lower order solution. The Gear solver therefore calculates the timestep based on an error estimation originating from the solution using different orders.

Author: Dirk Martin

Date: 01.08.17

See also: Byrne & Hindmarsh 1975, Longland et al. 2014, Martin, D. 2017

Edited:

OK 16.08.17: privatised local fields and methods
MR 18.04.19: fixed gear for termination criterion 0
MR 15.01.21: fixed gear for all termination criterions in timestep_module
MR 22.01.21: Extented comments
MR 17.07.22: Fixed a bug related to switching the orders. This bug was only visible with the compiler flag check_bounds which lead to slightly different results
MR 17.02.23: Introduced "revert_timestep" to revert the timestep if the solution did not converge

Function/Subroutine Documentation

◆ cfun()

subroutine gear_module::cfun	(	integer, intent(in)	qp1,
		real(r_kind), dimension(histsize), intent(in)	xi_in,
		real(r_kind), intent(in)	helpprod,
		real(r_kind), intent(out)	cval
	)

Calculate helper function C.

The function is defined as:

\[ C_{n+1} = \frac{\prod \limits _{i=1}^{q} \xi _{i}}{(q+1)!}\left[1+ \prod \limits _{i=2}^{q} \frac{t_{n+1} - t_{n+1-i}}{t_n -t_{n+1-i}} \right]. \]

Author: D. Martin

Date: 01.08.17

Definition at line 747 of file gear_module.f90.

◆ errorq()

subroutine gear_module::errorq

Calculate the errors for the current order.

The error is defined as:

\[ \vec{E}_{n+1}(q) = -\frac{1}{\ell _1}\left[ 1+ \prod _{i=2}^{q}\left( \frac{t_{n+1}-t_{n+1-i}}{t_{n}-t_{n+1-i}} \right) \right] ^{-1} \vec{e}_{n+1} \]

Author: D. Martin

Date: 01.08.17

Definition at line 597 of file gear_module.f90.

◆ errorqm1()

subroutine, private gear_module::errorqm1

private

Calculate the errors for the order (q-1)

The error is defined as:

\[ \vec{E}_{n+1}(q-1) = - \frac{\prod \limits _{i=1}^{q-1} \xi _i } {\ell _1 (q-1)} \frac{h^q \vec{y}_{n+1}^{(q)}}{q!} \]

Author: D. Martin

Date: 01.08.17

Definition at line 632 of file gear_module.f90.

◆ errorqp1()

subroutine, private gear_module::errorqp1

private

Calculate the errors for the order (q+1)

The error is defined as:

\[ \vec{E}_{n+1}(q+1) = \frac{-\xi_{q+1}(\vec{e}_{n+1}-Q_{n+1}\vec{e}_n) } {(q+2)\ell _1 (q+1)\left[1+ \prod \limits_{i=2}^{q} \frac{t_{n+1}-t_{n+1-i}} {t_n-t_{n+1-i}} \right] } \]

Author: D. Martin

Date: 01.08.17

Definition at line 665 of file gear_module.f90.

Here is the call graph for this function:

◆ get_cfactor()

real(r_kind) function, private gear_module::get_cfactor ( integer, intent(in) q_val )

private

Function to set the conservative timestep factor depending on the current order q.

Note: This function is not used!

Author: D. Martin

Date: 01.08.17

Definition at line 835 of file gear_module.f90.

◆ get_l1()

real(r_kind) function, public gear_module::get_l1

The current l_1 (in Fortran: l(2))

See also: set_l

Author: D. Martin

Date: 01.08.17

Definition at line 280 of file gear_module.f90.

◆ get_predictor_dydt()

real(r_kind) function, dimension(net_size), public gear_module::get_predictor_dydt

Determines the predicted dydt_next.

This returns the predicted dydt from the predictor step from z_new.

Author: D. Martin

Date: 01.08.17

Definition at line 331 of file gear_module.f90.

◆ get_predictor_y()

real(r_kind) function, dimension(net_size), public gear_module::get_predictor_y

Determines the predicted y_next.

The predicted abundance is stored in z_new(1,:)

Author: D. Martin

Date: 01.08.17

Definition at line 313 of file gear_module.f90.

◆ get_solution()

real(r_kind) function, dimension(net_size), public gear_module::get_solution

Gives a copy of the solution at the current time.

The solution abundance is stored at z(1,:).

Author: D. Martin

Date: 01.08.17

Definition at line 296 of file gear_module.f90.

◆ get_solution_at()

subroutine, public gear_module::get_solution_at	(	real(r_kind), intent(in)	time_inter,
		real(r_kind), dimension(net_size), intent(out)	Y_inter
	)

Determines the solution at a given point in time for [ti-h, ti].

Author: D. Martin

Date: 01.08.17

Definition at line 345 of file gear_module.f90.

◆ get_time()

real(r_kind) function, public gear_module::get_time

The current time.

Author: D. Martin

Date: 01.08.17

Definition at line 195 of file gear_module.f90.

◆ get_timestep()

real(r_kind) function, public gear_module::get_timestep

The current timestep.

Author: D. Martin

Date: 01.08.17

Definition at line 209 of file gear_module.f90.

◆ geterrors()

subroutine, private gear_module::geterrors

private

Calculate errors (helper functions below)

Calculates the errors (E, E(q+1), and E(q-1)) by calling the respective subroutines errorq, errorqm1, errorqp1. Also estimates the resulting timestep from each order.

Author: D. Martin

Date: 01.08.17

Edited:

M. Reichert, 15.07.22: Set time step of qp1 to zero if already in maximum order.

Definition at line 545 of file gear_module.f90.

Here is the call graph for this function:

◆ init_dydt()

subroutine, public gear_module::init_dydt ( real(r_kind), dimension(net_size), intent(in) dYdt )

Initialize dYdt component of the Nordsieck vector.

Stores h * dYdt in the Nordsieck vector z(2,:).

Author: D. Martin

Date: 01.08.17

Definition at line 180 of file gear_module.f90.

◆ init_gear_solver()

subroutine, public gear_module::init_gear_solver	(	real(r_kind), dimension(net_size), intent(in)	Y,
		real(r_kind), intent(in)	t_init,
		real(r_kind), intent(in)	h_init
	)

Initialize iteration variables.

Here, all necessary allocations are done. And the initial abundance is written into the abundance history array (z(1,:)).

Author: D. Martin

Date: 01.08.17

Definition at line 130 of file gear_module.f90.

Here is the call graph for this function:

◆ nordsieck_product()

subroutine, public gear_module::nordsieck_product

Nordsieck vector-matrix product.

This calculates:

\[ \vec{z}_{n+1}^{(0)} = \mathbf{A} \cdot \vec{z}_n \]

and stores the result in z_new.

Author: D. Martin

Date: 01.08.17

Definition at line 392 of file gear_module.f90.

◆ prepare_next_step()

subroutine, public gear_module::prepare_next_step

Stepsize and order control to pepare the next step.

Author: D. Martin

Date: 01.08.17

Definition at line 502 of file gear_module.f90.

Here is the call graph for this function:

◆ qfun()

subroutine, private gear_module::qfun	(	real(r_kind), intent(in)	helpprod,
		real(r_kind), intent(out)	qval
	)

private

Calculate Helper function Q.

Q is defined as:

\[ Q_{n+1} =\frac{C_{n+1}}{C_n}\left( \frac{h_{n+1}}{h_n}\right)^{q+1} \\ \]

Author: D. Martin

Date: 01.08.17

Definition at line 706 of file gear_module.f90.

Here is the call graph for this function:

◆ revert_timestep()

subroutine, public gear_module::revert_timestep ( real(r_kind), intent(in) timestep )

Revert a timestep.

This routine reverts the timestep of Gear to the previous timestep. It is used in case the result did not converge

See also: timestep_module::advance_gear

Author: M. Reichert

Definition at line 261 of file gear_module.f90.

Here is the call graph for this function:

◆ set_l()

subroutine, public gear_module::set_l

Function to calculate \( \ell \) (needed for the corrector step)

\[ \begin{align} \ell _0 (q) &=1, \nonumber \\ \quad \ell_1(q)&=\sum \limits _{i=1}^{q} \left( \xi _i ^{-1} \right), \nonumber\\ \quad \ell_j (q) &= \ell_j(q-1)+\ell_{j-1}(q-1)/\xi _q, \nonumber\\ \quad \ell _q (q) &= \left( \prod \limits_{i=1}^{q}\xi _i \right) ^{-1}. \nonumber \end{align} \]

Author: D. Martin

Date: 01.08.17

Definition at line 447 of file gear_module.f90.

◆ set_new_result()

subroutine, public gear_module::set_new_result ( real(r_kind), dimension(net_size), intent(in) ydiff )

Updates the solution (corrector step)

This routine calculates

\[ \vec{z}_{n+1}=\vec{z}_{n+1}^{(0)}+\vec{e}_{n+1}\vec{\ell} \]

See also: set_l, timestep_module::advance_gear

Author: D. Martin

Date: 01.08.17

Definition at line 475 of file gear_module.f90.

◆ set_timestep()

subroutine, public gear_module::set_timestep ( real(r_kind), intent(in) timestep )

Set the current timestep.

Warning: When changing the timestep of Gear, the Nordsieck vector (z) has to be rescaled because it depends on the timestep. The timestep of gear should therefore only changed by this routine!

See also: timestep_module::restrict_timestep, network_init_module::switch_evolution_mode.

Author: D. Martin

Date: 01.08.17

Definition at line 232 of file gear_module.f90.

◆ set_xi()

subroutine, public gear_module::set_xi

Function to calculate xi (needed for getting l in the corrector step)

xi is defined as:

\[ \xi _i = \frac{t_{n+1} - t_{n+1-i}}{h} \]

Author: D. Martin

Date: 01.08.17

Definition at line 420 of file gear_module.f90.

◆ shift_tau()

subroutine, public gear_module::shift_tau

Shifts tau vector and controls the current time.

Author: D. Martin

Date: 01.08.17

Definition at line 370 of file gear_module.f90.

◆ shiftorder()

subroutine, private gear_module::shiftorder

private

Adjusts the order, if possible and necessary. This is done every qth step, or earlier in case the time step was modified.

Author: D. Martin

Date: 01.08.17

Definition at line 782 of file gear_module.f90.

Variable Documentation

◆ alloc_stat

integer, private gear_module::alloc_stat

private

Definition at line 95 of file gear_module.f90.

◆ amat

real(r_kind), dimension(6,6), private gear_module::amat = reshape( (/ 1.0d0, 1.0d0, 1.0d0, 1.0d0, 1.0d0, 1.0d0, 0.0d0, 1.0d0, 2.0d0, 3.0d0, 4.0d0, 5.0d0, 0.0d0, 0.0d0, 1.0d0, 3.0d0, 6.0d0,10.0d0, 0.0d0, 0.0d0, 0.0d0, 1.0d0, 4.0d0,10.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 1.0d0, 5.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 1.0d0 /), (/6, 6/) )

private

Helper Matrix.

This Matrix is defined by:

\[ A_{ij}=\begin{cases}0 & \text{if } i<j\\\begin{pmatrix}i\\j\end{pmatrix} = \frac{i!}{j!(i-j)!} & \text{if } i\ge j\end{cases} \qquad \text{with } i,j \in [0,1,...,q] \]

The Taylor-series of \(\vec{y}_n\) truncated at order q, can be calculated as the product of the Nordsieck vector (z).

Definition at line 48 of file gear_module.f90.

◆ e

real(r_kind), dimension(:,:), allocatable, private gear_module::e

private

Definition at line 91 of file gear_module.f90.

◆ el

real(r_kind), dimension(histsize), private gear_module::el

private

Definition at line 92 of file gear_module.f90.

◆ en

real(r_kind), dimension(:), allocatable, private gear_module::en

private

Definition at line 90 of file gear_module.f90.

◆ errq

real(r_kind), dimension(:), allocatable, private gear_module::errq

private

Definition at line 93 of file gear_module.f90.

◆ errqm1

real(r_kind), dimension(:), allocatable, private gear_module::errqm1

private

Definition at line 93 of file gear_module.f90.

◆ errqp1

real(r_kind), dimension(:), allocatable, private gear_module::errqp1

private

Definition at line 93 of file gear_module.f90.

◆ factorial

integer, dimension(qmax+1), private gear_module::factorial

private

Definition at line 65 of file gear_module.f90.

◆ gear_nr_count

integer, public gear_module::gear_nr_count = 0

counter for NR steps

Definition at line 68 of file gear_module.f90.

◆ h

real(r_kind), private gear_module::h

private

Definition at line 72 of file gear_module.f90.

◆ helperarr

integer, dimension(qmax+1), private gear_module::helperarr = (/ (ifac, ifac=1,qmax+1) /)

private

Definition at line 64 of file gear_module.f90.

◆ histsize

integer, parameter, private gear_module::histsize = 13

private

length of history in tau, e, ...

Definition at line 45 of file gear_module.f90.

◆ hq

real(r_kind), private gear_module::hq

private

Definition at line 96 of file gear_module.f90.

◆ hqm1

real(r_kind), private gear_module::hqm1

private

Definition at line 96 of file gear_module.f90.

◆ hqp1

real(r_kind), private gear_module::hqp1

private

Definition at line 96 of file gear_module.f90.

◆ ifac

integer, private gear_module::ifac

private

Definition at line 63 of file gear_module.f90.

◆ l

real(r_kind), dimension(histsize), private gear_module::l

private

Definition at line 92 of file gear_module.f90.

◆ n

integer, private gear_module::n

private

Definition at line 71 of file gear_module.f90.

◆ n_ini

integer, private gear_module::n_ini = 1

private

iteration step

Definition at line 76 of file gear_module.f90.

◆ nq

integer, private gear_module::nq

private

Definition at line 71 of file gear_module.f90.

◆ nq_ini

integer, private gear_module::nq_ini = 1

private

iteration step at current order

Definition at line 77 of file gear_module.f90.

◆ q

integer, private gear_module::q

private

Definition at line 71 of file gear_module.f90.

◆ q_ini

integer, private gear_module::q_ini = 1

private

start with order 1

Definition at line 75 of file gear_module.f90.

◆ qmax

integer, parameter, private gear_module::qmax = 5

private

maximum order of q

Definition at line 46 of file gear_module.f90.

◆ tau

real(r_kind), dimension(histsize), private gear_module::tau

private

hlist inverted

Definition at line 80 of file gear_module.f90.

◆ ti

real(r_kind), private gear_module::ti

private

Definition at line 72 of file gear_module.f90.

◆ xi

real(r_kind), dimension(histsize), private gear_module::xi

private

Definition at line 89 of file gear_module.f90.

◆ z

real(r_kind), dimension(:,:), allocatable, private gear_module::z

private

Nordsieck vector It is defined as.

\[ \vec{z}_n = \left[\vec{y}_n,h\dot{\vec{y}}_n,\frac{h^2 \ddot{\vec{y}}_n}{2!} ,...,\frac{h^q \vec{y}_n^{(q)}}{q!} \right] \]

Definition at line 81 of file gear_module.f90.

◆ znew

real(r_kind), dimension(:,:), allocatable, private gear_module::znew

private

Predictor step.

Definition at line 88 of file gear_module.f90.

Functions/Subroutines

Variables

Detailed Description

Function/Subroutine Documentation

◆ cfun()

◆ errorq()

◆ errorqm1()

◆ errorqp1()

◆ get_cfactor()

◆ get_l1()

◆ get_predictor_dydt()

◆ get_predictor_y()

◆ get_solution()

◆ get_solution_at()

◆ get_time()

◆ get_timestep()

◆ geterrors()

◆ init_dydt()

◆ init_gear_solver()

◆ nordsieck_product()

◆ prepare_next_step()

◆ qfun()

◆ revert_timestep()

◆ set_l()

◆ set_new_result()

◆ set_timestep()

◆ set_xi()

◆ shift_tau()

◆ shiftorder()

Variable Documentation

◆ alloc_stat

◆ amat

◆ e

◆ el

◆ en

◆ errq

◆ errqm1

◆ errqp1

◆ factorial

◆ gear_nr_count

◆ h

◆ helperarr

◆ histsize

◆ hq

◆ hqm1

◆ hqp1

◆ ifac

◆ l

◆ n

◆ n_ini

◆ nq

◆ nq_ini

◆ q

◆ q_ini

◆ qmax

◆ tau

◆ ti

◆ xi

◆ z

◆ znew