diff --git a/.zenodo.json b/.zenodo.json index 168f5f8..aa9d768 100644 --- a/.zenodo.json +++ b/.zenodo.json @@ -1,5 +1,5 @@ { - "description": "RK-Opt: Optimization of Runge-Kutta and related time integration methods", + "description": "RK-Opt: A Package for the Design of Numerical ODE Solvers", "license": "bsd-license", "title": "ketch/RK-Opt", "upload_type": "software", @@ -21,7 +21,7 @@ }, { "affiliation": "Capital One", - "name": "aron@ahmadia.net", + "name": "Aron Ahmadia", "orcid": "0000-0002-2573-2481" } ], diff --git a/README.md b/README.md index 48c8c99..7820260 100644 --- a/README.md +++ b/README.md @@ -1,29 +1,34 @@ -# RK-opt +# RK-opt: A Package for the Design of Numerical ODE solvers See the full documentation [here](http://numerics.kaust.edu.sa/RK-opt). -RK-opt is a collection of MATLAB code for designing optimized Runge-Kutta methods. +RK-opt is a collection of MATLAB codes for designing optimized numerical ODE solvers. +The main emphasis is on Runge-Kutta methods, but some routines deal with other classes of methods. It is primarily developed and used by the [KAUST Numerical Mathematics Group](http://numerics.kaust.edu.sa). It includes the following sub-packages: + - **polyopt**: Find optimal stability polynomials of a given degree and order of + accuracy for a specified spectrum. - **RK-coeff-opt**: Find optimal Runge-Kutta method coefficients, for a prescribed order of accuracy and number of stages. - **am_rad-opt**: Find stability functions with optimal radius of absolute monotonicity. Includes capabilities for both multistep and multistage methods. - - **polyopt**: Find optimal stability polynomials of a given degree and order of - accuracy for a specified spectrum. - **RKtools**: A collection of routines for analyzing or computing various properties of Runge-Kutta methods. For a much more extensive package along these lines, see [NodePy](http://numerics.kaust.edu.sa/nodepy). +A common workflow for designing Runge-Kutta methods is to use **polyopt** to find an +appropriate stability function and then **RK-coeff-opt** to determine the Runge-Kutta +method coefficients. # Authors -The following people have contributed code to RK-opt (listed alphabetically): +The code is primarily developed and maintained by David Ketcheson. +The following people have also made important constributions to RK-opt (listed alphabetically): - - Aron Ahmadia: **polyopt** routines - - David Ketcheson: Principal author and maintainer - - Matteo Parsani: **RK-coeff-opt** routines and organization + - Aron Ahmadia: Co-developer of **polyopt** algorithm and routines. - Zachary Grant: Extension of order conditions to multistep RK with more than two stages and addition of order conditions for orders 9-11. + - Matteo Parsani: Many improvements to **RK-coeff-opt** routines and organization. + - Hendrik Ranocha: General improvements and updates, including updating the test routines. diff --git a/lgo/LGOMAIN.F90 b/lgo/LGOMAIN.F90 deleted file mode 100755 index 7ccb6a5..0000000 --- a/lgo/LGOMAIN.F90 +++ /dev/null @@ -1,579 +0,0 @@ - -! SSP model formulation received from David Ketcheson, KAUST, October 16, 2012 - - module globals - implicit none - integer, parameter :: s=5 - integer, parameter :: p=4 - character, parameter :: RKclass*3 = 'erk' - real*8, allocatable :: order_conditions(:) - - contains - subroutine alloc_oc(n_order_conditions) - integer n_order_conditions - - allocate(order_conditions(n_order_conditions)) - end subroutine - - end module globals - - program main -! LGOMAIN.F90 --- Main program file for SSP optimization - -! Basic program structure: -! LGOMAIN calls LGO_RUN; the latter calls USER_FCT (in LGOFCT.FOR) -! iteratively. -! ------------------------------------------------------------------ - -! This main program includes the explicit definition of all LGO I/O -! arguments. The results are optionally written also to output files. - -! The detailed example below illustrates the input information setup. - -! --- Declarations --- - - use globals - - IMPLICIT REAL*8 (a-h,o-z) -! Maximal currently handled model size: maxdim variables, maxcon con - PARAMETER (maxdim=5000,maxcon=2000) - CHARACTER*20 modname, varname(maxdim), conname(maxcon), objname - DIMENSION varlb(maxdim), varnom(maxdim), varub(maxdim) - DIMENSION varopt(maxdim), conopt(maxcon) - INTEGER opmode, g_maxfct, tlimit, sstat, ctype(maxcon) - REAL*8 kt_tol, l_target - common / modsize/ nvars, ncons -! CHARACTER*40 sst(100),mst(100) - - integer ncons, nvars - - integer n_consistency_conditions, n_order_conditions - integer n_am_conditions - -! --- Model descriptors --- - -! Model name (character) - modname='SSP RK' - - -! Number of variables (positive integer) - call set_nvars(nvars) - -! Number of general constraints (non-negative integer) - call set_ncons(n_order_conditions,n_consistency_conditions,n_am_conditions) - neqcons = n_order_conditions + n_consistency_conditions - ncons=neqcons + n_am_conditions - - - call alloc_oc(n_order_conditions) -! Default variable bound, set for variables without explicit bounds - defbnd=5.d0 - - varlb(:) = 0.d0 - varub(:) = 0.d0 - varnom(:) = 0.d0 - varopt(:) = 0.d0 - conopt(:) = 0.d0 - ctype(:) = 0.d0 - varname(:) = '' - conname(:) = '' - -! Variable lower bounds (real) - do i=1,nvars - varlb(i)=0.d0 - end do - -! Variable upper bounds (real) - do i=1,nvars-1 - varub(i)=1.d0 - end do - varub(nvars) = s-p+1.d0 - -! Variable nominal values (real) - do i=1,nvars - varnom(i)=0.5d0*(varlb(i)+varub(i)) -! varnom(i)=dble(i) -! if (i.eq.nvars) varnom(i)=1. - end do - -! Objective function name (character) - objname='r' - - do i=1,nvars - varname(i)='' - end do - -! Constraint names (character) - do i=1,ncons - conname(i)='' - end do - -! Constraint types: equations 0, <= inequalities -1 - do i=1,neqcons - ctype(i)=0 - end do - - do i=neqcons+1, ncons - ctype(i)=-1 - end do - -! write(*,*) n_order_conditions, n_consistency_conditions, n_am_conditions -! write(*,*) nvars, ncons, ctype -! stop - -! --- Solver options and parameters --- - -! Operational mode selected: 0: LS; 1: B&B+LS; 2: GARS+LS; 3: MS+LS -! Suggested default for global search: 3; for local search only: 0 -! opmode=0 -! opmode=1 -! opmode=2 - opmode=3 - -! Maximal no. of model function evaluations in global search phase (integer) -! Suggested optional default settings - !g_maxfct=100 - g_maxfct=10000*(nvars+ncons+1)**2 - -! Maximal no. of model function evaluations in each executed local search -! phase (integer) -! Suggested optional default setting - l_maxfct=10000*(nvars+ncons+1)**2 - -! Maximal no. of model function evaluations in global search phase, w/o improvement (integer) -! Suggested default settings: g_maxfct/10 <= max_nosuc <= g_maxfct -! Suggested optional default setting - max_nosuc=5000*(nvars+ncons+1)**2 - -! Constraint penalty multiplier (positive real, used in merit function evaluation) -! Suggested default setting: 1.; one can incresase it to enforce feasibility - pen_mult=100.d0 - -! Target objective function value in global search phase (real) -! Suggested default setting: "close" to optimum value when known -! Otherwise can be set to an "unrealistic" value such as e.g. -1.d30 - g_target=-8. - -! Target objective function value in local search phase (real) -! Suggested default setting: "even closer" to optimum value when known -! l_target <= g_target settings are expected -! Otherwise can be set to an "unrealistic" value such as e.g. -1.d20 - l_target=-9. - -! Required local search precision parameter (positive real) -! Suggested default setting: 1.d-8 - fi_tol=1.d-8 - -! Constraint violation tolerance in final solution (positive real) -! Suggested default setting: 1.d-8 -! con_tol=1.d-8 -! in this model we want high precision to avoid the many pseudo-solutions - con_tol=1.d-8 - -! Kuhn-Tucker local optimality conditions tolerance in local search phase -! (positive real) -! Suggested default setting: 1.d-8 - kt_tol=1.d-8 - -! Seed value for built-in random number generator (non-negative integer) -! Suggested default setting: 0 - irngs=0 - -! Program execution time limit (seconds, positive integer) -! Suggested default setting: 300, but large-size and/or difficult models -! may require longer runs - tlimit=3000 - -! Optional result and logfile printout level -! 0: no printout of results -! 1: summary printout to file LGO.SUM (unit 10) -! 2: additional more detailed printout to file LGO.OUT (unit 9) -! 3: additional printout in all iterations: input argument and return -! function values, to file LGO.LOG (unit 11) -! Suggested default setting: 1; use 2 for more information, or 3 in -! checking and troubleshooting sessions - iprl=1 - -! --- LGO solver call --- - - write(*,*) nvars, ncons - CALL lgo_run (modname, nvars, ncons, varname, objname, conname, & - ctype, defbnd, varlb, varnom, varub, opmode, g_maxfct, l_maxfct, & - max_nosuc, pen_mult, g_target, l_target, fi_tol, con_tol, kt_tol, & - irngs, tlimit, iprl, varopt, conopt, objopt, resmax, nevals, & - runtime, sstat, mstat) - -! --- LGO return values and status indicators --- - -! varopt: optimized variable values (real array) -! conopt: constraint function values at varopt (real array) -! objopt: objective function value at varopt (real) -! resmax: maximal constraint violation at varopt (real) -! nevals: total no. of model function evaluations (integer) -! runtime: LGO runtime, given to 0.01 sec precision (real) - -! sstat: solver status (integer) -! Solver status return values -! sst(1)='Normal completion' ! Solver terminated with (expectedly) good result -! sst(2)='Iteration interrupt' ! Interrupt by user -! sst(3)='Resource interrupt' ! Preset no. of FCT evals or time limit (both set by user exceeded) -! sst(4)='Run terminated by solver' ! Due to set solver parameters or solver limitations -! sst(7)='Licensing error' ! Built-in size limitations are exceeded by the user's model - -! mstat: model status (integer) -! Model status return values -! mst(1)='Global search based feasible solution' ! Solution is generated by global+local search -! mst(2)='Local search based feasible solution' ! Solution is generated by local search only -! mst(3)='Unbounded solution (error)' ! In presence of variable bounds, this implies model error -! mst(4)='Infeasible solution' ! Solution found does not meet the set constraint violation tolerance -! mst(7)='Intermediate solution' ! This value is returned when sstat=2, 3, 4, or 7 - -! Verification printout to file -! This can (and typically will) be omitted in LGO solver versions -! linked to other modeling environments: AIMMS, GAMS, AMPL, MPL, ... - OPEN(unit=12, file="LGOresults", status="unknown") - n=nvars - m=ncons - WRITE (12,250) modname,nevals,objname,objopt, & - (i,varname(i),varopt(i),i=1,n) -250 FORMAT (/,5x,'--- LGO Solver Results Summary ---',//,5x, & - 'Model name ',a20,//,5x,'Number of function evaluations', & - 15x,i15,//,5x,'Objective function name ',a20,//,5x, & - 'Objective function value',16x,f20.10,//,5x, & - 'Variables',/,(5x,i9,3x,a20,8x,f20.10)) - IF (m.gt.0) THEN - WRITE (12,270) (j,conname(j),conopt(j),j=1,m) -270 FORMAT (/,5x,'Constraints', & - /,(5x,i9,3x,a20,8x,f20.10)) - WRITE (12,290) resmax -290 FORMAT (/,5x,'Maximal constraint violation',12x,f20.10) - END IF -! WRITE( 9,103) sstat, sst(sstat), mstat, mst(mstat) -!103 FORMAT(/,5x,'Solver status indicator value ',i1, 1x, a40, & -! //,5x,'Model status indicator value ',i1, 1x, a40) - WRITE(12,300) sstat, mstat -300 FORMAT(/,5x,'Solver status indicator value ',i1, & - //,5x,'Model status indicator value ',i1) - WRITE (12,310) runtime -310 FORMAT (/,5x,'LGO solver system execution time (seconds)',3x, & - f15.2) - CLOSE (12) - - STOP - end program main - -subroutine set_ncons(n_order_conditions,n_consistency_conditions, n_am_conditions) - - use globals - implicit none - integer n_order_conditions, n_consistency_conditions, n_am_conditions - - ! Number of order conditions - select case(p) - case(1) - n_order_conditions = 1 - case(2) - n_order_conditions = 2 - case(3) - n_order_conditions = 4 - case(4) - n_order_conditions = 8 - case(5) - n_order_conditions = 17 - case(6) - n_order_conditions = 37 - case(7) - n_order_conditions = 85 - end select - - n_consistency_conditions = s - - n_am_conditions = (s+1)*(s+2) - -end subroutine set_ncons - - -subroutine USER_FCT(x, obj, con) - use globals - implicit none - - real*8 x(*) - real*8 obj - real*8 con(*) - real*8 A(s,s), b(s), c(s) - real*8 K(s+1,s+1), kiprkinv(s+1,s+1), Kpow(s+1,s+1) - real*8 stage_consistency(s) - real*8 r - real*8 am_conditions_1((s+1)**2), am_conditions_2(s+1) - integer ncons - - integer n_consistency_conditions, n_order_conditions - integer n_am_conditions - - integer nvars, i - - call set_nvars(nvars) - call set_ncons(n_order_conditions,n_consistency_conditions,n_am_conditions) - - - ncons = n_order_conditions + n_consistency_conditions + n_am_conditions - - call unpack_rk(x,A,b,c) - - r = x(nvars) - obj = -r - - ! Compute order conditions - call oc_butcher(A,b,c,order_conditions) - - ! Compute consistency conditions - stage_consistency = c - sum(A,2) - - ! Compute abs. mon. conditions - K(:,:) = 0.d0 - K(1:s,1:s) = A - K(s+1,1:s) = b - - kiprkinv(:,:) = 0.d0 - Kpow(:,:) = 0.d0 - ! (I + rK)^(-1) = I - rK + (rK)^2 - (rK^3) + ... - do i = 1,s+1 - kiprkinv(i,i) = 1.d0 - Kpow(i,i) = 1.d0 - end do - do i = 1, s - Kpow = matmul(K,Kpow) - kiprkinv = kiprkinv + (-r)**i*Kpow - end do - kiprkinv = matmul(K,kiprkinv) - - ! K(I+rK)^{-1} >= 0 -! am_conditions_1 = -reshape(kiprkinv, [(s+1)**2]) - am_conditions_1 = -reshape( kiprkinv, (/ (s+1)**2 /) ) - ! rK(I+rK)^{-1}e <= 1 - am_conditions_2 = r*sum(kiprkinv,2)-1.d0 - - ! Concatenate all constraints, equalities first -! con(1:ncons) = [order_conditions,stage_consistency,am_conditions_1, & -! am_conditions_2] - con(1:ncons) = (/ order_conditions,stage_consistency,am_conditions_1, & - am_conditions_2 /) -end subroutine ! objective_and_constraints - - -subroutine set_nvars(n) - ! Set total number of decision variables - ! Note that (unlike the MATLAB code) we store all s - ! abscissae for explicit methods - - use globals - implicit none - integer n - - select case(RKclass) - !===================== - ! RK classes - !===================== - case ('erk') ! Explicit (A is strictly lower triangular) - n=2*s+s*(s-1)/2 + 1; - end select - !case 'irk' %Fully implicit - ! n=2*s+s^2+1; - !case 'irk5' %Implicit, p>=5 (so first row of A is zero) - ! n=2*s+s*(s-1); - !case 'dirk' %Diagonally Implicit (A is lower triangular) - ! n=2*s+s*(s+1)/2+1; - !case 'dirk5' %Diagonally Implicit, p>=5 (lower tri. and first row of A is zero) - ! n=2*s+s*(s+1)/2-1; - !case 'sdirk' %Singly Diagonally Implicit - ! n=2*s+s*(s+1)/2+1-s+1; - return -end subroutine set_nvars - -subroutine unpack_rk(X,A,b,c) -! Extracts the coefficient arrays from the vector of decision variables. -! -! The coefficients are tored in a single vector x as:: -! -! x=[c' b' A] -! -! A is stored row-by-row. -! Note that we always store all of c (unlike the MATLAB code). - - use globals - implicit none - - real*8 A(s,s) - real*8 b(s,1) - real*8 c(s,1) - real*8 X(*) - integer i - real*8 K(s+1,s+1) - - A(:,:) = 0.d0 - - ! Generate Butcher coefficients (for RK) or Shu-Osher coefficients (for lsRK) - select case(RKclass) - case ('erk') -! c = reshape( [X(1:s)], [s,1]) -! b = reshape( X(s+1:2*s) , [s,1]) - c = reshape( (/ X(1:s) /), (/ s,1 /) ) - b = reshape( X(s+1:2*s) , (/ s,1 /) ) - do i = 1, s - A(i,1:i-1)=X(2*s+1+(i-2)*(i-1)/2:2*s+i*(i-1)/2); - end do - end select - - K(:,:) = 0.d0 - K(1:s,1:s) = A - K(s+1,1:s) = b(:,1) - return -end -! case ('irk') -! c=X(1:s)'; b=X(s+1:2*s)'; -! A=reshape(X(2*s+1:2*s+s^2),s,s)'; -! case ('irk5') -! c=[0 X(1:s-1)]'; b=X(s:2*s-1)'; -! for i=2:s -! A(i,:)=X(i*s:(i+1)*s-1); -! end -! case ('dirk') -! c=X(1:s)'; b=X(s+1:2*s)'; -! for i=1:s -! A(i,1:i)=X(2*s+1+i*(i-1)/2:2*s+i*(i+1)/2); -! end -! case ('dirk5') -! c=[0 X(1:s-1)]'; b=X(s:2*s-1)'; -! for i=2:s -! A(i,1:i)=X(2*s-1+i*(i-1)/2:2*s-2+i*(i+1)/2); -! end -! case ('sdirk') -! c=X(1:s)'; b=X(s+1:2*s)'; -! A(1,1)=X(2*s+1); -! for i=2:s -! A(i,1:i-1)=X(2*s+2+(i-2)*(i-1)/2:2*s+1+i*(i-1)/2); -! A(i,i)=X(2*s+1); -! end - - -subroutine oc_butcher(A,b,c,coneq) - - use globals - real*8 A(s,s), b(s), c(s) - real*8 coneq(*) - if (p>=1) then - ! order 1 conditions: - coneq(1)=sum(b)-1.d0 - endif - - if (p>=2) then - ! order 2 conditions: - coneq(2)=dot_product(b,c)-1/2.d0 - endif - - if (p>=3) then - ! order 3 conditions: - coneq(3)=dot_product(b,matmul(A,c))-1/6.d0 - coneq(4)=dot_product(b,c**2)-1/3.d0 - endif - - if (p>=4) then - ! order 4 conditions: - coneq(5)=dot_product(b,matmul(A,c**2))-1/12.d0 - coneq(6)=dot_product(b,matmul(A,matmul(A,c)))-1/24.d0 - coneq(7)=dot_product(b,matmul(A,c)*c)-1/8.d0 - coneq(8)=dot_product(b,c**3)-1/4.d0 - endif - - if (p>=5) then - ! order 5 conditions: - coneq(9)=dot_product(b,matmul(A,c**3))-1/20.d0 - coneq(10)=dot_product(b,matmul(A,matmul(A,c**2)))-1/60.d0 - coneq(11)=dot_product(b,matmul(A,matmul(A,matmul(A,c))))-1/120.d0 - coneq(12)=dot_product(b,matmul(A,c*matmul(A,c)))-1/40.d0 - coneq(13)=dot_product(b,matmul(A,c**2)*c)-1/15.d0 - coneq(14)=dot_product(b,matmul(A,matmul(A,c))*c)-1/30.d0 - coneq(15)=dot_product(b,matmul(A,c)*c**2)-1/10.d0 - coneq(16)=dot_product(b,matmul(A,c)*matmul(A,c))-1/20.d0 - coneq(17)=dot_product(b,c**4)-1/5.d0 - endif - - if (p>=6) then - ! order 6 conditions: - coneq(18)=dot_product(b,matmul(A,c**4))-1/30.d0 - coneq(19)=dot_product(b,matmul(A,matmul(A,c**3)))-1/120.d0 - coneq(20)=dot_product(b,matmul(A,matmul(A,matmul(A,c**2))))-1/360.d0 - coneq(21)=dot_product(b,matmul(A,matmul(A,matmul(A,matmul(A,c)))))-1/720.d0 - coneq(22)=dot_product(b,matmul(A,matmul(A,c*matmul(A,c))))-1/240.d0 - coneq(23)=dot_product(b,matmul(A,c*matmul(A,c**2)))-1/90.d0 - coneq(24)=dot_product(b,matmul(A,c*matmul(A,matmul(A,c))))-1/180.d0 - coneq(25)=dot_product(b,matmul(A,c**2*matmul(A,c)))-1/60.d0 - coneq(26)=dot_product(b,matmul(A,matmul(A,c)*matmul(A,c)))-1/120.d0 - coneq(27)=dot_product(b,matmul(A,c**3)*c)-1/24.d0 - coneq(28)=dot_product(b,matmul(A,matmul(A,c**2))*c)-1/72.d0 - coneq(29)=dot_product(b,matmul(A,matmul(A,matmul(A,c)))*c)-1/144.d0 - coneq(30)=dot_product(b,matmul(A,c*matmul(A,c))*c)-1/48.d0 - coneq(31)=dot_product(b,matmul(A,c**2)*c**2)-1/18.d0 - coneq(32)=dot_product(b,matmul(A,matmul(A,c))*c**2)-1/36.d0 - coneq(33)=dot_product(b,matmul(A,c)*c**3)-1/12.d0 - coneq(34)=dot_product(b,matmul(A,c**2)*matmul(A,c))-1/36.d0 - coneq(35)=dot_product(b,matmul(A,matmul(A,c))*matmul(A,c))-1/72.d0 - coneq(36)=dot_product(b,matmul(A,c)*matmul(A,c)*c)-1/24.d0 - coneq(37)=dot_product(b,c**5)-1/6.d0 - endif - - if (p>=7) then - ! order 7 conditions: - coneq(38)=dot_product(b,matmul(A,c**5))-1/42.d0 - coneq(39)=dot_product(b,matmul(A,matmul(A,c**4)))-1/210.d0 - coneq(40)=dot_product(b,matmul(A,matmul(A,matmul(A,c**3))))-1/840.d0 - coneq(41)=dot_product(b,matmul(A,matmul(A,matmul(A,matmul(A,c**2)))))-1/2520.d0 - coneq(42)=dot_product(b,matmul(A,matmul(A,matmul(A,matmul(A,matmul(A,c))))))-1/5040.d0 - coneq(43)=dot_product(b,matmul(A,matmul(A,matmul(A,c*matmul(A,c)))))-1/1680.d0 - coneq(44)=dot_product(b,matmul(A,matmul(A,c*matmul(A,c**2))))-1/630.d0 - coneq(45)=dot_product(b,matmul(A,matmul(A,c*matmul(A,matmul(A,c)))))-1/1260.d0 - coneq(46)=dot_product(b,matmul(A,matmul(A,c**2*matmul(A,c))))-1/420.d0 - coneq(47)=dot_product(b,matmul(A,matmul(A,matmul(A,c)*matmul(A,c))))-1/840.d0 - coneq(48)=dot_product(b,matmul(A,c*matmul(A,c**3)))-1/168.d0 - coneq(49)=dot_product(b,matmul(A,c*matmul(A,matmul(A,c**2))))-1/504.d0 - coneq(50)=dot_product(b,matmul(A,c*matmul(A,matmul(A,matmul(A,c)))))-1/1008.d0 - coneq(51)=dot_product(b,matmul(A,c*matmul(A,c*matmul(A,c))))-1/336.d0 - coneq(52)=dot_product(b,matmul(A,c**2*matmul(A,c**2)))-1/126.d0 - coneq(53)=dot_product(b,matmul(A,c**2*matmul(A,matmul(A,c))))-1/252.d0 - coneq(54)=dot_product(b,matmul(A,c**3*matmul(A,c)))-1/84.d0 - coneq(55)=dot_product(b,matmul(A,matmul(A,c)*matmul(A,c**2)))-1/252.d0 - coneq(56)=dot_product(b,matmul(A,matmul(A,c)*matmul(A,matmul(A,c))))-1/504.d0 - coneq(57)=dot_product(b,matmul(A,c*matmul(A,c)*matmul(A,c)))-1/168.d0 - coneq(58)=dot_product(b,matmul(A,c**4)*c)-1/35.d0 - coneq(59)=dot_product(b,matmul(A,matmul(A,c**3))*c)-1/140.d0 - coneq(60)=dot_product(b,matmul(A,matmul(A,matmul(A,c**2)))*c)-1/420.d0 - coneq(61)=dot_product(b,matmul(A,matmul(A,matmul(A,matmul(A,c))))*c)-1/840.d0 - coneq(62)=dot_product(b,matmul(A,matmul(A,c*matmul(A,c)))*c)-1/280.d0 - coneq(63)=dot_product(b,matmul(A,c*matmul(A,c**2))*c)-1/105.d0 - coneq(64)=dot_product(b,matmul(A,c*matmul(A,matmul(A,c)))*c)-1/210.d0 - coneq(65)=dot_product(b,matmul(A,c**2*matmul(A,c))*c)-1/70.d0 - coneq(66)=dot_product(b,matmul(A,matmul(A,c)*matmul(A,c))*c)-1/140.d0 - coneq(67)=dot_product(b,matmul(A,c**3)*c**2)-1/28.d0 - coneq(68)=dot_product(b,matmul(A,matmul(A,c**2))*c**2)-1/84.d0 - coneq(69)=dot_product(b,matmul(A,matmul(A,matmul(A,c)))*c**2)-1/168.d0 - coneq(70)=dot_product(b,matmul(A,c*matmul(A,c))*c**2)-1/56.d0 - coneq(71)=dot_product(b,matmul(A,c**2)*c**3)-1/21.d0 - coneq(72)=dot_product(b,matmul(A,matmul(A,c))*c**3)-1/42.d0 - coneq(73)=dot_product(b,matmul(A,c)*c**4)-1/14.d0 - coneq(74)=dot_product(b,matmul(A,c**2)*matmul(A,c**2))-1/63.d0 - coneq(75)=dot_product(b,matmul(A,matmul(A,c))*matmul(A,c**2))-1/126.d0 - coneq(76)=dot_product(b,matmul(A,matmul(A,c))*matmul(A,matmul(A,c)))-1/252.d0 - coneq(77)=dot_product(b,matmul(A,c**3)*matmul(A,c))-1/56.d0 - coneq(78)=dot_product(b,matmul(A,matmul(A,c**2))*matmul(A,c))-1/168.d0 - coneq(79)=dot_product(b,matmul(A,matmul(A,matmul(A,c)))*matmul(A,c))-1/336.d0 - coneq(80)=dot_product(b,matmul(A,c*matmul(A,c))*matmul(A,c))-1/112.d0 - coneq(81)=dot_product(b,matmul(A,c**2)*matmul(A,c)*c)-1/42.d0 - coneq(82)=dot_product(b,matmul(A,matmul(A,c))*matmul(A,c)*c)-1/84.d0 - coneq(83)=dot_product(b,matmul(A,c)*matmul(A,c)*c**2)-1/28.d0 - coneq(84)=dot_product(b,matmul(A,c)*matmul(A,c)*matmul(A,c))-1/56.d0 - coneq(85)=dot_product(b,c**6)-1/7.d0 - endif - -return -end subroutine - diff --git a/lgo/compile.sh b/lgo/compile.sh deleted file mode 100755 index a4a6ac6..0000000 --- a/lgo/compile.sh +++ /dev/null @@ -1,5 +0,0 @@ -#!/bin/sh - -gfortran -c -Wall LGOMAIN.F90 > messages.txt -gfortran -c -Wall LGO.FOR >> messages.txt -gfortran LGOMAIN.O LGO.O -o LGO.EXE >> messages.txt diff --git a/lgo/test2.F90 b/lgo/test2.F90 deleted file mode 100644 index 7352b92..0000000 --- a/lgo/test2.F90 +++ /dev/null @@ -1,343 +0,0 @@ -module globals - implicit none - integer, parameter :: s=2 - integer, parameter :: p=2 - character, parameter :: RKclass*3 = 'erk' -end module globals - -program main - use globals - real*8 obj - real*8, allocatable :: con(:), x(:) - integer n1, n2, n3, ncons, nvars - - call set_ncons(n1,n2,n3) - ncons = n1+n2+n3 - allocate(con(ncons)) - - call set_nvars(nvars) - allocate(x(nvars)) - x = [0.d0,1.d0,0.5d0,0.5d0,1.d0,1.d0] - - !construct some dummy data and call obj_and_constraints() - - call obj_and_constraints(x,obj,con) - call obj_and_constraints(x,obj,con) - -endprogram main - - -subroutine set_ncons(n_order_conditions,n_consistency_conditions, n_am_conditions) - - use globals - implicit none - integer n_order_conditions, n_consistency_conditions, n_am_conditions - - ! Number of order conditions - select case(p) - case(1) - n_order_conditions = 1 - case(2) - n_order_conditions = 2 - case(3) - n_order_conditions = 4 - case(4) - n_order_conditions = 8 - case(5) - n_order_conditions = 17 - case(6) - n_order_conditions = 37 - case(7) - n_order_conditions = 85 - end select - - n_consistency_conditions = s - - n_am_conditions = (s+1)*(s+2) - -end subroutine set_ncons - - -subroutine obj_and_constraints(x, obj, con) - use globals - implicit none - - real*8 x(*) - real*8 obj - real*8 con(*) - real*8 A(s,s), b(s), c(s) - real*8, allocatable :: order_conditions(:) - real*8 K(s+1,s+1), kiprkinv(s+1,s+1), Kpow(s+1,s+1) - real*8 stage_consistency(s) - real*8 r - real*8 am_conditions_1((s+1)**2), am_conditions_2(s+1) - integer ncons - - integer n_consistency_conditions, n_order_conditions - integer n_am_conditions - - integer nvars, i - - call set_nvars(nvars) - call set_ncons(n_order_conditions,n_consistency_conditions,n_am_conditions) - - - ncons = n_order_conditions + n_consistency_conditions + n_am_conditions - - call unpack_rk(x,A,b,c) - - r = x(nvars) - obj = -r - - allocate(order_conditions(n_order_conditions)) - - ! Compute order conditions - call oc_butcher(A,b,c,order_conditions) - - ! Compute consistency conditions - stage_consistency = c - sum(A,2) - - ! Compute abs. mon. conditions - K(:,:) = 0.d0 - K(1:s,1:s) = A - K(s+1,1:s) = b - - kiprkinv(:,:) = 0.d0 - Kpow(:,:) = 0.d0 - ! (I + rK)^(-1) = I - rK + (rK)^2 - (rK^3) + ... - do i = 1,s+1 - kiprkinv(i,i) = 1.d0 - Kpow(i,i) = 1.d0 - end do - do i = 1, s - Kpow = matmul(K,Kpow) - kiprkinv = kiprkinv + (-r)**i*Kpow - end do - kiprkinv = matmul(K,kiprkinv) - - ! K(I+rK)^{-1} >= 0 - am_conditions_1 = -reshape(kiprkinv, [(s+1)**2]) - ! rK(I+rK)^{-1}e <= 1 - am_conditions_2 = r*sum(kiprkinv,2)-1.d0 - - ! Concatenate all constraints, equalities first - con(1:ncons) = [order_conditions,stage_consistency,am_conditions_1, & - am_conditions_2] - -end subroutine obj_and_constraints - - -subroutine set_nvars(n) - ! Set total number of decision variables - ! Note that (unlike the MATLAB code) we store all s - ! abscissae for explicit methods - - use globals - implicit none - integer n - - select case(RKclass) - !===================== - ! RK classes - !===================== - case ('erk') ! Explicit (A is strictly lower triangular) - n=2*s+s*(s-1)/2 + 1; - end select - !case 'irk' %Fully implicit - ! n=2*s+s^2+1; - !case 'irk5' %Implicit, p>=5 (so first row of A is zero) - ! n=2*s+s*(s-1); - !case 'dirk' %Diagonally Implicit (A is lower triangular) - ! n=2*s+s*(s+1)/2+1; - !case 'dirk5' %Diagonally Implicit, p>=5 (lower tri. and first row of A is zero) - ! n=2*s+s*(s+1)/2-1; - !case 'sdirk' %Singly Diagonally Implicit - ! n=2*s+s*(s+1)/2+1-s+1; - return -end subroutine set_nvars - -subroutine unpack_rk(X,A,b,c) -! Extracts the coefficient arrays from the vector of decision variables. -! -! The coefficients are tored in a single vector x as:: -! -! x=[c' b' A] -! -! A is stored row-by-row. -! Note that we always store all of c (unlike the MATLAB code). - - use globals - implicit none - - real*8 A(s,s) - real*8 b(s,1) - real*8 c(s,1) - real*8 X(*) - integer i - real*8 K(s+1,s+1) - - A(:,:) = 0.d0 - - ! Generate Butcher coefficients (for RK) or Shu-Osher coefficients (for lsRK) - select case(RKclass) - case ('erk') - c = reshape( [X(1:s)], [s,1]) - b = reshape( X(s+1:2*s) , [s,1]) - do i = 1, s - A(i,1:i-1)=X(2*s+1+(i-2)*(i-1)/2:2*s+i*(i-1)/2); - end do - end select - - K(:,:) = 0.d0 - K(1:s,1:s) = A - K(s+1,1:s) = b(:,1) - return -end -! case ('irk') -! c=X(1:s)'; b=X(s+1:2*s)'; -! A=reshape(X(2*s+1:2*s+s^2),s,s)'; -! case ('irk5') -! c=[0 X(1:s-1)]'; b=X(s:2*s-1)'; -! for i=2:s -! A(i,:)=X(i*s:(i+1)*s-1); -! end -! case ('dirk') -! c=X(1:s)'; b=X(s+1:2*s)'; -! for i=1:s -! A(i,1:i)=X(2*s+1+i*(i-1)/2:2*s+i*(i+1)/2); -! end -! case ('dirk5') -! c=[0 X(1:s-1)]'; b=X(s:2*s-1)'; -! for i=2:s -! A(i,1:i)=X(2*s-1+i*(i-1)/2:2*s-2+i*(i+1)/2); -! end -! case ('sdirk') -! c=X(1:s)'; b=X(s+1:2*s)'; -! A(1,1)=X(2*s+1); -! for i=2:s -! A(i,1:i-1)=X(2*s+2+(i-2)*(i-1)/2:2*s+1+i*(i-1)/2); -! A(i,i)=X(2*s+1); -! end - - -subroutine oc_butcher(A,b,c,coneq) - - use globals - real*8 A(s,s), b(s), c(s) - real*8 coneq(*) - - if (p>=1) then - ! order 1 conditions: - coneq(1)=sum(b)-1d0 - endif - - if (p>=2) then - ! order 2 conditions: - coneq(2)=dot_product(b,c)-1/2d0 - endif - - if (p>=3) then - ! order 3 conditions: - coneq(3)=dot_product(b,matmul(A,c))-1/6d0 - coneq(4)=dot_product(b,c**2)-1/3d0 - endif - - if (p>=4) then - ! order 4 conditions: - coneq(5)=dot_product(b,matmul(A,c**2))-1/12d0 - coneq(6)=dot_product(b,matmul(A,matmul(A,c)))-1/24d0 - coneq(7)=dot_product(b,matmul(A,c)*c)-1/8d0 - coneq(8)=dot_product(b,c**3)-1/4d0 - endif - - if (p>=5) then - ! order 5 conditions: - coneq(9)=dot_product(b,matmul(A,c**3))-1/20d0 - coneq(10)=dot_product(b,matmul(A,matmul(A,c**2)))-1/60d0 - coneq(11)=dot_product(b,matmul(A,matmul(A,matmul(A,c))))-1/120d0 - coneq(12)=dot_product(b,matmul(A,c*matmul(A,c)))-1/40d0 - coneq(13)=dot_product(b,matmul(A,c**2)*c)-1/15d0 - coneq(14)=dot_product(b,matmul(A,matmul(A,c))*c)-1/30d0 - coneq(15)=dot_product(b,matmul(A,c)*c**2)-1/10d0 - coneq(16)=dot_product(b,matmul(A,c)*matmul(A,c))-1/20d0 - coneq(17)=dot_product(b,c**4)-1/5d0 - endif - - if (p>=6) then - ! order 6 conditions: - coneq(18)=dot_product(b,matmul(A,c**4))-1/30d0 - coneq(19)=dot_product(b,matmul(A,matmul(A,c**3)))-1/120d0 - coneq(20)=dot_product(b,matmul(A,matmul(A,matmul(A,c**2))))-1/360d0 - coneq(21)=dot_product(b,matmul(A,matmul(A,matmul(A,matmul(A,c)))))-1/720d0 - coneq(22)=dot_product(b,matmul(A,matmul(A,c*matmul(A,c))))-1/240d0 - coneq(23)=dot_product(b,matmul(A,c*matmul(A,c**2)))-1/90d0 - coneq(24)=dot_product(b,matmul(A,c*matmul(A,matmul(A,c))))-1/180d0 - coneq(25)=dot_product(b,matmul(A,c**2*matmul(A,c)))-1/60d0 - coneq(26)=dot_product(b,matmul(A,matmul(A,c)*matmul(A,c)))-1/120d0 - coneq(27)=dot_product(b,matmul(A,c**3)*c)-1/24d0 - coneq(28)=dot_product(b,matmul(A,matmul(A,c**2))*c)-1/72d0 - coneq(29)=dot_product(b,matmul(A,matmul(A,matmul(A,c)))*c)-1/144d0 - coneq(30)=dot_product(b,matmul(A,c*matmul(A,c))*c)-1/48d0 - coneq(31)=dot_product(b,matmul(A,c**2)*c**2)-1/18d0 - coneq(32)=dot_product(b,matmul(A,matmul(A,c))*c**2)-1/36d0 - coneq(33)=dot_product(b,matmul(A,c)*c**3)-1/12d0 - coneq(34)=dot_product(b,matmul(A,c**2)*matmul(A,c))-1/36d0 - coneq(35)=dot_product(b,matmul(A,matmul(A,c))*matmul(A,c))-1/72d0 - coneq(36)=dot_product(b,matmul(A,c)*matmul(A,c)*c)-1/24d0 - coneq(37)=dot_product(b,c**5)-1/6d0 - endif - - if (p>=7) then - ! order 7 conditions: - coneq(38)=dot_product(b,matmul(A,c**5))-1/42d0 - coneq(39)=dot_product(b,matmul(A,matmul(A,c**4)))-1/210d0 - coneq(40)=dot_product(b,matmul(A,matmul(A,matmul(A,c**3))))-1/840d0 - coneq(41)=dot_product(b,matmul(A,matmul(A,matmul(A,matmul(A,c**2)))))-1/2520d0 - coneq(42)=dot_product(b,matmul(A,matmul(A,matmul(A,matmul(A,matmul(A,c))))))-1/5040d0 - coneq(43)=dot_product(b,matmul(A,matmul(A,matmul(A,c*matmul(A,c)))))-1/1680d0 - coneq(44)=dot_product(b,matmul(A,matmul(A,c*matmul(A,c**2))))-1/630d0 - coneq(45)=dot_product(b,matmul(A,matmul(A,c*matmul(A,matmul(A,c)))))-1/1260d0 - coneq(46)=dot_product(b,matmul(A,matmul(A,c**2*matmul(A,c))))-1/420d0 - coneq(47)=dot_product(b,matmul(A,matmul(A,matmul(A,c)*matmul(A,c))))-1/840d0 - coneq(48)=dot_product(b,matmul(A,c*matmul(A,c**3)))-1/168d0 - coneq(49)=dot_product(b,matmul(A,c*matmul(A,matmul(A,c**2))))-1/504d0 - coneq(50)=dot_product(b,matmul(A,c*matmul(A,matmul(A,matmul(A,c)))))-1/1008d0 - coneq(51)=dot_product(b,matmul(A,c*matmul(A,c*matmul(A,c))))-1/336d0 - coneq(52)=dot_product(b,matmul(A,c**2*matmul(A,c**2)))-1/126d0 - coneq(53)=dot_product(b,matmul(A,c**2*matmul(A,matmul(A,c))))-1/252d0 - coneq(54)=dot_product(b,matmul(A,c**3*matmul(A,c)))-1/84d0 - coneq(55)=dot_product(b,matmul(A,matmul(A,c)*matmul(A,c**2)))-1/252d0 - coneq(56)=dot_product(b,matmul(A,matmul(A,c)*matmul(A,matmul(A,c))))-1/504d0 - coneq(57)=dot_product(b,matmul(A,c*matmul(A,c)*matmul(A,c)))-1/168d0 - coneq(58)=dot_product(b,matmul(A,c**4)*c)-1/35d0 - coneq(59)=dot_product(b,matmul(A,matmul(A,c**3))*c)-1/140d0 - coneq(60)=dot_product(b,matmul(A,matmul(A,matmul(A,c**2)))*c)-1/420d0 - coneq(61)=dot_product(b,matmul(A,matmul(A,matmul(A,matmul(A,c))))*c)-1/840d0 - coneq(62)=dot_product(b,matmul(A,matmul(A,c*matmul(A,c)))*c)-1/280d0 - coneq(63)=dot_product(b,matmul(A,c*matmul(A,c**2))*c)-1/105d0 - coneq(64)=dot_product(b,matmul(A,c*matmul(A,matmul(A,c)))*c)-1/210d0 - coneq(65)=dot_product(b,matmul(A,c**2*matmul(A,c))*c)-1/70d0 - coneq(66)=dot_product(b,matmul(A,matmul(A,c)*matmul(A,c))*c)-1/140d0 - coneq(67)=dot_product(b,matmul(A,c**3)*c**2)-1/28d0 - coneq(68)=dot_product(b,matmul(A,matmul(A,c**2))*c**2)-1/84d0 - coneq(69)=dot_product(b,matmul(A,matmul(A,matmul(A,c)))*c**2)-1/168d0 - coneq(70)=dot_product(b,matmul(A,c*matmul(A,c))*c**2)-1/56d0 - coneq(71)=dot_product(b,matmul(A,c**2)*c**3)-1/21d0 - coneq(72)=dot_product(b,matmul(A,matmul(A,c))*c**3)-1/42d0 - coneq(73)=dot_product(b,matmul(A,c)*c**4)-1/14d0 - coneq(74)=dot_product(b,matmul(A,c**2)*matmul(A,c**2))-1/63d0 - coneq(75)=dot_product(b,matmul(A,matmul(A,c))*matmul(A,c**2))-1/126d0 - coneq(76)=dot_product(b,matmul(A,matmul(A,c))*matmul(A,matmul(A,c)))-1/252d0 - coneq(77)=dot_product(b,matmul(A,c**3)*matmul(A,c))-1/56d0 - coneq(78)=dot_product(b,matmul(A,matmul(A,c**2))*matmul(A,c))-1/168d0 - coneq(79)=dot_product(b,matmul(A,matmul(A,matmul(A,c)))*matmul(A,c))-1/336d0 - coneq(80)=dot_product(b,matmul(A,c*matmul(A,c))*matmul(A,c))-1/112d0 - coneq(81)=dot_product(b,matmul(A,c**2)*matmul(A,c)*c)-1/42d0 - coneq(82)=dot_product(b,matmul(A,matmul(A,c))*matmul(A,c)*c)-1/84d0 - coneq(83)=dot_product(b,matmul(A,c)*matmul(A,c)*c**2)-1/28d0 - coneq(84)=dot_product(b,matmul(A,c)*matmul(A,c)*matmul(A,c))-1/56d0 - coneq(85)=dot_product(b,c**6)-1/7d0 - endif - -return -end subroutine diff --git a/paper.bib b/paper.bib new file mode 100644 index 0000000..51994d9 --- /dev/null +++ b/paper.bib @@ -0,0 +1,211 @@ +@article{2012_optimal_stability_polynomials, + Author = {Ketcheson, David I. and Ahmadia, Aron J.}, + Journal = {Communications in Applied Mathematics and Computational Science}, + Number = {2}, + Pages = {247--271}, + Title = {Optimal stability polynomials for numerical integration of initial value problems}, + Volume = {7}, + Year = {2012} +} + + +@article{2009_monotonicity-glms, + Author = {Ketcheson, David I.}, + Journal = {Mathematics of Computation}, + Month = sep, + Number = {267}, + Pages = {1497--1513}, + Title = {{Computation of optimal monotonicity preserving general linear methods}}, + Volume = {78}, + Year = {2009}} + + + +@inproceedings{parsani-eccomas, + Author = {Parsani, M and Ketcheson, D I and Deconinck, W}, + Booktitle = {Proceedings of the European Congress on Computational Methods in Applied Sciences and Engineering}, + Month = {September}, + Title = {Explicit Runge-Kutta schemes for the spectral difference method applied to wave propagation problems}, + Year = {2012}} + + +@inproceedings{Parsani_finnish, + Author = {Parsani, Matteo and Ketcheson, David I}, + Booktitle = {Proceedings of the 11th Finnish Mechanics Days}, + Month = {November}, + Pages = {1--5}, + Title = {{Optimized low-order explicit Runge-Kutta schemes for high- order spectral difference method}}, + Year = {2012}} + +@article{2013_sd_erk, + Author = {Parsani, Matteo and Ketcheson, David I. and Deconinck, W}, + Journal = {SIAM Journal on Scientific Computing}, + Month = jul, + Number = {2}, + Pages = {A957-A986}, + Title = {Optimized explicit {R}unge-{K}utta schemes for the spectral difference method applied to wave propagation problems}, + Volume = {35}, + Year = {2013}} + +@article{2014_ssp_rkdg, + Author = {Kubatko, Ethan J. and Yeager, Benjamin A. and Ketcheson, David I.}, + Journal = {Journal of Scientific Computing}, + Number = {2}, + Pages = {313-344}, + Title = {Optimal strong-stability-preserving {R}unge--{K}utta time discretizations for discontinuous {G}alerkin methods}, + Volume = {60}, + Year = {2014}} + +@article{conde2017implicit, + title={Implicit and implicit--explicit strong stability preserving Runge--Kutta methods with high linear order}, + author={Conde, Sidafa and Gottlieb, Sigal and Grant, Zachary J and Shadid, John N}, + journal={Journal of Scientific Computing}, + volume={73}, + number={2-3}, + pages={667--690}, + year={2017}, + publisher={Springer} +} + + +@inproceedings{2018_wso, + Author = {Ketcheson, David and Seibold, Benjamin and Shirokoff, David and Zhou, Dong}, + Booktitle = {Spectral and High Order Methods for Partial Differential Equations ICOSAHOM 2018}, + Editor = {Sherwin, Spencer et. al.w}, + Title = {DIRK Schemes with High Weak Stage Order}, + Year = {2018}} + +@article{christlieb2016explicit, + title={Explicit strong stability preserving multistage two-derivative time-stepping schemes}, + author={Christlieb, Andrew J and Gottlieb, Sigal and Grant, Zachary and Seal, David C}, + journal={Journal of Scientific Computing}, + volume={68}, + number={3}, + pages={914--942}, + year={2016}, + publisher={Springer} +} + +@article{gottlieb2015optimal, + title={Optimal explicit strong stability preserving Runge--Kutta methods with high linear order and optimal nonlinear order}, + author={Gottlieb, Sigal and Grant, Zachary and Higgs, Daniel}, + journal={Mathematics of Computation}, + volume={84}, + number={296}, + pages={2743--2761}, + year={2015} +} + +@article{grant2019strong, + title={A strong stability preserving analysis for explicit multistage two-derivative time-stepping schemes based on Taylor series conditions}, + author={Grant, Zachary and Gottlieb, Sigal and Seal, David C}, + journal={Communications on Applied Mathematics and Computation}, + volume={1}, + number={1}, + pages={21--59}, + year={2019}, + publisher={Springer} +} + +@inproceedings{reynoso2017strong, + title={Strong stability preserving sixth order two-derivative Runge--Kutta methods}, + author={Reynoso, Gustavo Franco and Gottlieb, Sigal and Grant, Zachary J}, + booktitle={AIP Conference Proceedings}, + volume={1863}, + number={1}, + pages={560068}, + year={2017}, + organization={AIP Publishing LLC} +} + +@book{hairer1993, + Address = {Berlin}, + Author = {Hairer, Ernst and N\o rsett, Syvert P. and Wanner, G.}, + Edition = {Second}, + Publisher = {Springer}, + Series = {Springer Series in Computational Mathematics}, + Title = {Solving ordinary differential equations {I}: Nonstiff Problems}, + Year = {1993} +} + +@book{Hairer:ODEs2, + Address = {Berlin}, + Author = {Hairer, Ernst and Wanner, G.}, + Edition = {Second}, + Publisher = {Springer}, + Series = {Springer Series in Computational Mathematics}, + Title = {Solving ordinary differential equations {II}: Stiff and differential-algebraic problems}, + Volume = {14}, + Year = {1996} +} + + +@article{2008_explicit_ssp, + Author = {Ketcheson, David I.}, + Journal = {SIAM Journal on Scientific Computing}, + Number = {4}, + Pages = {2113--2136}, + Title = {Highly Efficient Strong Stability Preserving {R}unge-{K}utta Methods with Low-Storage Implementations}, + Volume = {30}, + Year = {2008}} + + +@article{2009_implicit_ssp, + Author = {Ketcheson, David I. and Macdonald, Colin B. and Gottlieb, Sigal}, + Journal = {Applied Numerical Mathematics}, + Month = feb, + Number = {2}, + Pages = {373--392}, + Title = {{Optimal implicit strong stability preserving {R}unge--{K}utta methods}}, + Volume = {59}, + Year = {2009}} + + +@article{2010_LSRK, + Author = {Ketcheson, David I.}, + Journal = {Journal of Computational Physics}, + Month = mar, + Number = {5}, + Pages = {1763--1773}, + Title = {{R}unge--{K}utta methods with minimum storage implementations}, + Volume = {229}, + Year = {2010}} + + +@article{2011_dwssp, + Author = {Ketcheson, David I.}, + Journal = {SIAM Journal on Numerical Analysis}, + Month = may, + Number = {4}, + Pages = {1649}, + Title = {{Step Sizes for Strong Stability Preservation with Downwind-Biased Operators}}, + Volume = {49}, + Year = {2011}} + + +@article{2011_tsrk, + Author = {Ketcheson, David I. and Gottlieb, Sigal and Macdonald, Colin B.}, + Journal = {SIAM Journal on Numerical Analysis}, + Month = jun, + Number = {6}, + Pages = {2618-2639}, + Title = {Strong Stability Preserving Two-step {R}unge--{K}utta Methods}, + Volume = {49}, + Year = {2011}} + + +@article{2017_msrk, + Author = {Bresten, Christopher and Gottlieb, Sigal and Grant, Zachary and Higgs, Daniel and Ketcheson, David I. and N{\'e}meth, Adrian}, + Journal = {Math. of Comp.}, + Pages = {747-769}, + Title = {Strong Stability Preserving Multistep {R}unge-{K}utta Methods}, + Volume = {86}, + Year = {2017}} + + +@article{conde2018embedded, + title={Embedded error estimation and adaptive step-size control for optimal explicit strong stability preserving Runge--Kutta methods}, + author={Conde, Sidafa and Fekete, Imre and Shadid, John N}, + journal={arXiv preprint arXiv:1806.08693}, + year={2018} +} diff --git a/paper.md b/paper.md new file mode 100644 index 0000000..41612ab --- /dev/null +++ b/paper.md @@ -0,0 +1,162 @@ +--- +title: '`RK-Opt`: A package for the design of numerical ODE solvers' +tags: + - Python + - numerical analysis + - differential equations + - Runge-Kutta method + - linear multistep method + - strong stability preservation + - absolute stability +authors: + - name: David I. Ketcheson^[Corresponding author.] + orcid: 0000-0002-1212-126X + affiliation: 1 + - name: Hendrik Ranocha + orcid: 0000-0002-3456-2277 + affiliation: 1 + - name: Matteo Parsani + orcid: 0000-0001-7300-1280 + affiliation: 1 + - name: Aron Ahmadia + orcid: 0000-0002-2573-2481 + affiliation: 2 +affiliations: + - name: King Abdullah University of Science and Technology + index: 1 + - name: Capital One + index: 2 +date: 9 July 2020 +bibliography: paper.bib +--- + +# Summary + +Ordinary and partial differential equations (ODEs and PDEs) are used to model +many important phenomena. In most cases, solutions of these models must be +approximated by numerical methods. Most of the relevant algorithms fall within +a few classes of methods, with the properties of individual methods determined +by their coefficients. The choice of appropriate coefficients in the design of +methods for specific applications is an important area of research. +`RK-Opt` is a software package for designing numerical ODE solvers with +coefficients optimally chosen to provide desired properties. +The original (and still primary) focus of the package is on the design of +Runge-Kutta methods, but some routines for designing other classes of methods +are also included. + +# Statement of need + +Over the last several decades, a great deal of work has gone into the design +of numerical ODE solvers. Initially this work was aimed at developing general +purpose solvers, but over time the emphasis shifted increasingly toward +development of optimized methods for specific applications. Different +accuracy, stability, performance, and other properties may be relevant or +essential depending on the nature of the equations to be solved. + +`RK-Opt` provides code that can enforce desired properties and/or objective +functions. The constraints and objective are then used within an optimization +framework, to determine coefficients of methods that best achieve the desired +goal. Thus `RK-Opt` is a sort of meta-software, consisting of algorithms whose +purpose is to create other algorithms. + +Typically, the most obvious formulation of the corresponding +optimization is intractable. Therefore, these optimization problems +are reformulated in ways that make them amenable to available techniques. +These reformulations include, for instance, turning a nonconvex problem into +a sequence of convex problems or even linear programs. The resulting algorithms +can often guarantee optimality of their output. However, for the +general problem of determining Runge-Kutta coefficients, the nonconvex problem +must be attacked directly and optimality cannot be guaranteed. + +`RK-Opt` is written entirely in MATLAB, and leverages the MATLAB Optimization +Toolbox as well as the Global Optimization Toolbox. +Its development has been motivated largely by research needs and +it has been used in a number of papers (see below). + +# Features + +`RK-Opt` includes the following subpackages. + +## `polyopt` + +This package computes optimal stability functions for Runge-Kutta methods. +Here optimal means that the stable step size is maximized for a given ODE +spectrum. The corresponding optimization problem is intractable under a +direct implementation. The package uses the algorithm developed in +[@2012_optimal_stability_polynomials], which transforms the problem into a +sequence of convex problems and typically yields a solution in a few seconds or +less. This package is usually used as the first step in designing a +Runge-Kutta method. + +## `RK-Coeff-Opt` + +This package computes optimal Runge-Kutta coefficients based on a desired +set of constraints and an objective. Available constraints include: + + - The number of stages and order of accuracy + - The class of method (explicit, implicit, diagonally implicit, low-storage) + - The coefficients of the stability polynomial (usually determined using `polyopt`) + +Two objective functions are provided; methods can be optimized for the +strong stability preserving (SSP) coefficient or the principal error norm +(a measure of the leading-order truncation error coefficients). +In addition to standard Runge-Kutta methods, various classes of multistep +Runge-Kutta methods can also be optimized. + +The optimization problem in question is highly nonconvex and the available +solvers may fail to find a solution or converge to a non-optimal solution. +For this reason, the implementation is based on doing many local optimizations +in parallel from different random initial points, using MATLAB's Global +Optimization Toolbox. + +The packages `dwrk-opt` and `low-storage` are specialized but less full-featured +versions of `RK-Coeff-Opt` that were developed for specific research projects +involving downwind Runge-Kutta methods and low-storage Runge-Kutta methods. + +## `am_radius-opt` + +Whereas the previous two subpackages are fairly general-purpose tools, +this package solves a very specific set of problems described in +[@]. Specifically, the provided routines determine the cofficients of +multistep methods (including classes of general linear methods) with the +largest possible SSP coefficient (also known +as radius of absolute monotonicity). The corresponding optimization problem +had previously been attacked using brute force search, but this limited +its solvability to methods with very few steps. In this package the +problem is reformulated as a sequence of linear programming problems, +enabling its efficient solution for methods with many steps. + + +# Related research and software + +`RK-Opt` development has proceeded in close connection to the `NodePy` package. +Whereas `RK-Opt` is focused on the design of numerical methods, `NodePy` is focused +more on their analysis. A common workflow involves generating new methods with +`RK-Opt` and then studying their properties in more detail using `NodePy`. + +Some of the research projects that have made use of `RK-Opt` include development of: + + - SSP Runge-Kutta methods + [@2008_explicit_ssp;@2009_implicit_ssp;@2013_effective_order_ssp;gottlieb2015optimal] + - SSP linear multistep methods [@2009_monotonicity] + - SSP general linear methods [@2011_tsrk;@2017_msrk] + - SSP IMEX Runge-Kutta methods [@conde2017implicit] + - Low-storage Runge-Kutta methods [@2010_LSRK] + - Optimal Runge-Kutta stability polynomials [@2012_optimal_stability_polynomials] + - Additive and downwind SSP Runge-Kutta methods [@2011_dwssp;@2018_perturbations] + - Optimal Runge-Kutta methods for specific PDE semi-discretizations [@parsani-eccomas;Parsani_finnish;2013_sd_erk;2014_ssp_rkdg] + - Embedded pairs for Runge-Kutta methods [@horvathembedded;@conde2018embedded] + - Runge-Kutta methods with high weak stage order [@2018_wso] + - SSP multistage, multiderivative methods [@christlieb2016explicit;grant2019strong;reynoso2017strong] + +As can be seen from this list, applications have mostly stemmed from the +work of the main developer's research group, but have since expanded +beyond that. + +# Acknowledgements + +Much of the initial `RK-Opt` development was performed by D. Ketcheson while +he was supported by a DOE Computational Science Graduate Fellowship. Development +has also been supported by funding from King Abdullah University of Science and Technology. + +# References