A Radau2A IRK implementation for C Adapted from Hairer and Wanner's RADAU5 code and the FATODE ODE integration library. More...

#include "header.h"
#include "solver_props.h"
#include "solver_options.h"
#include "lapack_dfns.h"
#include "dydt.h"
#include "jacob.h"
#include <complex.h>
#include <stdio.h>
#include <stdbool.h>
#include <string.h>
#include <math.h>

Include dependency graph for radau2a.c:

Go to the source code of this file.

Namespaces
	radau2a

Macros
#define	Max_no_steps (200000)
	Maximum number of allowed internal timesteps before error. More...

#define	NewtonMaxit (8)
	Maximum number of allowed Newton iteration steps before error. More...

#define	StartNewton (true)
	Use quadratic interpolation from previous step if possible. More...

#define	Gustafsson
	Use gustafsson time stepping control. More...

#define	Roundoff (EPS)
	Smallist representable double precision number. More...

#define	FacMin (0.2)
	Controls maximum decrease in timestep size. More...

#define	FacMax (8)
	Controls maximum increase in timestep size. More...

#define	FacSafe (0.9)
	Safety factor for Gustafsson time stepping control. More...

#define	FacRej (0.1)
	Time step factor on rejected step. More...

#define	ThetaMin (0.001)
	Minimum Newton convergence rate. More...

#define	NewtonTol (0.03)
	Newton convergence tolerance. More...

#define	Qmin (1.0)
	Min Timestep ratio to skip LU decomposition. More...

#define	Qmax (1.2)
	Max Timestep ratio to skip LU decomposition. More...

#define	Max_consecutive_errs (5)
	Error allowed on this many consecutive internal timesteps before exit. More...

Functions
static void	radau2a::scale (const double __restrict__ y0, const double __restrict__ y, double *__restrict__ sc)
	Computes error weight scaling from initial and current state. More...

static void	radau2a::scale_init (const double __restrict__ y0, double __restrict__ sc)
	Computes error weight scaling from initial state. More...

static void	radau2a::RK_Decomp (const double H, double __restrict__ E1, double complex __restrict__ E2, const double __restrict__ Jac, int __restrict__ ipiv1, int __restrict__ ipiv2, int __restrict__ info)
	Compute E1 & E2 matricies and their LU Decomposition. More...

static void	radau2a::RK_Make_Interpolate (const double __restrict__ Z1, const double __restrict__ Z2, const double __restrict__ Z3, double __restrict__ CONT)
	Compute Quadaratic interpolate. More...

static void	radau2a::RK_Interpolate (const double H, const double Hold, double __restrict__ Z1, double __restrict__ Z2, double __restrict__ Z3, const double __restrict__ CONT)
	Apply quadaratic interpolate to get initial values. More...

static void	radau2a::WADD (const double __restrict__ X, const double __restrict__ Y, double *__restrict__ Z)
	Performs \(Z:= X + Y\) with unrolled (or at least bounds known at compile time) loops. More...

static void	radau2a::DAXPY3 (const double DA1, const double DA2, const double DA3, const double __restrict__ DX, double __restrict__ DY1, double __restrict__ DY2, double __restrict__ DY3)
	Sepcialization of DAXPY with unrolled (or at least bounds known at compile time) loops. More...

static void	radau2a::RK_PrepareRHS (const double t, const double pr, const double H, const double __restrict__ Y, const double __restrict__ Z1, const double __restrict__ Z2, const double __restrict__ Z3, double __restrict__ R1, double __restrict__ R2, double *__restrict__ R3)
	Prepare the right-hand side for Newton iterations: \(R = Z - hA * F\). More...

static void	radau2a::RK_Solve (const double H, double __restrict__ E1, double complex __restrict__ E2, double __restrict__ R1, double __restrict__ R2, double __restrict__ R3, int __restrict__ ipiv1, int *__restrict__ ipiv2)
	Solves for the RHS values in the Newton iteration. More...

static double	radau2a::RK_ErrorNorm (const double __restrict__ scale, double __restrict__ DY)
	Computes the scaled error norm from the given `scale` and `DY` vectors. More...

static double	radau2a::RK_ErrorEstimate (const double H, const double t, const double pr, const double __restrict__ Y, const double __restrict__ F0, const double __restrict__ Z1, const double __restrict__ Z2, const double __restrict__ Z3, const double __restrict__ scale, double __restrict__ E1, int __restrict__ ipiv1, const bool FirstStep, const bool Reject)
	Computes and returns the error estimate for this step. More...

int	radau2a::integrate (const double t_start, const double t_end, const double pr, double *y)
	5th-order Radau2A CPU implementation More...

Variables
static const double	radau2a::rkA [3][3]

static const double	radau2a::rkB [3]

static const double	radau2a::rkC [3]

static const double	radau2a::rkE [4]

static const double	radau2a::rkTheta [3]

static const double	radau2a::rkGamma = 3.637834252744495732208418513577775

static const double	radau2a::rkAlpha = 2.681082873627752133895790743211112

static const double	radau2a::rkBeta = 3.050430199247410569426377624787569

static const double	radau2a::rkT [3][3]

static const double	radau2a::rkTinv [3][3]

static const double	radau2a::rkTinvAinv [3][3]

static const double	radau2a::rkAinvT [3][3]

static const double	radau2a::rkELO = 4

static char	radau2a::TRANS = 'N'
	Lapack - non-transpose. More...

static int	radau2a::NRHS = 1
	Lapack - 1 RHS solve. More...

static int	radau2a::ARRSIZE = NSP
	Lapack - Array size. More...

Detailed Description

A Radau2A IRK implementation for C Adapted from Hairer and Wanner's RADAU5 code and the FATODE ODE integration library.

Author: Nicholas J. Curtis

Date: 03/16/2015

For full reference see:
G. Wanner, E. Hairer, Solving Ordinary Differential Equations II: Stiff and DifferentialAlgebraic Problems, 2nd Edition, Springer-Verlag, Berlin, 1996. doi:10.1007/978-3-642- 05221-7.

NOTE: all matricies stored in column major format!

Definition in file radau2a.c.

Macro Definition Documentation

◆ FacMax

#define FacMax (8)

Controls maximum increase in timestep size.

Definition at line 50 of file radau2a.c.

◆ FacMin

#define FacMin (0.2)

Controls maximum decrease in timestep size.

Definition at line 48 of file radau2a.c.

◆ FacRej

#define FacRej (0.1)

Time step factor on rejected step.

Definition at line 54 of file radau2a.c.

◆ FacSafe

#define FacSafe (0.9)

Safety factor for Gustafsson time stepping control.

Definition at line 52 of file radau2a.c.

◆ Gustafsson

#define Gustafsson

Use gustafsson time stepping control.

Definition at line 44 of file radau2a.c.

◆ Max_consecutive_errs

#define Max_consecutive_errs (5)

Error allowed on this many consecutive internal timesteps before exit.

Definition at line 65 of file radau2a.c.

◆ Max_no_steps

#define Max_no_steps (200000)

Maximum number of allowed internal timesteps before error.

Definition at line 38 of file radau2a.c.

◆ NewtonMaxit

#define NewtonMaxit (8)

Maximum number of allowed Newton iteration steps before error.

Definition at line 40 of file radau2a.c.

◆ NewtonTol

#define NewtonTol (0.03)

Newton convergence tolerance.

Definition at line 58 of file radau2a.c.

◆ Qmax

#define Qmax (1.2)

Max Timestep ratio to skip LU decomposition.

Definition at line 62 of file radau2a.c.

◆ Qmin

#define Qmin (1.0)

Min Timestep ratio to skip LU decomposition.

Definition at line 60 of file radau2a.c.

◆ Roundoff

#define Roundoff (EPS)

Smallist representable double precision number.

Definition at line 46 of file radau2a.c.

◆ StartNewton

#define StartNewton (true)

Use quadratic interpolation from previous step if possible.

Definition at line 42 of file radau2a.c.

◆ ThetaMin

#define ThetaMin (0.001)

Minimum Newton convergence rate.

Definition at line 56 of file radau2a.c.

Namespaces

Macros

Functions

Variables