From 6bc74a0f3e6d54569d60a744b9ecdaf3c40c18a2 Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Thu, 13 Jun 2024 14:25:40 -0500 Subject: [PATCH 01/28] Update function prototypes to avoid compiler warnings with GCC 13 src/code.h - Add "char* rootFileName" as the argument to Use_C, Use_F, Use_F90, and Use_MATLAB function prototypes src/code_c.c - Add "char* rootFileName" as the argument to Use_C src/code_f77.c - Add "char* rootFileName" as the argument to Use_F src/code_f90.c - Add "char *rootFileName" as the argument to Use_F90 src/code_matlab.c - Add "char *rootFileName" as the argument to Use_MATLAB src/gdata.h - Add "char* rootFileName" as the argument to Generate function prototype src/gen.c - Add "char* rootFileName" as the argument to Generate Signed-off-by: Bob Yantosca --- CHANGELOG.md | 4 ++++ src/code.h | 8 ++++---- src/code_c.c | 2 +- src/code_f77.c | 2 +- src/code_f90.c | 2 +- src/code_matlab.c | 2 +- src/gdata.h | 2 +- src/gen.c | 2 +- 8 files changed, 14 insertions(+), 10 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 11fb6ca..210ea08 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,10 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [Unreleased] - TBD +### Fixed +- Add `char* rootFileName` to functions and function prototypes for `Use_C`, `Use_F`, `Use_F90`, `Use_MATLAB`, and `Generate` + ## [3.1.1] - 2024-04-30 ### Changed - Updated Python package versions for ReadTheDocs in `docs/requirements.txt` diff --git a/src/code.h b/src/code.h index 0dbaf5b..e99acfe 100644 --- a/src/code.h +++ b/src/code.h @@ -187,10 +187,10 @@ void CommentFncBegin( int f, int *vars ); void CommentFunctionBegin( int f, ... ); void CommentFunctionEnd( int f ); -void Use_C(); -void Use_F(); -void Use_F90(); -void Use_MATLAB(); +void Use_C( char* rootFileName ); +void Use_F( char* rootFileName); +void Use_F90( char* rootFileName ); +void Use_MATLAB( char* rootFileName ); extern void (*WriteElm)( NODE *n ); extern void (*WriteSymbol)( int op ); diff --git a/src/code_c.c b/src/code_c.c index 9056067..f6dc99c 100644 --- a/src/code_c.c +++ b/src/code_c.c @@ -514,7 +514,7 @@ char buf[ MAX_K ]; FlushBuf(); } -void Use_C() +void Use_C( char* rootFileName ) { WriteElm = C_WriteElm; WriteSymbol = C_WriteSymbol; diff --git a/src/code_f77.c b/src/code_f77.c index 0106f32..f0cf699 100644 --- a/src/code_f77.c +++ b/src/code_f77.c @@ -536,7 +536,7 @@ char buf[ MAX_K ]; } /*************************************************************************************************/ -void Use_F() +void Use_F( char* rootFileName ) { WriteElm = F77_WriteElm; WriteSymbol = F77_WriteSymbol; diff --git a/src/code_f90.c b/src/code_f90.c index 86b36ff..d1d7b35 100644 --- a/src/code_f90.c +++ b/src/code_f90.c @@ -750,7 +750,7 @@ char buf[ MAX_K ]; } /*************************************************************************************************/ -void Use_F90() +void Use_F90( char* rootFileName ) { // Temporary string variable char buf[ MAX_K ]; diff --git a/src/code_matlab.c b/src/code_matlab.c index c4cde27..e157abc 100644 --- a/src/code_matlab.c +++ b/src/code_matlab.c @@ -660,7 +660,7 @@ char buf[ MAX_K ]; } /*************************************************************************************************/ -void Use_MATLAB() +void Use_MATLAB( char *rootFileName ) { WriteElm = MATLAB_WriteElm; WriteSymbol = MATLAB_WriteSymbol; diff --git a/src/gdata.h b/src/gdata.h index a732f0c..333412b 100644 --- a/src/gdata.h +++ b/src/gdata.h @@ -257,7 +257,7 @@ void CmdIntegrator( char *cmd ); void CmdDriver( char *cmd ); void CmdStochastic( char *cmd ); void CmdFlux( char *cmd ); -void Generate(); +void Generate( char* rootFileName ); char * FileName( char *name, char* env, char *dir, char *ext ); diff --git a/src/gen.c b/src/gen.c index d491362..3114c26 100644 --- a/src/gen.c +++ b/src/gen.c @@ -3472,7 +3472,7 @@ case 't': /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ -void Generate() +void Generate( char* rootFileName ) { int i; int n; From 506098c9817f1507956655a892cadf048ec67f89 Mon Sep 17 00:00:00 2001 From: Simon Rosanka Date: Thu, 11 Jul 2024 12:49:48 -0700 Subject: [PATCH 02/28] Use Y instead of C in Update_RCONST - #93 Update_RCONST now uses variable concentrations for rate constants that depend on species concentrations. To ensure backwards compatibility, YIN is an optional input. If YIN is not passed to Update_RCONST, concentrations from the global C array are used instead. --- CHANGELOG.md | 3 +++ int/rosenbrock.f90 | 4 ++-- int/rosenbrock_adj.f90 | 6 +++--- int/rosenbrock_autoreduce.f90 | 4 ++-- int/rosenbrock_tlm.f90 | 6 +++--- src/gen.c | 21 +++++++++++++++++++-- 6 files changed, 32 insertions(+), 12 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 210ea08..1421ad2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,9 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). ## [Unreleased] - TBD +### Changed +- Use `Y` instead of `C` in `Update_RCONST` to account for updated variable species concentrations + ### Fixed - Add `char* rootFileName` to functions and function prototypes for `Use_C`, `Use_F`, `Use_F90`, `Use_MATLAB`, and `Generate` diff --git a/int/rosenbrock.f90 b/int/rosenbrock.f90 index f302e40..6b6f138 100644 --- a/int/rosenbrock.f90 +++ b/int/rosenbrock.f90 @@ -1326,7 +1326,7 @@ SUBROUTINE FunTemplate( T, Y, Ydot, P_VAR, D_VAR ) Told = TIME TIME = T IF ( Do_Update_SUN ) CALL Update_SUN() - IF ( Do_Update_RCONST ) CALL Update_RCONST() + IF ( Do_Update_RCONST ) CALL Update_RCONST(Y) CALL KPP_FUN_OR_FUN_SPLIT TIME = Told @@ -1364,7 +1364,7 @@ SUBROUTINE JacTemplate( T, Y, Jcb ) Told = TIME TIME = T IF ( Do_Update_SUN ) CALL Update_SUN() - IF ( Do_Update_RCONST ) CALL Update_RCONST() + IF ( Do_Update_RCONST ) CALL Update_RCONST(Y) #ifdef FULL_ALGEBRA CALL Jac_SP(Y, FIX, RCONST, JV) DO j=1,NVAR diff --git a/int/rosenbrock_adj.f90 b/int/rosenbrock_adj.f90 index b95d75a..90c6e2b 100644 --- a/int/rosenbrock_adj.f90 +++ b/int/rosenbrock_adj.f90 @@ -2523,7 +2523,7 @@ SUBROUTINE FunTemplate( T, Y, Ydot ) Told = TIME TIME = T IF ( Do_Update_SUN ) CALL Update_SUN() - IF ( Do_Update_RCONST ) CALL Update_RCONST() + IF ( Do_Update_RCONST ) CALL Update_RCONST(Y) CALL Fun( Y, FIX, RCONST, Ydot ) TIME = Told @@ -2554,7 +2554,7 @@ SUBROUTINE JacTemplate( T, Y, Jcb ) Told = TIME TIME = T IF ( Do_Update_SUN ) CALL Update_SUN() - IF ( Do_Update_RCONST ) CALL Update_RCONST() + IF ( Do_Update_RCONST ) CALL Update_RCONST(Y) #ifdef FULL_ALGEBRA CALL Jac_SP(Y, FIX, RCONST, JV) DO j=1,NVAR @@ -2591,7 +2591,7 @@ SUBROUTINE HessTemplate( T, Y, Hes ) Told = TIME TIME = T IF ( Do_Update_SUN ) CALL Update_SUN() - IF ( Do_Update_RCONST ) CALL Update_RCONST() + IF ( Do_Update_RCONST ) CALL Update_RCONST(Y) CALL Hessian( Y, FIX, RCONST, Hes ) TIME = Told diff --git a/int/rosenbrock_autoreduce.f90 b/int/rosenbrock_autoreduce.f90 index 42db25a..a7a8e07 100644 --- a/int/rosenbrock_autoreduce.f90 +++ b/int/rosenbrock_autoreduce.f90 @@ -2471,7 +2471,7 @@ SUBROUTINE FunTemplate( T, Y, Ydot, P_VAR, D_VAR ) Told = TIME TIME = T IF ( Do_Update_SUN ) CALL Update_SUN() - IF ( Do_Update_RCONST ) CALL Update_RCONST() + IF ( Do_Update_RCONST ) CALL Update_RCONST(Y) CALL KPP_FUN_OR_FUN_SPLIT TIME = Told @@ -2568,7 +2568,7 @@ SUBROUTINE JacTemplate( T, Y, Jcb ) Told = TIME TIME = T IF ( Do_Update_SUN ) CALL Update_SUN() - IF ( Do_Update_RCONST ) CALL Update_RCONST() + IF ( Do_Update_RCONST ) CALL Update_RCONST(Y) #ifdef FULL_ALGEBRA CALL Jac_SP(Y, FIX, RCONST, JV) DO j=1,NVAR diff --git a/int/rosenbrock_tlm.f90 b/int/rosenbrock_tlm.f90 index 1386503..4710de3 100644 --- a/int/rosenbrock_tlm.f90 +++ b/int/rosenbrock_tlm.f90 @@ -1348,7 +1348,7 @@ SUBROUTINE FunTemplate( T, Y, Ydot ) Told = TIME TIME = T IF ( Do_Update_SUN ) CALL Update_SUN() - IF ( Do_Update_RCONST ) CALL Update_RCONST() + IF ( Do_Update_RCONST ) CALL Update_RCONST(Y) CALL Fun( Y, FIX, RCONST, Ydot ) TIME = Told @@ -1373,7 +1373,7 @@ SUBROUTINE JacTemplate( T, Y, Jcb ) Told = TIME TIME = T IF ( Do_Update_SUN ) CALL Update_SUN() - IF ( Do_Update_RCONST ) CALL Update_RCONST() + IF ( Do_Update_RCONST ) CALL Update_RCONST(Y) CALL Jac_SP( Y, FIX, RCONST, Jcb ) TIME = Told @@ -1398,7 +1398,7 @@ SUBROUTINE HessTemplate( T, Y, Hes ) Told = TIME TIME = T IF ( Do_Update_SUN ) CALL Update_SUN() - IF ( Do_Update_RCONST ) CALL Update_RCONST() + IF ( Do_Update_RCONST ) CALL Update_RCONST(Y) CALL Hessian( Y, FIX, RCONST, Hes ) TIME = Told diff --git a/src/gen.c b/src/gen.c index 3114c26..92dd365 100644 --- a/src/gen.c +++ b/src/gen.c @@ -2167,15 +2167,32 @@ void GenerateUpdateRconst() { int i; int UPDATE_RCONST; +int YIN,Y; UseFile( rateFile ); - UPDATE_RCONST = DefFnc( "Update_RCONST", 0, "function to update rate constants"); + YIN = DefvElmO( "YIN", real, -NVAR, "Optional input concentrations of variable species" ); + UPDATE_RCONST = DefFnc( "Update_RCONST", 1, "function to update rate constants"); + + FunctionBegin( UPDATE_RCONST, YIN ); + + Y = DefvElm( "Y", real, -NSPEC, "Concentrations of species (local)" ); + Declare(Y); + NewLines(1); - FunctionBegin( UPDATE_RCONST ); F77_Inline(" INCLUDE '%s_Global.h'", rootFileName); MATLAB_Inline("global SUN TEMP RCONST"); + switch( useLang ){ + case F90_LANG: + WriteComment("Ensure local Y array is filled with variable and constant concentrations"); + bprintf(" Y(1:NSPEC) = C(1:NSPEC)\n"); + NewLines(1); + WriteComment("Update local Y array if variable concentrations are provided"); + bprintf(" if (present(YIN)) Y(1:NVAR) = YIN(1:NVAR)\n"); + break; + } + if ( useLang==F77_LANG ) IncludeCode( "%s/util/UserRateLaws_FcnHeader", Home ); From 03ff07814fc6fff1a2e9a6c61c4a4f5f88b2591f Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Wed, 17 Jul 2024 10:44:00 -0400 Subject: [PATCH 03/28] Only declare UPDATE_RCONST with optional Y argument for F90 src/gen.c - Wrap the code to declare UPDATE_RCONST with YIN in an if block that only executes when useLang==F90 - Place code to declare UPDATE_RCONST without any variables in the else block. Signed-off-by: Bob Yantosca --- src/gen.c | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/src/gen.c b/src/gen.c index 92dd365..8b728df 100644 --- a/src/gen.c +++ b/src/gen.c @@ -2171,14 +2171,20 @@ int YIN,Y; UseFile( rateFile ); - YIN = DefvElmO( "YIN", real, -NVAR, "Optional input concentrations of variable species" ); - UPDATE_RCONST = DefFnc( "Update_RCONST", 1, "function to update rate constants"); - - FunctionBegin( UPDATE_RCONST, YIN ); - - Y = DefvElm( "Y", real, -NSPEC, "Concentrations of species (local)" ); - Declare(Y); - NewLines(1); + if (useLang==F90_LANG) { + // F90: Declare function with optional YIN argument and local Y variable. + YIN = DefvElmO( "YIN", real, -NVAR, "Optional input concentrations of variable species" ); + UPDATE_RCONST = DefFnc( "Update_RCONST", 1, "function to update rate constants"); + + FunctionBegin( UPDATE_RCONST, YIN ); + Y = DefvElm( "Y", real, -NSPEC, "Concentrations of species (local)" ); + Declare(Y); + NewLines(1); + } else { + // For other languages, declare function w/o any arguments + UPDATE_RCONST = DefFnc( "Update_RCONST", 0, "function to update rate constants"); + FunctionBegin( UPDATE_RCONST ); + } F77_Inline(" INCLUDE '%s_Global.h'", rootFileName); MATLAB_Inline("global SUN TEMP RCONST"); From 45cda946044d423786eed19e215e5484035cf103 Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Wed, 17 Jul 2024 10:48:53 -0400 Subject: [PATCH 04/28] C-I tests now print a header with the compiler information .ci-pipelinoes/ci-common.defs.sh - Add a function "print_compiler_versions" to print the gcc and gfortran compiler versions that are used to run the tests. .ci-pipelines/ci-testing-script.sh - Call "print_compiler_versions" function before running tests. This will print the compiler versions at the top of stdout output. CHANGELOG.md - Updated accordingly Signed-off-by: Bob Yantosca --- .ci-pipelines/ci-common-defs.sh | 12 ++++++++++++ .ci-pipelines/ci-testing-script.sh | 3 +++ CHANGELOG.md | 1 + 3 files changed, 16 insertions(+) diff --git a/.ci-pipelines/ci-common-defs.sh b/.ci-pipelines/ci-common-defs.sh index bac2526..62340d3 100755 --- a/.ci-pipelines/ci-common-defs.sh +++ b/.ci-pipelines/ci-common-defs.sh @@ -46,6 +46,18 @@ function get_ci_test_path() { return } +# Prints a headeer with the compiler versions +function print_compiler_versions() { + echo \ +"###########################################################################" + echo " KPP CONTINUOUS INTEGRATION TESTS, USING THESE COMPILERS:" + echo "" + gcc --version + gfortran --version + echo \ +"###########################################################################" +} + # Run a C-I test function run_ci_test() { diff --git a/.ci-pipelines/ci-testing-script.sh b/.ci-pipelines/ci-testing-script.sh index 81b39e9..d8143af 100755 --- a/.ci-pipelines/ci-testing-script.sh +++ b/.ci-pipelines/ci-testing-script.sh @@ -10,6 +10,9 @@ # Current directory cwd=$(pwd -P) +# Print a header with the compiler versions +print_compiler_versions + # Run C-I tests with various mechanism + integrator combinations for this_test in ${GENERAL_TESTS}; do run_ci_test "${this_test}" "${cwd}" "" diff --git a/CHANGELOG.md b/CHANGELOG.md index 1421ad2..e2af5f0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] - TBD ### Changed - Use `Y` instead of `C` in `Update_RCONST` to account for updated variable species concentrations +- C-I tests now print the compiler version that are used ### Fixed - Add `char* rootFileName` to functions and function prototypes for `Use_C`, `Use_F`, `Use_F90`, `Use_MATLAB`, and `Generate` From 915b6f8a597429cb1b45950de73a5f99c91a983f Mon Sep 17 00:00:00 2001 From: Simon Rosanka Date: Tue, 23 Jul 2024 12:09:59 -0700 Subject: [PATCH 05/28] Update documentation related to #93 --- docs/source/using_kpp/04_input_for_kpp.rst | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/docs/source/using_kpp/04_input_for_kpp.rst b/docs/source/using_kpp/04_input_for_kpp.rst index 427b311..e235ba8 100644 --- a/docs/source/using_kpp/04_input_for_kpp.rst +++ b/docs/source/using_kpp/04_input_for_kpp.rst @@ -224,6 +224,18 @@ These reactions must be merged by adding the rate coefficients: N2O5 + H2O = 2 HNO3 : k_gas + k_aerosol; +If the reaction rate coefficient depends on the concentration of a specific +species, :code:`Y` should be referenced instead of :code:`C` in the reaction +rate coefficient declaration. This ensures that the species concentration at +the specific integration time is used when the reaction rate coefficient is +updated within the integrator. In the following example, the reaction rate +coefficient depends on the concentration of the hydrogen ion, whose +concentration at that specific integration time is given by :code:`Y(ind_Hp)`: + +.. code-block:: console + + HSO3m + HSO5m + Hp = 2 HSO4m + Hp : k_aqueous( Y(ind_Hp) ); + .. _families: #FAMILIES From d459c4e0e386fc578a65492879f6ea1975c394ac Mon Sep 17 00:00:00 2001 From: Rolf Sander Date: Thu, 25 Jul 2024 09:09:27 +0200 Subject: [PATCH 06/28] documentation adjusted --- docs/source/citations/09_acknowledgments.rst | 3 +++ docs/source/using_kpp/04_input_for_kpp.rst | 24 ++++++++++---------- 2 files changed, 15 insertions(+), 12 deletions(-) diff --git a/docs/source/citations/09_acknowledgments.rst b/docs/source/citations/09_acknowledgments.rst index 8ade7e3..69ba095 100644 --- a/docs/source/citations/09_acknowledgments.rst +++ b/docs/source/citations/09_acknowledgments.rst @@ -41,4 +41,7 @@ Stuart Lacy wrote an export function for the `Master Chemical Mechanism used out-of-the-box for the small model in the :file:`examples/mcm` directory. +Simon Rosanka improved the :code:`Update_RCONST` subroutine by providing +time-dependent concentration values via an optional parameter. + Parts of this user manual are based on :cite:t:`Damian-Iordache_1996`. diff --git a/docs/source/using_kpp/04_input_for_kpp.rst b/docs/source/using_kpp/04_input_for_kpp.rst index e235ba8..c1cdf75 100644 --- a/docs/source/using_kpp/04_input_for_kpp.rst +++ b/docs/source/using_kpp/04_input_for_kpp.rst @@ -224,18 +224,6 @@ These reactions must be merged by adding the rate coefficients: N2O5 + H2O = 2 HNO3 : k_gas + k_aerosol; -If the reaction rate coefficient depends on the concentration of a specific -species, :code:`Y` should be referenced instead of :code:`C` in the reaction -rate coefficient declaration. This ensures that the species concentration at -the specific integration time is used when the reaction rate coefficient is -updated within the integrator. In the following example, the reaction rate -coefficient depends on the concentration of the hydrogen ion, whose -concentration at that specific integration time is given by :code:`Y(ind_Hp)`: - -.. code-block:: console - - HSO3m + HSO5m + Hp = 2 HSO4m + Hp : k_aqueous( Y(ind_Hp) ); - .. _families: #FAMILIES @@ -1232,6 +1220,18 @@ ICNTRL :code:`Update_SUN` and :code:`Update_PHOTO`, but not to :code:`Update_RCONST`. + Calling :code:`Update_RCONST` may be necessary when reaction rate + coefficients depend on the concentration of a specific species, e. + g.: + + .. code-block:: console + + HSO3m + HSO5m + Hp = 2 HSO4m + Hp : k_aqueous( C(ind_Hp) ); + + This ensures that the concentration :code:`C(ind_Hp)` at the specific + integration time is used when the reaction rate coefficient is + updated within the integrator. + .. option:: ICNTRL(16) Treatment of negative concentrations: From 56242f80239a2fa5a6775a4b2268b4c8b9291fa0 Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Wed, 4 Sep 2024 10:20:59 -0400 Subject: [PATCH 07/28] Update jinja2 to 3.1.4 in docs/requirements.txt docs/requirements.txt - Use jinja 3.1.4 as there is a dependabot security warning for the previous version 3.1.3. CHANGELOG.md - Updated accordingly Signed-off-by: Bob Yantosca --- CHANGELOG.md | 7 ++++--- docs/requirements.txt | 2 +- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index e2af5f0..e4c5ef1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,11 +10,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] - TBD ### Changed -- Use `Y` instead of `C` in `Update_RCONST` to account for updated variable species concentrations -- C-I tests now print the compiler version that are used +- Updated `Update_RCONST` to use `Y` instead of `C` to account for updated variable species concentrations +- Updated C-I tests to print the compiler versions that are used ### Fixed -- Add `char* rootFileName` to functions and function prototypes for `Use_C`, `Use_F`, `Use_F90`, `Use_MATLAB`, and `Generate` +- Added `char* rootFileName` to functions and function prototypes for `Use_C`, `Use_F`, `Use_F90`, `Use_MATLAB`, and `Generate` +- Updated `docs/requirements.txt` to use `jinja2==3.1.4` (fixes a security issue) ## [3.1.1] - 2024-04-30 ### Changed diff --git a/docs/requirements.txt b/docs/requirements.txt index 60b3987..71e269d 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -12,4 +12,4 @@ sphinxcontrib-bibtex==2.6.2 sphinx-autobuild==2021.3.14 recommonmark==0.7.1 docutils==0.20.1 -jinja2==3.1.3 +jinja2==3.1.4 From ac3615980a23084d3f377097452311dd1dfac61f Mon Sep 17 00:00:00 2001 From: Rolf Sander Date: Wed, 4 Sep 2024 16:31:40 +0200 Subject: [PATCH 08/28] bug fix in error message --- int/rosenbrock.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/int/rosenbrock.f90 b/int/rosenbrock.f90 index f302e40..cf93a1c 100644 --- a/int/rosenbrock.f90 +++ b/int/rosenbrock.f90 @@ -460,7 +460,7 @@ SUBROUTINE ros_ErrorMsg(Code,T,H,IERR) CASE (-6) PRINT * , '--> No of steps exceeds maximum bound' CASE (-7) - PRINT * , '--> Step size too small: T + 10*H = T', & + PRINT * , '--> Step size too small: T + 0.1*H = T', & ' or H < Roundoff' CASE (-8) PRINT * , '--> Matrix is repeatedly singular' From a91404f07486bf0e746279355026805c56b02c52 Mon Sep 17 00:00:00 2001 From: Rolf Sander Date: Fri, 13 Sep 2024 15:01:40 +0200 Subject: [PATCH 09/28] comment needs a space after // --- docs/source/using_kpp/04_input_for_kpp.rst | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/docs/source/using_kpp/04_input_for_kpp.rst b/docs/source/using_kpp/04_input_for_kpp.rst index c1cdf75..92502a4 100644 --- a/docs/source/using_kpp/04_input_for_kpp.rst +++ b/docs/source/using_kpp/04_input_for_kpp.rst @@ -24,8 +24,9 @@ description file: :ref:`kpp-commands`, and :ref:`inlined-code`. Their syntax is presented in :ref:`bnf-description`. -- Comments are either enclosed between the curly braces :code:`{` and - :code:`}`, or written in a line starting with two slashes :code:`//`. +- Comments are either enclosed between the curly braces ":code:`{`" and + ":code:`}`", or written in a line starting with two slashes and a + space ":code:`// `". - Any name given by the user to denote an atom or a species is restricted to be less than 32 character in length and can only From 2a6b3fa2cc2b3016dcf27f530e04995ddd85fc73 Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Thu, 19 Dec 2024 16:42:04 -0500 Subject: [PATCH 10/28] Now inline "INLINED RCONST" at the top of routine "UPDATE_RCONST" src/gen.c - Modified routine GenerateUpdateRconst so that the INLINED RCONST contents will be placed immedately after the start of the UPDATE_RCONST subroutine, and to place F90 variable declarations etc. immediately after that. This will prevent compile-time errors if a USE statement is included via "INLINED RCONST", as USE statements must come before variable declarationas and other statements. - Updated comments CHANGELOG.md - Updated accordingly Signed-off-by: Bob Yantosca --- CHANGELOG.md | 1 + src/gen.c | 91 ++++++++++++++++++++++++++++++++++++++++------------ 2 files changed, 71 insertions(+), 21 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index e4c5ef1..45177ec 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -16,6 +16,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Fixed - Added `char* rootFileName` to functions and function prototypes for `Use_C`, `Use_F`, `Use_F90`, `Use_MATLAB`, and `Generate` - Updated `docs/requirements.txt` to use `jinja2==3.1.4` (fixes a security issue) +- Updated `gen.c` to write the `INLINED RCONST` section at the top of routine `UPDATE_RCONST`. This will prevent compilation errors when `INLINED_RCONST` contains F90 `USE` statements (as these must precede other F90 statements). ## [3.1.1] - 2024-04-30 ### Changed diff --git a/src/gen.c b/src/gen.c index 8b728df..7e3d172 100644 --- a/src/gen.c +++ b/src/gen.c @@ -2167,21 +2167,38 @@ void GenerateUpdateRconst() { int i; int UPDATE_RCONST; -int YIN,Y; +int YIN, Y; UseFile( rateFile ); if (useLang==F90_LANG) { - // F90: Declare function with optional YIN argument and local Y variable. - YIN = DefvElmO( "YIN", real, -NVAR, "Optional input concentrations of variable species" ); - UPDATE_RCONST = DefFnc( "Update_RCONST", 1, "function to update rate constants"); - - FunctionBegin( UPDATE_RCONST, YIN ); - Y = DefvElm( "Y", real, -NSPEC, "Concentrations of species (local)" ); - Declare(Y); + // + // SPECIAL HANDLING FOR F90: + // ------------------------- + // Write the INLINED RCONST section immediately after the start of the + // subroutine UPDATE_RCONST. This will avoid compile-time errors if a + // "USE" statement is included via INLINED RCONST. Recall that F90 + // "USE" statements must precede variable declarations or any other + // executable statements. + // + // The FunctionBegin() routine writes the variable declaration + // statements immediately after the subroutine declaration line. + // Therefore, we will not be able to use FunctionBegin() to declare + // the UPDATE_RCONST subroutine. Instead, we will manually write + // "SUBROUTINE UPDATE_RCONST ( YIN )" here. + // + // Also note, because we are not using FunctionBegin, we do not + // need to declare the "int UPDATE_RCONST;" variable if we are + // generating F90 output. Move that to the else block. + // + // -- Bob Yantosca (19 Dec 2024) + // + bprintf("SUBROUTINE Update_RCONST ( YIN )"); NewLines(1); } else { + // // For other languages, declare function w/o any arguments + // UPDATE_RCONST = DefFnc( "Update_RCONST", 0, "function to update rate constants"); FunctionBegin( UPDATE_RCONST ); } @@ -2189,16 +2206,6 @@ int YIN,Y; F77_Inline(" INCLUDE '%s_Global.h'", rootFileName); MATLAB_Inline("global SUN TEMP RCONST"); - switch( useLang ){ - case F90_LANG: - WriteComment("Ensure local Y array is filled with variable and constant concentrations"); - bprintf(" Y(1:NSPEC) = C(1:NSPEC)\n"); - NewLines(1); - WriteComment("Update local Y array if variable concentrations are provided"); - bprintf(" if (present(YIN)) Y(1:NVAR) = YIN(1:NVAR)\n"); - break; - } - if ( useLang==F77_LANG ) IncludeCode( "%s/util/UserRateLaws_FcnHeader", Home ); @@ -2223,6 +2230,35 @@ int YIN,Y; NewLines(1); WriteComment("End INLINED RCONST"); NewLines(1); + // + // SPECIAL HANDLING FOR F90: + // ------------------------- + // Write F90 variable declarations after the INLINED RCONST section. + // This will avoid compile-time errors if a USE statement is placed + // into the INLINED RCONST section, as described above. + // + // -- Bob Yantosca (19 Dec 2024) + // + if (useLang == F90_LANG) { + + // Declare optional YIN argument + YIN = DefvElmO( "YIN", real, -NVAR, "Optional input concentrations of variable species" ); + Declare(YIN); + NewLines(1); + + // Declare local Y variable + Y = DefvElm( "Y", real, -NSPEC, "Concentrations of species (local)" ); + Declare(Y); + NewLines(1); + + // Copy values of YIN to Y if YIN is present + WriteComment("Ensure local Y array is filled with variable and constant concentrations"); + bprintf(" Y(1:NSPEC) = C(1:NSPEC)\n"); + NewLines(1); + WriteComment("Update local Y array if variable concentrations are provided"); + bprintf(" if (present(YIN)) Y(1:NVAR) = YIN(1:NVAR)\n"); + NewLines(2); + } for( i = 0; i < EqnNr; i++) { /* mz_rs_20220701+ */ @@ -2242,9 +2278,22 @@ int YIN,Y; } MATLAB_Inline(" RCONST = RCONST(:);"); - - FunctionEnd( UPDATE_RCONST ); - FreeVariable( UPDATE_RCONST ); + // + // SPECIAL HANDLING FOR F90: + // ------------------------- + // Manually write the "END SUBROUTINE UPDATE_RCONST" line when + // generating F90 output. But if generating C, F77, MatLab output, + // then close the UPDATE_CONST routine as we normally would. + // -- Bob Yantosca (19 Dec 2024) + // + if (useLang == F90_LANG) { + NewLines(1); + bprintf("END SUBROUTINE UPDATE_RCONST"); + NewLines(1); + } else { + FunctionEnd( UPDATE_RCONST ); + FreeVariable( UPDATE_RCONST ); + } } From 0da69839b67bb64b8db2e11b877dc2980938bd24 Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Thu, 19 Dec 2024 17:01:21 -0500 Subject: [PATCH 11/28] Updated comments in src/gen.c src/gen.c - Removed leftover comment about "int UPDATE_RCONST;" variable - Fixed typo Signed-off-by: Bob Yantosca --- src/gen.c | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/src/gen.c b/src/gen.c index 7e3d172..b19c819 100644 --- a/src/gen.c +++ b/src/gen.c @@ -2187,10 +2187,6 @@ int YIN, Y; // the UPDATE_RCONST subroutine. Instead, we will manually write // "SUBROUTINE UPDATE_RCONST ( YIN )" here. // - // Also note, because we are not using FunctionBegin, we do not - // need to declare the "int UPDATE_RCONST;" variable if we are - // generating F90 output. Move that to the else block. - // // -- Bob Yantosca (19 Dec 2024) // bprintf("SUBROUTINE Update_RCONST ( YIN )"); @@ -2283,7 +2279,7 @@ int YIN, Y; // ------------------------- // Manually write the "END SUBROUTINE UPDATE_RCONST" line when // generating F90 output. But if generating C, F77, MatLab output, - // then close the UPDATE_CONST routine as we normally would. + // then close the UPDATE_RCONST routine as we normally would. // -- Bob Yantosca (19 Dec 2024) // if (useLang == F90_LANG) { From bd4745fff76100ecdff3ffdf7cabb4bc2f065e5c Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Fri, 20 Dec 2024 12:00:12 -0500 Subject: [PATCH 12/28] Added new #INLINE F90_RCONST_USE option for F90 USE statements src/gdata.h - Declare F90_RCONST_USE as a variable of type "enum inl_code" src/scanner.c - Added F90_RCONST_USE to InlineKeys[] src/gen.c - In routine GenerateUpdateRconst we now proceed in this order: 1. Manually write the "SUBROUTINE UPDATE_RCONST ( YIN )" line 2. Inline code from the "#INLINE F90_RCONST_USE" section 3. Declare YIN optional argument and Y local variable 4. Copy values of YIN to Y (if necessary) 5. Inline code from the "#INLINE F90_RCONST" section 6. Manually write the "END SUBROUTINE UPDATE_RCONST" line .gitignore - Now ignore executable files in MCM example folders CHANGELOG.md - Updated accordingly Signed-off-by: Bob Yantosca --- .gitignore | 3 ++- CHANGELOG.md | 7 ++++- src/gdata.h | 10 +++---- src/gen.c | 75 ++++++++++++++++++++++++--------------------------- src/scanner.c | 1 + 5 files changed, 49 insertions(+), 47 deletions(-) diff --git a/.gitignore b/.gitignore index 00cec10..da5f2ba 100644 --- a/.gitignore +++ b/.gitignore @@ -67,5 +67,6 @@ docs/build/* # Other files/dirs to exclude *.pdf -/examples/mcm/__pycache__ +/examples/mcm*/__pycache__ +/examples/mcm*/*.exe diff --git a/CHANGELOG.md b/CHANGELOG.md index 45177ec..d372552 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,14 +9,19 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). ## [Unreleased] - TBD +### Added +- Added new inline key `F90_RCONST_USE` in `src/gdata.h` and `src/scanner.c` + ### Changed - Updated `Update_RCONST` to use `Y` instead of `C` to account for updated variable species concentrations - Updated C-I tests to print the compiler versions that are used +- Updated routine `GenerateUpdateRconst` to manually write the `SUBROUTINE` and `END SUBROUTINE` lines (F90 only) +- Updated routine `GenerateUpdateRconst` to inline code from `#INLINE F90_RCONST_USE` before any other F90 variable declarations or statements +- Updated `.gitignore` to ignore executables in the MCM example folders ### Fixed - Added `char* rootFileName` to functions and function prototypes for `Use_C`, `Use_F`, `Use_F90`, `Use_MATLAB`, and `Generate` - Updated `docs/requirements.txt` to use `jinja2==3.1.4` (fixes a security issue) -- Updated `gen.c` to write the `INLINED RCONST` section at the top of routine `UPDATE_RCONST`. This will prevent compilation errors when `INLINED_RCONST` contains F90 `USE` statements (as these must precede other F90 statements). ## [3.1.1] - 2024-04-30 ### Changed diff --git a/src/gdata.h b/src/gdata.h index 333412b..435337b 100644 --- a/src/gdata.h +++ b/src/gdata.h @@ -92,11 +92,11 @@ enum krtypes { NUMBER, EXPRESION, PHOTO }; enum table_modes { F_TEXT, FC_TEXT, C_TEXT, S_TEXT }; enum lang { NO_LANG, C_LANG, F77_LANG, F90_LANG, MATLAB_LANG }; enum inl_code { - F77_GLOBAL, F77_INIT, F77_DATA, F77_UTIL, F77_RATES, - F77_RCONST, F90_GLOBAL, F90_INIT, F90_DATA, F90_UTIL, - F90_RATES, F90_RCONST, C_GLOBAL, C_INIT, C_DATA, - C_UTIL, C_RATES, C_RCONST, MATLAB_GLOBAL, MATLAB_INIT, - MATLAB_DATA, MATLAB_UTIL, MATLAB_RATES, MATLAB_RCONST, + F77_GLOBAL, F77_INIT, F77_DATA, F77_UTIL, F77_RATES, + F77_RCONST, F90_GLOBAL, F90_INIT, F90_DATA, F90_UTIL, + F90_RATES, F90_RCONST, F90_RCONST_USE, C_GLOBAL, C_INIT, + C_DATA, C_UTIL, C_RATES, C_RCONST, MATLAB_GLOBAL, + MATLAB_INIT, MATLAB_DATA, MATLAB_UTIL, MATLAB_RATES, MATLAB_RCONST, INLINE_OPT }; diff --git a/src/gen.c b/src/gen.c index b19c819..eeeb5c3 100644 --- a/src/gen.c +++ b/src/gen.c @@ -2175,22 +2175,45 @@ int YIN, Y; // // SPECIAL HANDLING FOR F90: // ------------------------- - // Write the INLINED RCONST section immediately after the start of the - // subroutine UPDATE_RCONST. This will avoid compile-time errors if a - // "USE" statement is included via INLINED RCONST. Recall that F90 - // "USE" statements must precede variable declarations or any other - // executable statements. - // - // The FunctionBegin() routine writes the variable declaration - // statements immediately after the subroutine declaration line. - // Therefore, we will not be able to use FunctionBegin() to declare - // the UPDATE_RCONST subroutine. Instead, we will manually write - // "SUBROUTINE UPDATE_RCONST ( YIN )" here. + // Manually write the "SUBROUTINE RCONST ( YIN )" line instead of + // using the FunctionBegin( UPDATE_RCONST ). This will allow us + // to add any inlined F90 use statements immediately afterwards. + // Recall that F90 "USE" statements must precede variable + // declarations or any other executable statements, or else a + // compilation error will be generated. // // -- Bob Yantosca (19 Dec 2024) // bprintf("SUBROUTINE Update_RCONST ( YIN )"); + NewLines(2); + + // Inline USE statements right after the subroutine declaration + WriteComment("Begin INLINED RCONST - F90 USE STATEMENTS"); + NewLines(1); + bprintf( InlineCode[ F90_RCONST_USE ].code ); + FlushBuf(); NewLines(1); + WriteComment("End INLINED RCONST - F90 USE STATEMENTS"); + NewLines(1); + + // Declare optional YIN argument + YIN = DefvElmO( "YIN", real, -NVAR, "Optional input concentrations of variable species" ); + Declare(YIN); + NewLines(1); + + // Declare local Y variable + Y = DefvElm( "Y", real, -NSPEC, "Concentrations of species (local)" ); + Declare(Y); + NewLines(1); + + // Copy values of YIN to Y if YIN is present + WriteComment("Ensure local Y array is filled with variable and constant concentrations"); + bprintf(" Y(1:NSPEC) = C(1:NSPEC)\n"); + NewLines(1); + WriteComment("Update local Y array if variable concentrations are provided"); + bprintf(" if (present(YIN)) Y(1:NVAR) = YIN(1:NVAR)\n"); + NewLines(2); + } else { // // For other languages, declare function w/o any arguments @@ -2207,6 +2230,7 @@ int YIN, Y; NewLines(1); + // Inline code from {C,F77,F90_MATLAB}_RCONST inline keys NewLines(1); WriteComment("Begin INLINED RCONST"); NewLines(1); @@ -2226,35 +2250,6 @@ int YIN, Y; NewLines(1); WriteComment("End INLINED RCONST"); NewLines(1); - // - // SPECIAL HANDLING FOR F90: - // ------------------------- - // Write F90 variable declarations after the INLINED RCONST section. - // This will avoid compile-time errors if a USE statement is placed - // into the INLINED RCONST section, as described above. - // - // -- Bob Yantosca (19 Dec 2024) - // - if (useLang == F90_LANG) { - - // Declare optional YIN argument - YIN = DefvElmO( "YIN", real, -NVAR, "Optional input concentrations of variable species" ); - Declare(YIN); - NewLines(1); - - // Declare local Y variable - Y = DefvElm( "Y", real, -NSPEC, "Concentrations of species (local)" ); - Declare(Y); - NewLines(1); - - // Copy values of YIN to Y if YIN is present - WriteComment("Ensure local Y array is filled with variable and constant concentrations"); - bprintf(" Y(1:NSPEC) = C(1:NSPEC)\n"); - NewLines(1); - WriteComment("Update local Y array if variable concentrations are provided"); - bprintf(" if (present(YIN)) Y(1:NVAR) = YIN(1:NVAR)\n"); - NewLines(2); - } for( i = 0; i < EqnNr; i++) { /* mz_rs_20220701+ */ diff --git a/src/scanner.c b/src/scanner.c index abeaf52..06ca034 100644 --- a/src/scanner.c +++ b/src/scanner.c @@ -79,6 +79,7 @@ INLINE_KEY InlineKeys[] = { { F77_GLOBAL, APPEND, "F77_GLOBAL" }, { F90_UTIL, APPEND, "F90_UTIL" }, { F90_RATES, APPEND, "F90_RATES" }, { F90_RCONST, APPEND, "F90_RCONST" }, + { F90_RCONST_USE, APPEND, "F90_RCONST_USE"}, { C_GLOBAL, APPEND, "C_GLOBAL" }, { C_INIT, APPEND, "C_INIT" }, { C_DATA, APPEND, "C_DATA" }, From 2f978250e23f346a710dfafb659f1105a89933c0 Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Fri, 20 Dec 2024 15:26:42 -0500 Subject: [PATCH 13/28] Add minor fixes to PR #120 suggested by @RolfSander .gitignore - Ignore all *.exe files src/gen.c - Updated commments where inlined code is added to be more descriptive (referencing F90_RCONST_USE and F90_RCONST) docs/source/using_kpp/04_input_for_kpp.rst - Added documentation for F90_RCONST_USE CHANGELOG.md - Updated accordingly Signed-off-by: Bob Yantosca --- .gitignore | 2 +- CHANGELOG.md | 5 +- docs/source/using_kpp/04_input_for_kpp.rst | 81 ++++++++++++---------- src/gen.c | 10 +-- 4 files changed, 56 insertions(+), 42 deletions(-) diff --git a/.gitignore b/.gitignore index da5f2ba..6a87b70 100644 --- a/.gitignore +++ b/.gitignore @@ -15,6 +15,7 @@ # from compiler: *.mod *.o +*.exe # results from ci-tests: /ci-tests/**/*.m @@ -68,5 +69,4 @@ docs/build/* # Other files/dirs to exclude *.pdf /examples/mcm*/__pycache__ -/examples/mcm*/*.exe diff --git a/CHANGELOG.md b/CHANGELOG.md index d372552..6244891 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,13 +11,16 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] - TBD ### Added - Added new inline key `F90_RCONST_USE` in `src/gdata.h` and `src/scanner.c` +- Added documentation about `F90_RCONST_USE` for ReadTheDocs ### Changed - Updated `Update_RCONST` to use `Y` instead of `C` to account for updated variable species concentrations - Updated C-I tests to print the compiler versions that are used - Updated routine `GenerateUpdateRconst` to manually write the `SUBROUTINE` and `END SUBROUTINE` lines (F90 only) - Updated routine `GenerateUpdateRconst` to inline code from `#INLINE F90_RCONST_USE` before any other F90 variable declarations or statements -- Updated `.gitignore` to ignore executables in the MCM example folders +- Updated `.gitignore` to ignore all executable files +- Changed `Begin INLINED RCONST - F90 USE STATEMENTS` to `Begin inlined code from F90_RCONST_USE` in `src/gen.c` +- Changed `Begin INLINED RCONST` to `Begin inlined code from F90_RCONST` in `src/gen.c` ### Fixed - Added `char* rootFileName` to functions and function prototypes for `Use_C`, `Use_F`, `Use_F90`, `Use_MATLAB`, and `Generate` diff --git a/docs/source/using_kpp/04_input_for_kpp.rst b/docs/source/using_kpp/04_input_for_kpp.rst index 92502a4..01a4559 100644 --- a/docs/source/using_kpp/04_input_for_kpp.rst +++ b/docs/source/using_kpp/04_input_for_kpp.rst @@ -790,26 +790,28 @@ by :code:`C`, or :code:`matlab`, respectively. .. table:: KPP inlined types :align: center - +-----------------+-------------------+---------------------+---------------------+ - | Inline_type | File | Placement | Usage | - +=================+===================+=====================+=====================+ - | **F90_DATA** | :ref:`Monitor` | specification | (obsolete) | - | | | section | | - +-----------------+-------------------+---------------------+---------------------+ - | **F90_GLOBAL** | :ref:`Global` | specification | global variables | - | | | section | | - +-----------------+-------------------+---------------------+---------------------+ - | **F90_INIT** | :ref:`Initialize` | subroutine | integration | - | | | | parameters | - +-----------------+-------------------+---------------------+---------------------+ - | **F90_RATES** | :ref:`Rates` | executable section | rate law functions | - +-----------------+-------------------+---------------------+---------------------+ - | **F90_RCONST** | :ref:`Rates` | subroutine | statements and | - | | | | definitions of rate | - | | | | coefficients | - +-----------------+-------------------+---------------------+---------------------+ - | **F90_UTIL** | :ref:`Util` | executable section | utility functions | - +-----------------+-------------------+---------------------+---------------------+ + +--------------------+-------------------+---------------------+---------------------+ + | Inline_type | File | Placement | Usage | + +====================+===================+=====================+=====================+ + | **F90_DATA** | :ref:`Monitor` | specification | (obsolete) | + | | | section | | + +--------------------+-------------------+---------------------+---------------------+ + | **F90_GLOBAL** | :ref:`Global` | specification | global variables | + | | | section | | + +--------------------+-------------------+---------------------+---------------------+ + | **F90_INIT** | :ref:`Initialize` | subroutine | integration | + | | | | parameters | + +--------------------+-------------------+---------------------+---------------------+ + | **F90_RATES** | :ref:`Rates` | executable section | rate law functions | + +--------------------+-------------------+---------------------+---------------------+ + | **F90_RCONST** | :ref:`Rates` | subroutine | rate coefficient | + | | | | definitions | + +--------------------+-------------------+---------------------+---------------------+ + | **F90_RCONST_USE** | :ref:`Rates` | subroutine | rate coefficient | + | | | | definitions | + +--------------------+-------------------+---------------------+---------------------+ + | **F90_UTIL** | :ref:`Util` | executable section | utility functions | + +--------------------+-------------------+---------------------+---------------------+ .. _f90-data: @@ -909,29 +911,38 @@ F90_RCONST ---------- This inline type can be used to define time-dependent values of rate -coefficients. You may inline :code:`USE` statements that reference -modules where rate coefficients are computed, e.g.: +coefficients. You may inline variables directly, e.g.: .. code-block:: fortran #INLINE F90_RCONST - USE MyRateFunctionModule + k_DMS_OH = 1.E-9*EXP(5820./temp)*C(ind_O2)/ & + (1.E30+5.*EXP(6280./temp)*C(ind_O2)) #ENDINLINE -or define variables directly, e.g.: +The inlined code will be placed directly into the subroutines +:code:`UPDATE_RCONST` and :code:`UPDATE_PHOTO` in the :ref:`Rates` file. + +.. _f90-rconst-use: + +F90_RCONST_USE +-------------- + +Similar to :ref:`f90-rconst`, but allows you to inline Fortran-90 +:code:`USE` statements referencing modules where rate coefficients are +computed, such as: .. code-block:: fortran #INLINE F90_RCONST - k_DMS_OH = 1.E-9*EXP(5820./temp)*C(ind_O2)/ & - (1.E30+5.*EXP(6280./temp)*C(ind_O2)) + USE MyRateFunctionModule #ENDINLINE - -Note that the :code:`USE` statements must precede any variable -definitions. - + The inlined code will be placed directly into the subroutines -:code:`UPDATE_RCONST` and :code:`UPDATE_PHOTO` in the :ref:`Rates` file. +:code:`UPDATE_RCONST` in the :ref:`Rates` file. :code:`USE` +statements will be placed before Fortran variable definitions and +executable statements, as is required by the Fortran-90 language +standard. .. _f90-util: @@ -1078,7 +1089,7 @@ List of symbols replaced by the substitution preprocessor +--------------------------+-------------------------------+----------------------------+ .. _icntrl-rcntrl: - + ================================================================= Controlling the Integrator with :code:`ICNTRL` and :code:`RCNTRL` ================================================================= @@ -1226,13 +1237,13 @@ ICNTRL g.: .. code-block:: console - + HSO3m + HSO5m + Hp = 2 HSO4m + Hp : k_aqueous( C(ind_Hp) ); - + This ensures that the concentration :code:`C(ind_Hp)` at the specific integration time is used when the reaction rate coefficient is updated within the integrator. - + .. option:: ICNTRL(16) Treatment of negative concentrations: diff --git a/src/gen.c b/src/gen.c index eeeb5c3..98630c3 100644 --- a/src/gen.c +++ b/src/gen.c @@ -2188,12 +2188,12 @@ int YIN, Y; NewLines(2); // Inline USE statements right after the subroutine declaration - WriteComment("Begin INLINED RCONST - F90 USE STATEMENTS"); + WriteComment("Begin inlined code from F90_RCONST_USE"); NewLines(1); bprintf( InlineCode[ F90_RCONST_USE ].code ); FlushBuf(); NewLines(1); - WriteComment("End INLINED RCONST - F90 USE STATEMENTS"); + WriteComment("End inlined code from F90_RCONST_USE"); NewLines(1); // Declare optional YIN argument @@ -2212,7 +2212,7 @@ int YIN, Y; NewLines(1); WriteComment("Update local Y array if variable concentrations are provided"); bprintf(" if (present(YIN)) Y(1:NVAR) = YIN(1:NVAR)\n"); - NewLines(2); + NewLines(1); } else { // @@ -2232,7 +2232,7 @@ int YIN, Y; // Inline code from {C,F77,F90_MATLAB}_RCONST inline keys NewLines(1); - WriteComment("Begin INLINED RCONST"); + WriteComment("Begin inlined code from F90_RCONST"); NewLines(1); switch( useLang ) { @@ -2248,7 +2248,7 @@ int YIN, Y; FlushBuf(); NewLines(1); - WriteComment("End INLINED RCONST"); + WriteComment("End inlined code from F90_RCONST"); NewLines(1); for( i = 0; i < EqnNr; i++) { From ca50468b57a96f94f1e10a10ad71f7c85f1f0bbe Mon Sep 17 00:00:00 2001 From: Rolf Sander Date: Fri, 20 Dec 2024 21:50:45 +0100 Subject: [PATCH 14/28] consistent comments for inlined Code --- src/gen.c | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/src/gen.c b/src/gen.c index 98630c3..7b819eb 100644 --- a/src/gen.c +++ b/src/gen.c @@ -498,7 +498,7 @@ int dim; InitDeclare( FAM_NAMES, FamilyNr, (void*)sfam ); NewLines(1); - WriteComment("INLINED global variables"); + WriteComment("Begin inlined code from F90_DATA"); switch( useLang ) { case C_LANG: bprintf( InlineCode[ C_DATA ].code ); @@ -513,7 +513,7 @@ int dim; FlushBuf(); NewLines(1); - WriteComment("End INLINED global variables"); + WriteComment("End inlined code from F90_DATA"); NewLines(1); F77_Inline( "%6sEND\n\n", " " ); @@ -2122,7 +2122,7 @@ void GenerateRateLaws() NewLines(1); NewLines(1); - WriteComment("Begin INLINED Rate Law Functions"); + WriteComment("Begin inlined code from F90_RATES"); NewLines(1); switch( useLang ) { @@ -2138,7 +2138,7 @@ void GenerateRateLaws() FlushBuf(); NewLines(1); - WriteComment("End INLINED Rate Law Functions"); + WriteComment("End inlined code from F90_RATES"); NewLines(1); @@ -2308,7 +2308,7 @@ int UPDATE_PHOTO; MATLAB_Inline("global SUN TEMP RCONST"); NewLines(1); - WriteComment("Begin INLINED RCONST"); + WriteComment("Begin inlined code from F90_RCONST"); NewLines(1); switch( useLang ) { @@ -2324,7 +2324,7 @@ int UPDATE_PHOTO; FlushBuf(); NewLines(1); - WriteComment("End INLINED RCONST"); + WriteComment("End inlined code from F90_RCONST"); NewLines(1); for( i = 0; i < EqnNr; i++) { @@ -2389,7 +2389,7 @@ int UTIL; UseFile( utilFile ); NewLines(1); - WriteComment("User INLINED Utility Functions"); + WriteComment("Begin inlined code from F90_UTIL"); switch( useLang ) { case C_LANG: bprintf( InlineCode[ C_UTIL ].code ); @@ -2404,10 +2404,10 @@ int UTIL; FlushBuf(); NewLines(1); - WriteComment("End INLINED Utility Functions"); + WriteComment("End inlined code from F90_UTIL"); NewLines(1); - WriteComment("Utility Functions from KPP_HOME/util/util"); + WriteComment("Begin Utility Functions from KPP_HOME/util/util"); UTIL = DefFnc( "UTIL", 0, "Utility functions"); CommentFunctionBegin( UTIL); @@ -2684,7 +2684,7 @@ void GenerateGlobalHeader() } NewLines(1); - WriteComment("INLINED global variable declarations"); + WriteComment("Begin inlined code from F90_GLOBAL"); switch( useLang ) { case C_LANG: bprintf( InlineCode[ C_GLOBAL ].code ); @@ -2699,7 +2699,7 @@ void GenerateGlobalHeader() FlushBuf(); NewLines(1); - WriteComment("INLINED global variable declarations"); + WriteComment("End inlined code from F90_GLOBAL"); NewLines(1); } @@ -3014,15 +3014,15 @@ int INITVAL; } } - WriteComment("constant rate coefficients"); + WriteComment("Begin constant rate coefficients"); for( i = 0; i < EqnNr; i++) { if ( kr[i].type == NUMBER ) Assign( Elm( RCONST, i ), Const( kr[i].val.f ) ); } - WriteComment("END constant rate coefficients"); + WriteComment("End constant rate coefficients"); NewLines(1); - WriteComment("INLINED initializations"); + WriteComment("Begin inlined code from F90_INIT"); switch( useLang ) { case C_LANG: bprintf( InlineCode[ C_INIT ].code ); @@ -3037,7 +3037,7 @@ int INITVAL; FlushBuf(); NewLines(1); - WriteComment("End INLINED initializations"); + WriteComment("End inlined code from F90_INIT"); NewLines(1); MATLAB_Inline(" VAR = VAR(:);\n FIX = FIX(:);\n" ); From f191024130fee5e06aead00d875dfb24ada6999d Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Thu, 2 Jan 2025 10:16:50 -0500 Subject: [PATCH 15/28] Edited CHANGELOG.md for PR #122 CHANGELOG.md - Changed note about F90_RCONST_USE to be more general, as all inlined code comments have been updated to be more consistent in PR #122. Signed-off-by: Bob Yantosca --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 6244891..2064bfe 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -20,7 +20,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Updated routine `GenerateUpdateRconst` to inline code from `#INLINE F90_RCONST_USE` before any other F90 variable declarations or statements - Updated `.gitignore` to ignore all executable files - Changed `Begin INLINED RCONST - F90 USE STATEMENTS` to `Begin inlined code from F90_RCONST_USE` in `src/gen.c` -- Changed `Begin INLINED RCONST` to `Begin inlined code from F90_RCONST` in `src/gen.c` +- Changed inlined code comments to be more precise (e.g. `Begin inlined code from F90_RCONST`) in `src/gen.c` ### Fixed - Added `char* rootFileName` to functions and function prototypes for `Use_C`, `Use_F`, `Use_F90`, `Use_MATLAB`, and `Generate` From d12bb5c7e8b1ae58a74650c477061b7efcb6212c Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Tue, 14 Jan 2025 14:06:59 -0500 Subject: [PATCH 16/28] RTD update: Instruct users to set LD_LIBRARY_PATH for flex docs/source/getting_started/01_installation.html - Add a line to add the path to the flex library to LD_LIBRARY_PATH in the Flex installation instructions. CHANGELOG.md - Updated accordingly Signed-off-by: Bob Yantosca --- CHANGELOG.md | 1 + docs/source/getting_started/01_installation.rst | 3 ++- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2064bfe..1463ac0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -21,6 +21,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Updated `.gitignore` to ignore all executable files - Changed `Begin INLINED RCONST - F90 USE STATEMENTS` to `Begin inlined code from F90_RCONST_USE` in `src/gen.c` - Changed inlined code comments to be more precise (e.g. `Begin inlined code from F90_RCONST`) in `src/gen.c` +- Updated Flex library installation example on ReadTheDocs ### Fixed - Added `char* rootFileName` to functions and function prototypes for `Use_C`, `Use_F`, `Use_F90`, `Use_MATLAB`, and `Generate` diff --git a/docs/source/getting_started/01_installation.rst b/docs/source/getting_started/01_installation.rst index e7d59a8..aa3ebf9 100644 --- a/docs/source/getting_started/01_installation.rst +++ b/docs/source/getting_started/01_installation.rst @@ -178,8 +178,9 @@ the :envvar:`KPP_FLEX_LIB_DIR` environment variable in your .. code-block:: bash export KPP_FLEX_LIB_DIR=/usr/lib + export LD_LIBRARY_PATH="${LD_LIBRARY_PATH}:${KPP_FLEX_LIB_DIR}:" -And then apply the changes with: +Then apply the changes with: .. code-block:: console From 9887d0222b5eebb16d97d060d7b338eb5796a44e Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Mon, 27 Jan 2025 17:32:54 -0500 Subject: [PATCH 17/28] Bug fix: Inline F90_RCONST_USE code also to Update_Photo src/gen.c - In routine GenerateUpdatePhoto, we also now inline USE statements from inline block F90_RCONST_USE into routine Update_Photo. This fixes a bug that was found with the MCM mechanism. examples/mcm/mcm_isoprene.eqn - Moved "USE constants_mcm" from the F90_RCONST to the F90_RCONST_USE inline block. This will place it properly at the top of the routines. CHANGELOG.md - Updated accordingly Signed-off-by: Bob Yantosca --- CHANGELOG.md | 2 ++ examples/mcm/mcm_isoprene.eqn | 6 +++++- src/gen.c | 15 +++++++++++++++ 3 files changed, 22 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1463ac0..eaf5236 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added - Added new inline key `F90_RCONST_USE` in `src/gdata.h` and `src/scanner.c` - Added documentation about `F90_RCONST_USE` for ReadTheDocs +- Added `F90_RCONST_USE` inlined code to `Update_RConst`and `Update_Photo` routines ### Changed - Updated `Update_RCONST` to use `Y` instead of `C` to account for updated variable species concentrations @@ -26,6 +27,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Fixed - Added `char* rootFileName` to functions and function prototypes for `Use_C`, `Use_F`, `Use_F90`, `Use_MATLAB`, and `Generate` - Updated `docs/requirements.txt` to use `jinja2==3.1.4` (fixes a security issue) +- Moved `USE constants_mcm` from `F90_RCONST` to `F90_RCONST_USE` in `examples/mcm/mcm_isoprene.eqn` ## [3.1.1] - 2024-04-30 ### Changed diff --git a/examples/mcm/mcm_isoprene.eqn b/examples/mcm/mcm_isoprene.eqn index 022bfcc..b96b9d3 100644 --- a/examples/mcm/mcm_isoprene.eqn +++ b/examples/mcm/mcm_isoprene.eqn @@ -662,8 +662,12 @@ C536OOH = IGNORE ; C537O = IGNORE ; C537OOH = IGNORE ; -#INLINE F90_RCONST +#INLINE F90_RCONST_USE + ! Use statment to be inlined into Update_RCONST USE constants_mcm +#ENDINLINE + +#INLINE F90_RCONST ! Peroxy radicals ! WARNING: The following species do not have SMILES strings in the database. ! If any of these are peroxy radicals the RO2 sum will be wrong! diff --git a/src/gen.c b/src/gen.c index 7b819eb..13535b1 100644 --- a/src/gen.c +++ b/src/gen.c @@ -2300,6 +2300,21 @@ int UPDATE_PHOTO; UPDATE_PHOTO = DefFnc( "Update_PHOTO", 0, "function to update photolytical rate constants"); FunctionBegin( UPDATE_PHOTO ); + + if (useLang==F90_LANG) { + // + // SPECIAL HANDLING FOR F90: + // ------------------------- + // Inline USE statements right after the subroutine declaration + WriteComment("Begin inlined code from F90_RCONST_USE"); + NewLines(1); + bprintf( InlineCode[ F90_RCONST_USE ].code ); + FlushBuf(); + NewLines(1); + WriteComment("End inlined code from F90_RCONST_USE"); + NewLines(1); + } + F77_Inline(" INCLUDE '%s_Global.h'", rootFileName); /* mz_rs_20220212+ */ /* global is already used in the rates module, don't use it twice */ From 65615b921f45ba913aa75952657a875779d896fb Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Mon, 27 Jan 2025 17:41:24 -0500 Subject: [PATCH 18/28] Updated documentation about F90_RCONST_USE docs/source/using_kpp/04_input_for_kpp.rst - Now state that code specified in an #INLINE F90_RCONST_USE section will be placed into both Update_RCONST and Update_PHOTO. Signed-off-by: Bob Yantosca --- docs/source/using_kpp/04_input_for_kpp.rst | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/source/using_kpp/04_input_for_kpp.rst b/docs/source/using_kpp/04_input_for_kpp.rst index 01a4559..48713a2 100644 --- a/docs/source/using_kpp/04_input_for_kpp.rst +++ b/docs/source/using_kpp/04_input_for_kpp.rst @@ -939,10 +939,10 @@ computed, such as: #ENDINLINE The inlined code will be placed directly into the subroutines -:code:`UPDATE_RCONST` in the :ref:`Rates` file. :code:`USE` -statements will be placed before Fortran variable definitions and -executable statements, as is required by the Fortran-90 language -standard. +:code:`UPDATE_RCONST` and :code:`UPDATE_PHOTO` in the :ref:`Rates` +file. :code:`USE` statements will be placed before Fortran variable +definitions and executable statements, as is required by the +Fortran-90 language standard. .. _f90-util: From 0f15b2d17e90d3b0cd3974388e8ad43ed295ef56 Mon Sep 17 00:00:00 2001 From: Rolf Sander Date: Tue, 28 Jan 2025 08:35:28 +0100 Subject: [PATCH 19/28] bugfix in doc for F90_RCONST_USE --- docs/source/using_kpp/04_input_for_kpp.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/source/using_kpp/04_input_for_kpp.rst b/docs/source/using_kpp/04_input_for_kpp.rst index 48713a2..19d8ff7 100644 --- a/docs/source/using_kpp/04_input_for_kpp.rst +++ b/docs/source/using_kpp/04_input_for_kpp.rst @@ -934,7 +934,7 @@ computed, such as: .. code-block:: fortran - #INLINE F90_RCONST + #INLINE F90_RCONST_USE USE MyRateFunctionModule #ENDINLINE From bd2dfc9c4cd623c665bb66f7d921a2ddffb20691 Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Tue, 28 Jan 2025 13:58:24 -0500 Subject: [PATCH 20/28] Update documentation in advance of the 3.2.0 release docs/requirements.txt - Updated RTD dependency jinja2 to version 3.1.5 (security fix) docs/source/getting_started/00_revision_history.rst - Added what's new in KPP 3.2.0 section - Fixed tag position for KPP 3.1.0 anchor docs/source/reference/known-bugs.rst - Added links to both open and closed bug reports on GitHub - Added a blurb that LSODE is not thread-safe under known issues docs/source/tech_info/07_numerical_methods.rst - Added warning that LSODE is not thread-safe CHANGELOG.md - Updated accordingly Signed-off-by: Bob Yantosca --- CHANGELOG.md | 2 ++ docs/requirements.txt | 2 +- .../getting_started/00_revision_history.rst | 16 ++++++++- docs/source/reference/known-bugs.rst | 34 ++++++++++++++++--- .../source/tech_info/07_numerical_methods.rst | 9 +++++ 5 files changed, 56 insertions(+), 7 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1463ac0..0e08a01 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added - Added new inline key `F90_RCONST_USE` in `src/gdata.h` and `src/scanner.c` - Added documentation about `F90_RCONST_USE` for ReadTheDocs +- Added warning that LSODE is not thread-safe to ReadTheDocs documentation ### Changed - Updated `Update_RCONST` to use `Y` instead of `C` to account for updated variable species concentrations @@ -22,6 +23,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Changed `Begin INLINED RCONST - F90 USE STATEMENTS` to `Begin inlined code from F90_RCONST_USE` in `src/gen.c` - Changed inlined code comments to be more precise (e.g. `Begin inlined code from F90_RCONST`) in `src/gen.c` - Updated Flex library installation example on ReadTheDocs +- Updated ReadTheDocs dependency`jinja2` to version 3.1.5 (fixes a security issue) ### Fixed - Added `char* rootFileName` to functions and function prototypes for `Use_C`, `Use_F`, `Use_F90`, `Use_MATLAB`, and `Generate` diff --git a/docs/requirements.txt b/docs/requirements.txt index 71e269d..c0752da 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -12,4 +12,4 @@ sphinxcontrib-bibtex==2.6.2 sphinx-autobuild==2021.3.14 recommonmark==0.7.1 docutils==0.20.1 -jinja2==3.1.4 +jinja2==3.1.5 diff --git a/docs/source/getting_started/00_revision_history.rst b/docs/source/getting_started/00_revision_history.rst index c248f32..f4c342e 100644 --- a/docs/source/getting_started/00_revision_history.rst +++ b/docs/source/getting_started/00_revision_history.rst @@ -8,13 +8,25 @@ Only the major new features are listed here. For a detailed description of the changes, read `CHANGELOG.md `_. +.. _kpp320: + +========= +KPP 3.2.0 +========= + +- Added new inline key :literal:`F90_RCONST_USE` so that F90 + :literal:`USE` statements can be inlined into the + :code:`Update_RCONST` and :literal:`UPDATE_PHOTO` routines +- Updated code in :code:`src/gen.c` to generate the + :code:`UPDATE_RCONST` routine with an optional argument :code:`Y` +- Updated C-I tests to print the compiler versions that are used + .. _kpp311: ========= KPP 3.1.1 ========= -.. _kpp310: - Use newer Python packages to build ReadTheDocs documentation (see :file:`docs/requirements.txt`) - Increased :code:`MAX_NO_OF_LINES` and :code:`MAX_EQN` in order to @@ -22,6 +34,8 @@ KPP 3.1.1 - Now only add the extra `Aout` argument to `Fun` and `Fun_Split` for target language :literal:`Fortran90`. This fixes a Matlab build error. +.. _kpp310: + ========= KPP 3.1.0 ========= diff --git a/docs/source/reference/known-bugs.rst b/docs/source/reference/known-bugs.rst index 5668d7e..81d9ff5 100644 --- a/docs/source/reference/known-bugs.rst +++ b/docs/source/reference/known-bugs.rst @@ -1,6 +1,30 @@ -########## -Known Bugs -########## +##################### +Known bugs and issues +##################### -Bugs are discussed at the `KPP repository Github issues page -`_. +Please see our `Issue tracker on GitHub +`_ for a list of recent +bugs and fixes. + +=================== +Current bug reports +=================== + +These `bug reports (on GitHub) +`_ +are currently unresolved. We hope to fix these in future releases. + +LSODE integrator is not thread-safe +----------------------------------- + +We have discovered that the current implementation of the LSODE +integrator is not thread-safe for `OpenMP parallelization +`_. When LSODE is called from within an +OpenMP parallel loop, the integration will fail because key internal +variables in LSODE will be overwritten by concurrent threads. + +============================ +Bugs that have been resolved +============================ + +These `bugs (reported on GitHub) `_ have been resolved. diff --git a/docs/source/tech_info/07_numerical_methods.rst b/docs/source/tech_info/07_numerical_methods.rst index 1ac6a1d..2468ec2 100644 --- a/docs/source/tech_info/07_numerical_methods.rst +++ b/docs/source/tech_info/07_numerical_methods.rst @@ -534,6 +534,15 @@ LSODE, the Livermore ODE solver differentiation formula (BDF) methods for stiff problems. LSODE has been translated to Fortran90 for the incorporation into the KPP library. +.. attention:: + + We have discovered that the current implementation of the LSODE + integrator is not thread-safe for `OpenMP parallelization + `_. When LSODE is called from within + an OpenMP parallel loop, the integration will fail because key + internal variables in LSODE will be overwritten by concurrent + threads. + VODE ---- From 5c70d207f46cf63c0e04d6aa93503b7b2070925b Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Thu, 6 Feb 2025 15:46:58 -0500 Subject: [PATCH 21/28] Fixed MacOS-specific handling for x86_64 or arm64 src/Makefile - Now print the value of SYSTEM_M variable in "make list" target - Trimmed trailing whitespace src/Makefile.defs - Now use the $(strip) function when testing the values of SYSTEM and SYSTEM_M. This will make sure that trailing whitespace is removed. This was affecting the string equality tests. - Add option "-arch x86_64" if compiling on Intel x86_64 chipsets - Add option "-arch arm64" if compiling on Apple Silicon chipsets - Trimmed trailing whitespace CHANGELOG.md - Updated accordingly Signed-off-by: Bob Yantosca --- CHANGELOG.md | 1 + src/Makefile | 11 ++++++----- src/Makefile.defs | 24 +++++++++++++++++------- 3 files changed, 24 insertions(+), 12 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index eaf5236..b738c3c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -28,6 +28,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Added `char* rootFileName` to functions and function prototypes for `Use_C`, `Use_F`, `Use_F90`, `Use_MATLAB`, and `Generate` - Updated `docs/requirements.txt` to use `jinja2==3.1.4` (fixes a security issue) - Moved `USE constants_mcm` from `F90_RCONST` to `F90_RCONST_USE` in `examples/mcm/mcm_isoprene.eqn` +- Fixed MacOS-specific handling for x86_64 or arm64 in `src/Makefile.defs` ## [3.1.1] - 2024-04-30 ### Changed diff --git a/src/Makefile b/src/Makefile index 72622c8..1f6cae8 100644 --- a/src/Makefile +++ b/src/Makefile @@ -5,7 +5,7 @@ # # Copyright (C) 1995-1996 Valeriu Damian and Adrian Sandu # Copyright (C) 1997-2005 Adrian Sandu -# with contributions from: Mirela Damian, Rolf Sander +# with contributions from: Mirela Damian, Rolf Sander # # KPP is free software; you can redistribute it and/or modify it under the # terms of the GNU General Public License as published by the Free Software @@ -56,8 +56,8 @@ all: $(OBJS) @mv kpp $(PROG) .PHONY: clean -clean: - @rm -f *~ *.o +clean: + @rm -f *~ *.o .PHONY: distclean distclean: clean @@ -71,6 +71,7 @@ maintainer-clean: distclean .PHONY: list list: @echo "SYSTEM = $(SYSTEM)" + @echo "SYSTEM_M = $(SYSTEM_M)" @echo "HOST = $(HOST)" @echo "CC = $(CC)" @echo "CC_FLAGS = $(CC_FLAGS)" @@ -83,7 +84,7 @@ list: $(CC) $(CC_FLAGS) -I$(INCLUDE_DIR) -c $*.c lex.yy.c: scan.l scan.h gdata.h - $(FLEX) -olex.yy.c scan.l + $(FLEX) -olex.yy.c scan.l y.tab.c: scan.y scan.h $(BISON) -d -o y.tab.c scan.y @@ -101,6 +102,6 @@ debug.o: gdata.h scan.h gen.o: code.h gdata.h scan.h kpp.o: gdata.h scan.h lex.yy.o: gdata.h scan.h y.tab.h -scanner.o: gdata.h scan.h +scanner.o: gdata.h scan.h scanutil.o: gdata.h scan.h y.tab.o: scan.h y.tab.h diff --git a/src/Makefile.defs b/src/Makefile.defs index 1de3468..2199256 100644 --- a/src/Makefile.defs +++ b/src/Makefile.defs @@ -73,18 +73,18 @@ FLEX := flex FLEX_LIB_DIR := /usr/lib BISON := bison INCLUDE_DIR := /usr/include -SYSTEM := $(shell uname) # returns e.g. "Linux" -SYSTEM_M := $(shell uname -m) # returns e.g. "x86_64" -HOST := $(shell hostname) # returns name of machine +SYSTEM := $(shell uname) # e.g. "Linux", "Darwin" +SYSTEM_M := $(shell uname -m) # e.g. "x86_64", "arm64" +HOST := $(shell hostname) # name of machine ############################################################ ### Keep looking for the flex library until found ### ### Otherwise use the path specified in KPP_FLEX_LIB_DIR ### ### ### ### NOTE for MacOS X: IF you have installed flex with ### -### HomeBrew, then flex will be installed to a path such ### +### HomeBrew, then flex will be installed to a path such ### ### as /usr/local/Cellar/flex/X.Y.Z/lib, where X.Y.Z is ### -### the version number. -- Bob Yantosca (01 Nov 2021) ### +### the version number. -- Bob Yantosca (01 Nov 2021) ### ############################################################ # Try /usr/lib64 @@ -128,8 +128,18 @@ endif ############################################################ # Settings for MacOS -ifeq ($(SYSTEM),Darwin) +ifeq ($(strip $(SYSTEM)),Darwin) CC_FLAGS += -DMACOS -O + ifeq ($(strip $(SYSTEM_M)),x86_64) + # Build for Intel x86_64 chipsets + CC_FLAGS += -arch x86_64 + else + ifeq ($(strip $(SYSTEM_M)),arm64) + # Build for Apple Silicon chipsets + CC_FLAGS += -arch arm64 + else + $(warning Unknown architecture: $(strip $(SYSTEM_M))) + endif + endif endif - ############################################################# From 99fd9b5a889f896635cba3ef120e9bf142e1049c Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Tue, 18 Feb 2025 16:13:17 -0500 Subject: [PATCH 22/28] Replace int/sdirk.f90 module with int/beuler.f90 int/sdirk.f90 - Original version has been deleted - int/beuler.f90 is now renamed to int/sdirk.f90 int/beuler.f90 int/beuler.def - Removed docs/source/tech_info/07_numerical_methods.rst - Added documentation instructing user to select #INTEGRATOR sdirk with ICNTRL(3) = 6 for Backwards Euler integration CHANGELOG.md - Updated accordingly Signed-off-by: Bob Yantosca --- CHANGELOG.md | 6 + .../source/tech_info/07_numerical_methods.rst | 11 +- int/beuler.def | 1 - int/beuler.f90 | 1384 ----------------- int/sdirk.f90 | 162 +- 5 files changed, 132 insertions(+), 1432 deletions(-) delete mode 100644 int/beuler.def delete mode 100644 int/beuler.f90 diff --git a/CHANGELOG.md b/CHANGELOG.md index eaf5236..57381bc 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -23,12 +23,18 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Changed `Begin INLINED RCONST - F90 USE STATEMENTS` to `Begin inlined code from F90_RCONST_USE` in `src/gen.c` - Changed inlined code comments to be more precise (e.g. `Begin inlined code from F90_RCONST`) in `src/gen.c` - Updated Flex library installation example on ReadTheDocs +- Renamed `int/beuler.f90` to the `int/sdirk.f90`, as this is a newer version of the SDIRK integrator +- Updated documentation for Backwards Euler to instruct user to select `#INTEGRATOR sdirk` with `ICNTRL(3) = 6` ### Fixed - Added `char* rootFileName` to functions and function prototypes for `Use_C`, `Use_F`, `Use_F90`, `Use_MATLAB`, and `Generate` - Updated `docs/requirements.txt` to use `jinja2==3.1.4` (fixes a security issue) - Moved `USE constants_mcm` from `F90_RCONST` to `F90_RCONST_USE` in `examples/mcm/mcm_isoprene.eqn` +### Removed +- Removed `int/beuler.f90` +- Removed `int/beuler.def` + ## [3.1.1] - 2024-04-30 ### Changed - Updated Python package versions for ReadTheDocs in `docs/requirements.txt` diff --git a/docs/source/tech_info/07_numerical_methods.rst b/docs/source/tech_info/07_numerical_methods.rst index 1ac6a1d..fd79705 100644 --- a/docs/source/tech_info/07_numerical_methods.rst +++ b/docs/source/tech_info/07_numerical_methods.rst @@ -546,9 +546,16 @@ the KPP library uses directly the KPP sparse linear algebra routines. BEULER ------ -**Integrator file:** :file:`int/beuler.f90` +**Integrator file:** :file:`int/sdirk.f90` -Backward Euler integration method. +Backward Euler integration method. To request this method, make sure +you select + +.. code-block:: console + + #INTEGRATOR sdirk + +in your definition file, and then set :code:`ICNTRL(3) = 6`. .. _other-methods: diff --git a/int/beuler.def b/int/beuler.def deleted file mode 100644 index 4d45f49..0000000 --- a/int/beuler.def +++ /dev/null @@ -1 +0,0 @@ -#INTFILE beuler diff --git a/int/beuler.f90 b/int/beuler.f90 deleted file mode 100644 index 0e79fdb..0000000 --- a/int/beuler.f90 +++ /dev/null @@ -1,1384 +0,0 @@ -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! -! SDIRK - Singly-Diagonally-Implicit Runge-Kutta methods ! -! * Sdirk 2a, 2b: L-stable, 2 stages, order 2 ! -! * Sdirk 3a: L-stable, 3 stages, order 2, adj-invariant ! -! * Sdirk 4a, 4b: L-stable, 5 stages, order 4 ! -! * Backward Euler: L-stable, 1 stage, order 1 ! -! Fixed step size = RCNTRL(3) = Hstart ! -! By default the code employs the KPP sparse linear algebra routines ! -! Compile with -DFULL_ALGEBRA to use full linear algebra (LAPACK) ! -! ! -! (C) Adrian Sandu, July 2005 ! -! Virginia Polytechnic Institute and State University ! -! Contact: sandu@cs.vt.edu ! -! Revised by Philipp Miehe and Adrian Sandu, May 2006 ! -! Backward Euler added by Adrian Sandu, December 2012 ! -! This implementation is part of KPP - the Kinetic PreProcessor ! -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! - -MODULE KPP_ROOT_Integrator - - USE KPP_ROOT_Precision - USE KPP_ROOT_Global - USE KPP_ROOT_Parameters - USE KPP_ROOT_JacobianSP, ONLY : LU_DIAG - USE KPP_ROOT_LinearAlgebra, ONLY : KppDecomp, KppSolve, Set2zero, & - WLAMCH, WCOPY, WAXPY, & - WSCAL, WADD - - IMPLICIT NONE - PUBLIC - SAVE - -!~~~> Flags to determine if we should call the UPDATE_* routines from within -!~~~> the integrator. If using KPP in an external model, you might want to -!~~~> disable these calls (via ICNTRL(15)) to avoid excess computations. - LOGICAL, PRIVATE :: Do_Update_RCONST - LOGICAL, PRIVATE :: Do_Update_PHOTO - LOGICAL, PRIVATE :: Do_Update_SUN - -!~~~> Statistics on the work performed by the SDIRK method - INTEGER, PARAMETER :: Nfun=1, Njac=2, Nstp=3, Nacc=4, & - Nrej=5, Ndec=6, Nsol=7, Nsng=8, & - Ntexit=1, Nhexit=2, Nhnew=3 - - ! mz_rs_20171120+ - ! description of the error numbers IERR - CHARACTER(LEN=50), PARAMETER, DIMENSION(-1:1) :: IERR_NAMES = (/ & - 'DUMMY ERROR MESSAGE ', & ! -1 - ' ', & ! 0 (not used) - 'Success ' /) ! 1 - ! mz_rs_20171120- - -CONTAINS - -SUBROUTINE INTEGRATE( TIN, TOUT, ICNTRL_U, RCNTRL_U, & - ISTATUS_U, RSTATUS_U, Ierr_U ) - - USE KPP_ROOT_Util, ONLY : Integrator_Update_Options - - IMPLICIT NONE - - KPP_REAL, INTENT(IN) :: TIN ! Start Time - KPP_REAL, INTENT(IN) :: TOUT ! End Time - ! Optional input parameters and statistics - INTEGER, INTENT(IN), OPTIONAL :: ICNTRL_U(20) - KPP_REAL, INTENT(IN), OPTIONAL :: RCNTRL_U(20) - INTEGER, INTENT(OUT), OPTIONAL :: ISTATUS_U(20) - KPP_REAL, INTENT(OUT), OPTIONAL :: RSTATUS_U(20) - INTEGER, INTENT(OUT), OPTIONAL :: Ierr_U - - INTEGER, SAVE :: Ntotal = 0 - KPP_REAL :: RCNTRL(20), RSTATUS(20), T1, T2 - INTEGER :: ICNTRL(20), ISTATUS(20), Ierr - - !~~~> Zero input and output arrays for safety's sake - ICNTRL = 0 - RCNTRL = 0.0_dp - ISTATUS = 0 - RSTATUS = 0.0_dp - - !~~~> fine-tune the integrator: - ICNTRL(2) = 0 ! 0 - vector tolerances, 1 - scalar tolerances - ICNTRL(6) = 0 ! starting values of Newton iterations: - ! interpolated (0), zero (1) - ICNTRL(15) = 5 ! Call Update_SUN and Update_RCONST from w/in the int. - - !~~~> if optional parameters are given, and if they are /= 0, - ! then use them to overwrite default settings - IF ( PRESENT( ICNTRL_U ) ) THEN - WHERE( ICNTRL_U /= 0 ) ICNTRL = ICNTRL_U - ENDIF - IF ( PRESENT( RCNTRL_U ) ) THEN - WHERE( RCNTRL_U > 0 ) RCNTRL = RCNTRL_U - ENDIF - - !~~~> Determine the settings of the Do_Update_* flags, which determine - !~~~> whether or not we need to call Update_* routines in the integrator - !~~~> (or not, if we are calling them from a higher-level) - ! ICNTRL(15) = -1 ! Do not call Update_* functions within the integrator - ! = 0 ! Status quo - ! = 1 ! Call Update_RCONST from within the integrator - ! = 2 ! Call Update_PHOTO from within the integrator - ! = 3 ! Call Update_RCONST and Update_PHOTO from w/in the int. - ! = 4 ! Call Update_SUN from within the integrator - ! = 5 ! Call Update_SUN and Update_RCONST from within the int. - ! = 6 ! Call Update_SUN and Update_PHOTO from within the int. - ! = 7 ! Call Update_SUN, Update_PHOTO, Update_RCONST w/in int. - CALL Integrator_Update_Options( ICNTRL(15), & - Do_Update_RCONST, & - Do_Update_PHOTO, & - Do_Update_Sun ) - - !~~~> In order to remove the prior EQUIVALENCE statements (which - !~~~> are not thread-safe), we now have declared VAR and FIX as - !~~~> threadprivate pointer variables that can point to C. - VAR => C(1:NVAR ) - FIX => C(NVAR+1:NSPEC) - - !~~~> Call the integrator - T1 = TIN; T2 = TOUT - CALL SDIRK( NVAR,T1,T2,VAR,RTOL,ATOL, & - RCNTRL,ICNTRL,RSTATUS,ISTATUS,Ierr ) - - !~~~> Free pointers - VAR => NULL() - FIX => NULL() - - !~~~> Debug option: print number of steps - ! Ntotal = Ntotal + ISTATUS(Nstp) - ! PRINT*,'NSTEPS=',ISTATUS(Nstp), '(',Ntotal,')',' O3=',VAR(ind_O3), & - ! ' NO2=',VAR(ind_NO2) - - IF (Ierr < 0) THEN - PRINT *,'SDIRK: Unsuccessful exit at T=',TIN,' (Ierr=',Ierr,')' - ENDIF - - !~~~> if optional parameters are given for output - !~~~> use them to store information in them - IF ( PRESENT( ISTATUS_U ) ) ISTATUS_U = ISTATUS - IF ( PRESENT( RSTATUS_U ) ) RSTATUS_U = RSTATUS - IF ( PRESENT( IERR_U ) ) IERR_U = IERR - - END SUBROUTINE INTEGRATE - - -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - SUBROUTINE SDIRK(N, Tinitial, Tfinal, Y, RelTol, AbsTol, & - RCNTRL, ICNTRL, RSTATUS, ISTATUS, Ierr) -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -! -! Solves the system y'=F(t,y) using a Singly-Diagonally-Implicit -! Runge-Kutta (SDIRK) method. -! -! This implementation is based on the book and the code Sdirk4: -! -! E. Hairer and G. Wanner -! "Solving ODEs II. Stiff and differential-algebraic problems". -! Springer series in computational mathematics, Springer-Verlag, 1996. -! This code is based on the SDIRK4 routine in the above book. -! -! Methods: -! * Sdirk 2a, 2b: L-stable, 2 stages, order 2 -! * Sdirk 3a: L-stable, 3 stages, order 2, adjoint-invariant -! * Sdirk 4a, 4b: L-stable, 5 stages, order 4 -! * Backward Euler: L-stable, 1 stage, order 1 -! Fixed step size = RCNTRL(3) = Hstart -! -! (C) Adrian Sandu, July 2005 -! Virginia Polytechnic Institute and State University -! Contact: sandu@cs.vt.edu -! Revised by Philipp Miehe and Adrian Sandu, May 2006 -! Backward Euler added by Adrian Sandu, Dec. 2012 -! This implementation is part of KPP - the Kinetic PreProcessor -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -! -!~~~> INPUT ARGUMENTS: -! -!- Y(NVAR) = vector of initial conditions (at T=Tinitial) -!- [Tinitial,Tfinal] = time range of integration -! (if Tinitial>Tfinal the integration is performed backwards in time) -!- RelTol, AbsTol = user precribed accuracy -!- SUBROUTINE ode_Fun( T, Y, Ydot ) = ODE function, -! returns Ydot = Y' = F(T,Y) -!- SUBROUTINE ode_Fun( T, Y, Ydot ) = Jacobian of the ODE function, -! returns Jcb = dF/dY -!- ICNTRL(1:20) = integer inputs parameters -!- RCNTRL(1:20) = real inputs parameters -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -! -!~~~> OUTPUT ARGUMENTS: -! -!- Y(NVAR) -> vector of final states (at T->Tfinal) -!- ISTATUS(1:20) -> integer output parameters -!- RSTATUS(1:20) -> real output parameters -!- Ierr -> job status upon return -! success (positive value) or -! failure (negative value) -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -! -!~~~> INPUT PARAMETERS: -! -! Note: For input parameters equal to zero the default values of the -! corresponding variables are used. -! -! Note: For input parameters equal to zero the default values of the -! corresponding variables are used. -!~~~> -! ICNTRL(1) = not used -! -! ICNTRL(2) = 0: AbsTol, RelTol are NVAR-dimensional vectors -! = 1: AbsTol, RelTol are scalars -! -! ICNTRL(3) = Method -! -! ICNTRL(4) -> maximum number of integration steps -! For ICNTRL(4)=0 the default value of 200000 is used -! -! ICNTRL(5) -> maximum number of Newton iterations -! For ICNTRL(5)=0 the default value of 8 is used -! -! ICNTRL(6) -> starting values of Newton iterations: -! ICNTRL(6)=0 : starting values are interpolated (the default) -! ICNTRL(6)=1 : starting values are zero -! -! ICNTRL(15) -> Toggles calling of Update_* functions w/in the integrator -! = -1 : Do not call Update_* functions within the integrator -! = 0 : Status quo -! = 1 : Call Update_RCONST from within the integrator -! = 2 : Call Update_PHOTO from within the integrator -! = 3 : Call Update_RCONST and Update_PHOTO from w/in the int. -! = 4 : Call Update_SUN from within the integrator -! = 5 : Call Update_SUN and Update_RCONST from within the int. -! = 6 : Call Update_SUN and Update_PHOTO from within the int. -! = 7 : Call Update_SUN, Update_PHOTO, Update_RCONST from w/in the int. -! -!~~~> Real parameters -! -! RCNTRL(1) -> Hmin, lower bound for the integration step size -! It is strongly recommended to keep Hmin = ZERO -! RCNTRL(2) -> Hmax, upper bound for the integration step size -! RCNTRL(3) -> Hstart, starting value for the integration step size -! -! RCNTRL(4) -> FacMin, lower bound on step decrease factor (default=0.2) -! RCNTRL(5) -> FacMax, upper bound on step increase factor (default=6) -! RCNTRL(6) -> FacRej, step decrease factor after multiple rejections -! (default=0.1) -! RCNTRL(7) -> FacSafe, by which the new step is slightly smaller -! than the predicted value (default=0.9) -! RCNTRL(8) -> ThetaMin. If Newton convergence rate smaller -! than ThetaMin the Jacobian is not recomputed; -! (default=0.001) -! RCNTRL(9) -> NewtonTol, stopping criterion for Newton's method -! (default=0.03) -! RCNTRL(10) -> Qmin -! RCNTRL(11) -> Qmax. If Qmin < Hnew/Hold < Qmax, then the -! step size is kept constant and the LU factorization -! reused (default Qmin=1, Qmax=1.2) -! -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -! -!~~~> OUTPUT PARAMETERS: -! -! Note: each call to Rosenbrock adds the current no. of fcn calls -! to previous value of ISTATUS(1), and similar for the other params. -! Set ISTATUS(1:10) = 0 before call to avoid this accumulation. -! -! ISTATUS(1) = No. of function calls -! ISTATUS(2) = No. of jacobian calls -! ISTATUS(3) = No. of steps -! ISTATUS(4) = No. of accepted steps -! ISTATUS(5) = No. of rejected steps (except at the beginning) -! ISTATUS(6) = No. of LU decompositions -! ISTATUS(7) = No. of forward/backward substitutions -! ISTATUS(8) = No. of singular matrix decompositions -! -! RSTATUS(1) -> Texit, the time corresponding to the -! computed Y upon return -! RSTATUS(2) -> Hexit,last accepted step before return -! RSTATUS(3) -> Hnew, last predicted step before return -! For multiple restarts, use Hnew as Hstart in the following run -! -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - IMPLICIT NONE - -! Arguments - INTEGER, INTENT(IN) :: N, ICNTRL(20) - KPP_REAL, INTENT(IN) :: Tinitial, Tfinal, & - RelTol(NVAR), AbsTol(NVAR), RCNTRL(20) - KPP_REAL, INTENT(INOUT) :: Y(NVAR) - INTEGER, INTENT(OUT) :: Ierr - INTEGER, INTENT(INOUT) :: ISTATUS(20) - KPP_REAL, INTENT(OUT) :: RSTATUS(20) - -!~~~> SDIRK method coefficients, up to 5 stages - INTEGER, PARAMETER :: Smax = 5 - INTEGER, PARAMETER :: S2A=1, S2B=2, S3A=3, S4A=4, S4B=5, BEL=6 - KPP_REAL :: rkGamma, rkA(Smax,Smax), rkB(Smax), rkC(Smax), & - rkD(Smax), rkE(Smax), rkBhat(Smax), rkELO, & - rkAlpha(Smax,Smax), rkTheta(Smax,Smax) - INTEGER :: sdMethod, rkS ! The number of stages -! Local variables - LOGICAL :: StartNewton - KPP_REAL :: Hmin, Hmax, Hstart, Roundoff, & - FacMin, Facmax, FacSafe, FacRej, & - ThetaMin, NewtonTol, Qmin, Qmax - INTEGER :: ITOL, NewtonMaxit, Max_no_steps, i - KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0 - - - ISTATUS(1:20) = 0 - RSTATUS(1:20) = ZERO - Ierr = 0 - -!~~~> For Scalar tolerances (ICNTRL(2).NE.0) the code uses AbsTol(1) and RelTol(1) -! For Vector tolerances (ICNTRL(2) == 0) the code uses AbsTol(1:NVAR) and RelTol(1:NVAR) - IF (ICNTRL(2) == 0) THEN - ITOL = 1 - ELSE - ITOL = 0 - END IF - -!~~~> ICNTRL(3) - method selection - SELECT CASE (ICNTRL(3)) - CASE (0,1) - CALL Sdirk2a - CASE (2) - CALL Sdirk2b - CASE (3) - CALL Sdirk3a - CASE (4) - CALL Sdirk4a - CASE (5) - CALL Sdirk4b - CASE (6) - CALL BEuler - CASE DEFAULT - CALL Sdirk2a - END SELECT - -!~~~> The maximum number of time steps admitted - IF (ICNTRL(4) == 0) THEN - Max_no_steps = 200000 - ELSEIF (ICNTRL(4) > 0) THEN - Max_no_steps=ICNTRL(4) - ELSE - PRINT * ,'User-selected ICNTRL(4)=',ICNTRL(4) - CALL SDIRK_ErrorMsg(-1,Tinitial,ZERO,Ierr) - END IF - -!~~~> The maximum number of Newton iterations admitted - IF(ICNTRL(5) == 0)THEN - NewtonMaxit=8 - ELSE - NewtonMaxit=ICNTRL(5) - IF(NewtonMaxit <= 0)THEN - PRINT * ,'User-selected ICNTRL(5)=',ICNTRL(5) - CALL SDIRK_ErrorMsg(-2,Tinitial,ZERO,Ierr) - END IF - END IF - -!~~~> StartNewton: Use extrapolation for starting values of Newton iterations - IF (ICNTRL(6) == 0) THEN - StartNewton = .TRUE. - ELSE - StartNewton = .FALSE. - END IF - -!~~~> Unit roundoff (1+Roundoff>1) - Roundoff = WLAMCH('E') - -!~~~> Lower bound on the step size: (positive value) - IF (RCNTRL(1) == ZERO) THEN - Hmin = ZERO - ELSEIF (RCNTRL(1) > ZERO) THEN - Hmin = RCNTRL(1) - ELSE - PRINT * , 'User-selected RCNTRL(1)=', RCNTRL(1) - CALL SDIRK_ErrorMsg(-3,Tinitial,ZERO,Ierr) - END IF - -!~~~> Upper bound on the step size: (positive value) - IF (RCNTRL(2) == ZERO) THEN - Hmax = ABS(Tfinal-Tinitial) - ELSEIF (RCNTRL(2) > ZERO) THEN - Hmax = MIN(ABS(RCNTRL(2)),ABS(Tfinal-Tinitial)) - ELSE - PRINT * , 'User-selected RCNTRL(2)=', RCNTRL(2) - CALL SDIRK_ErrorMsg(-3,Tinitial,ZERO,Ierr) - END IF - -!~~~> Starting step size: (positive value) - IF (RCNTRL(3) == ZERO) THEN - Hstart = MAX(Hmin,101.*Roundoff) ! factor 101 added mz_rs_20171120 - ELSEIF (RCNTRL(3) > ZERO) THEN - Hstart = MIN(ABS(RCNTRL(3)),ABS(Tfinal-Tinitial)) - ELSE - PRINT * , 'User-selected Hstart: RCNTRL(3)=', RCNTRL(3) - CALL SDIRK_ErrorMsg(-3,Tinitial,ZERO,Ierr) - END IF - -!~~~> Step size can be changed s.t. FacMin < Hnew/Hexit < FacMax - IF (RCNTRL(4) == ZERO) THEN - FacMin = 0.2_dp - ELSEIF (RCNTRL(4) > ZERO) THEN - FacMin = RCNTRL(4) - ELSE - PRINT * , 'User-selected FacMin: RCNTRL(4)=', RCNTRL(4) - CALL SDIRK_ErrorMsg(-4,Tinitial,ZERO,Ierr) - END IF - IF (RCNTRL(5) == ZERO) THEN - FacMax = 10.0_dp - ELSEIF (RCNTRL(5) > ZERO) THEN - FacMax = RCNTRL(5) - ELSE - PRINT * , 'User-selected FacMax: RCNTRL(5)=', RCNTRL(5) - CALL SDIRK_ErrorMsg(-4,Tinitial,ZERO,Ierr) - END IF -!~~~> FacRej: Factor to decrease step after 2 succesive rejections - IF (RCNTRL(6) == ZERO) THEN - FacRej = 0.1_dp - ELSEIF (RCNTRL(6) > ZERO) THEN - FacRej = RCNTRL(6) - ELSE - PRINT * , 'User-selected FacRej: RCNTRL(6)=', RCNTRL(6) - CALL SDIRK_ErrorMsg(-4,Tinitial,ZERO,Ierr) - END IF -!~~~> FacSafe: Safety Factor in the computation of new step size - IF (RCNTRL(7) == ZERO) THEN - FacSafe = 0.9_dp - ELSEIF (RCNTRL(7) > ZERO) THEN - FacSafe = RCNTRL(7) - ELSE - PRINT * , 'User-selected FacSafe: RCNTRL(7)=', RCNTRL(7) - CALL SDIRK_ErrorMsg(-4,Tinitial,ZERO,Ierr) - END IF - -!~~~> ThetaMin: decides whether the Jacobian should be recomputed - IF(RCNTRL(8) == 0.D0)THEN - ThetaMin = 1.0d-3 - ELSE - ThetaMin = RCNTRL(8) - END IF - -!~~~> Stopping criterion for Newton's method - IF(RCNTRL(9) == ZERO)THEN - NewtonTol = 3.0d-2 - ELSE - NewtonTol = RCNTRL(9) - END IF - -!~~~> Qmin, Qmax: IF Qmin < Hnew/Hold < Qmax, STEP SIZE = CONST. - IF(RCNTRL(10) == ZERO)THEN - Qmin=ONE - ELSE - Qmin=RCNTRL(10) - END IF - IF(RCNTRL(11) == ZERO)THEN - Qmax=1.2D0 - ELSE - Qmax=RCNTRL(11) - END IF - -!~~~> Check if tolerances are reasonable - IF (ITOL == 0) THEN - IF (AbsTol(1) <= ZERO .OR. RelTol(1) <= 10.D0*Roundoff) THEN - PRINT * , ' Scalar AbsTol = ',AbsTol(1) - PRINT * , ' Scalar RelTol = ',RelTol(1) - CALL SDIRK_ErrorMsg(-5,Tinitial,ZERO,Ierr) - END IF - ELSE - DO i=1,N - IF (AbsTol(i) <= 0.D0.OR.RelTol(i) <= 10.D0*Roundoff) THEN - PRINT * , ' AbsTol(',i,') = ',AbsTol(i) - PRINT * , ' RelTol(',i,') = ',RelTol(i) - CALL SDIRK_ErrorMsg(-5,Tinitial,ZERO,Ierr) - END IF - END DO - END IF - -!~~~> Special treatment of backward Euler - IF (sdMethod == BEL) THEN - StartNewton = .FALSE. - FacMin = 1.0_dp - FacMax = 1.0_dp - FacSafe= 1.0_dp - FacRej = 1.0_dp - Qmin = 1.0_dp - Qmax = 1.0_dp - Hmax = Hstart - IF ( Hstart <= 100_dp*Roundoff ) THEN - PRINT * , ' Backward Euler calling error:' - PRINT * , ' it requires Hstart > 100*Roundoff' - STOP - END IF - END IF - - - IF (Ierr < 0) RETURN - - CALL SDIRK_Integrator( N,Tinitial,Tfinal,Y,Ierr ) - - - !END SUBROUTINE SDIRK - -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - CONTAINS ! PROCEDURES INTERNAL TO SDIRK -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Ierr ) -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - USE KPP_ROOT_Parameters - IMPLICIT NONE - -!~~~> Arguments: - INTEGER, INTENT(IN) :: N - KPP_REAL, INTENT(INOUT) :: Y(NVAR) - KPP_REAL, INTENT(IN) :: Tinitial, Tfinal - INTEGER, INTENT(OUT) :: Ierr - -!~~~> Local variables: - KPP_REAL :: Z(NVAR,Smax), G(NVAR), TMP(NVAR), & - NewtonRate, SCAL(NVAR), RHS(NVAR), & - T, H, Theta, Hratio, NewtonPredictedErr, & - Qnewton, Err, Fac, Hnew,Tdirection, & - NewtonIncrement, NewtonIncrementOld - INTEGER :: j, IER, istage, NewtonIter, IP(NVAR) - LOGICAL :: Reject, FirstStep, SkipJac, SkipLU, NewtonDone - -#ifdef FULL_ALGEBRA - KPP_REAL FJAC(NVAR,NVAR), E(NVAR,NVAR) -#else - KPP_REAL FJAC(LU_NONZERO), E(LU_NONZERO) -#endif - KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0 - - -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -!~~~> Initializations -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - T = Tinitial - Tdirection = SIGN(ONE,Tfinal-Tinitial) - H = MAX(ABS(Hmin),ABS(Hstart)) - IF (ABS(H) <= 10.D0*Roundoff) H=1.0D-6 - H=MIN(ABS(H),Hmax) - H=SIGN(H,Tdirection) - SkipLU =.FALSE. - SkipJac = .FALSE. - Reject=.FALSE. - FirstStep=.TRUE. - - CALL SDIRK_ErrorScale(N, ITOL, AbsTol, RelTol, Y, SCAL) - -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -!~~~> Time loop begins -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -Tloop: DO WHILE ( (Tfinal-T)*Tdirection - Roundoff > ZERO ) - - -!~~~> Compute E = 1/(h*gamma)-Jac and its LU decomposition - IF ( .NOT.SkipLU ) THEN ! This time around skip the Jac update and LU - CALL SDIRK_PrepareMatrix ( H, T, Y, FJAC, & - SkipJac, SkipLU, E, IP, Reject, IER ) - IF (IER /= 0) THEN - CALL SDIRK_ErrorMsg(-8,T,H,Ierr); RETURN - END IF - END IF - - IF (ISTATUS(Nstp) > Max_no_steps) THEN - CALL SDIRK_ErrorMsg(-6,T,H,Ierr); RETURN - END IF - IF ( (T+0.1d0*H == T) .OR. (ABS(H) <= Roundoff) ) THEN - CALL SDIRK_ErrorMsg(-7,T,H,Ierr); RETURN - END IF - -stages:DO istage = 1, rkS - -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -!~~~> Simplified Newton iterations -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -!~~~> Starting values for Newton iterations - CALL Set2zero(N,Z(1,istage)) - -!~~~> Prepare the loop-independent part of the right-hand side - CALL Set2zero(N,G) - IF (istage > 1) THEN - DO j = 1, istage-1 - ! Gj(:) = sum_j Theta(i,j)*Zj(:) = H * sum_j A(i,j)*Fun(Zj) - CALL WAXPY(N,rkTheta(istage,j),Z(1,j),1,G,1) - ! Zi(:) = sum_j Alpha(i,j)*Zj(:) - IF (StartNewton) THEN - CALL WAXPY(N,rkAlpha(istage,j),Z(1,j),1,Z(1,istage),1) - END IF - END DO - END IF - - !~~~> Initializations for Newton iteration - NewtonDone = .FALSE. - Fac = 0.5d0 ! Step reduction factor if too many iterations - -NewtonLoop:DO NewtonIter = 1, NewtonMaxit - -!~~~> Prepare the loop-dependent part of the right-hand side - CALL WADD(N,Y,Z(1,istage),TMP) ! TMP <- Y + Zi - CALL FUN_CHEM(T+rkC(istage)*H,TMP,RHS) ! RHS <- Fun(Y+Zi) - ISTATUS(Nfun) = ISTATUS(Nfun) + 1 -! RHS(1:N) = G(1:N) - Z(1:N,istage) + (H*rkGamma)*RHS(1:N) - CALL WSCAL(N, H*rkGamma, RHS, 1) - CALL WAXPY (N, -ONE, Z(1,istage), 1, RHS, 1) - CALL WAXPY (N, ONE, G,1, RHS,1) - -!~~~> Solve the linear system - CALL SDIRK_Solve ( H, N, E, IP, IER, RHS ) - -!~~~> Check convergence of Newton iterations - CALL SDIRK_ErrorNorm(N, RHS, SCAL, NewtonIncrement) - IF ( NewtonIter == 1 ) THEN - Theta = ABS(ThetaMin) - NewtonRate = 2.0d0 - ELSE - Theta = NewtonIncrement/NewtonIncrementOld - IF (Theta < 0.99d0) THEN - NewtonRate = Theta/(ONE-Theta) - ! Predict error at the end of Newton process - NewtonPredictedErr = NewtonIncrement & - *Theta**(NewtonMaxit-NewtonIter)/(ONE-Theta) - IF (NewtonPredictedErr >= NewtonTol) THEN - ! Non-convergence of Newton: predicted error too large - Qnewton = MIN(10.0d0,NewtonPredictedErr/NewtonTol) - Fac = 0.8d0*Qnewton**(-ONE/(1+NewtonMaxit-NewtonIter)) - EXIT NewtonLoop - END IF - ELSE ! Non-convergence of Newton: Theta too large - EXIT NewtonLoop - END IF - END IF - NewtonIncrementOld = NewtonIncrement - ! Update solution: Z(:) <-- Z(:)+RHS(:) - CALL WAXPY(N,ONE,RHS,1,Z(1,istage),1) - - ! Check error in Newton iterations - NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol) - IF (NewtonDone) EXIT NewtonLoop - - END DO NewtonLoop - - IF (.NOT.NewtonDone) THEN - !CALL RK_ErrorMsg(-12,T,H,Ierr); - H = Fac*H; Reject=.TRUE. - SkipJac = .TRUE.; SkipLU = .FALSE. - CYCLE Tloop - END IF - -!~~~> End of simplified Newton iterations - - - END DO stages - - -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -!~~~> Error estimation -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - ISTATUS(Nstp) = ISTATUS(Nstp) + 1 - - IF (sdMethod /= BEL) THEN ! All methods but Backward Euler - CALL Set2zero(N,TMP) - DO i = 1,rkS - IF (rkE(i)/=ZERO) CALL WAXPY(N,rkE(i),Z(1,i),1,TMP,1) - END DO - - CALL SDIRK_Solve( H, N, E, IP, IER, TMP ) - CALL SDIRK_ErrorNorm(N, TMP, SCAL, Err) - -!~~~> Computation of new step size Hnew - Fac = FacSafe*(Err)**(-ONE/rkELO) - Fac = MAX(FacMin,MIN(FacMax,Fac)) - Hnew = H*Fac - ELSE ! Backward Euler Method - Err = ZERO - Hnew = MIN( Hstart, 1.2_dp*H ) - END IF - -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -!~~~> Accept/Reject step -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -accept: IF ( Err < ONE ) THEN !~~~> Step is accepted - - FirstStep=.FALSE. - ISTATUS(Nacc) = ISTATUS(Nacc) + 1 - -!~~~> Update time and solution - T = T + H - ! Y(:) <-- Y(:) + Sum_j rkD(j)*Z_j(:) - DO i = 1,rkS - IF (rkD(i)/=ZERO) CALL WAXPY(N,rkD(i),Z(1,i),1,Y,1) - END DO - -!~~~> Update scaling coefficients - CALL SDIRK_ErrorScale(N, ITOL, AbsTol, RelTol, Y, SCAL) - -!~~~> Next time step - Hnew = Tdirection*MIN(ABS(Hnew),Hmax) - ! Last T and H - RSTATUS(Ntexit) = T - RSTATUS(Nhexit) = H - RSTATUS(Nhnew) = Hnew - ! No step increase after a rejection - IF (Reject) Hnew = Tdirection*MIN(ABS(Hnew),ABS(H)) - Reject = .FALSE. - IF ((T+Hnew/Qmin-Tfinal)*Tdirection > ZERO) THEN - H = Tfinal-T - ELSE - Hratio=Hnew/H - ! If step not changed too much keep Jacobian and reuse LU - SkipLU = ( (Theta <= ThetaMin) .AND. (Hratio >= Qmin) & - .AND. (Hratio <= Qmax) ) - IF (.NOT.SkipLU) H = Hnew - END IF - ! If convergence is fast enough, do not update Jacobian -! SkipJac = (Theta <= ThetaMin) - SkipJac = .FALSE. - - ELSE accept !~~~> Step is rejected - - IF (FirstStep .OR. Reject) THEN - H = FacRej*H - ELSE - H = Hnew - END IF - Reject = .TRUE. - SkipJac = .TRUE. - SkipLU = .FALSE. - IF (ISTATUS(Nacc) >= 1) ISTATUS(Nrej) = ISTATUS(Nrej) + 1 - - END IF accept - - END DO Tloop - - ! Successful return - Ierr = 1 - - END SUBROUTINE SDIRK_Integrator - - -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - SUBROUTINE SDIRK_ErrorScale(N, ITOL, AbsTol, RelTol, Y, SCAL) -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - IMPLICIT NONE - INTEGER :: i, N, ITOL - KPP_REAL :: AbsTol(NVAR), RelTol(NVAR), & - Y(NVAR), SCAL(NVAR) - IF (ITOL == 0) THEN - DO i=1,NVAR - SCAL(i) = ONE / ( AbsTol(1)+RelTol(1)*ABS(Y(i)) ) - END DO - ELSE - DO i=1,NVAR - SCAL(i) = ONE / ( AbsTol(i)+RelTol(i)*ABS(Y(i)) ) - END DO - END IF - END SUBROUTINE SDIRK_ErrorScale - - -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - SUBROUTINE SDIRK_ErrorNorm(N, Y, SCAL, Err) -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -! - INTEGER :: i, N - KPP_REAL :: Y(N), SCAL(N), Err - Err = ZERO - DO i=1,N - Err = Err+(Y(i)*SCAL(i))**2 - END DO - Err = MAX( SQRT(Err/DBLE(N)), 1.0d-10 ) -! - END SUBROUTINE SDIRK_ErrorNorm - - -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - SUBROUTINE SDIRK_ErrorMsg(Code,T,H,Ierr) -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -! Handles all error messages -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - KPP_REAL, INTENT(IN) :: T, H - INTEGER, INTENT(IN) :: Code - INTEGER, INTENT(OUT) :: Ierr - - Ierr = Code - PRINT * , & - 'Forced exit from SDIRK due to the following error:' - - SELECT CASE (Code) - CASE (-1) - PRINT * , '--> Improper value for maximal no of steps' - CASE (-2) - PRINT * , '--> Improper value for maximal no of Newton iterations' - CASE (-3) - PRINT * , '--> Hmin/Hmax/Hstart must be positive' - CASE (-4) - PRINT * , '--> FacMin/FacMax/FacRej must be positive' - CASE (-5) - PRINT * , '--> Improper tolerance values' - CASE (-6) - PRINT * , '--> No of steps exceeds maximum bound' - CASE (-7) - PRINT * , '--> Step size too small: T + 10*H = T', & - ' or H < Roundoff' - CASE (-8) - PRINT * , '--> Matrix is repeatedly singular' - CASE DEFAULT - PRINT *, 'Unknown Error code: ', Code - END SELECT - - PRINT *, "T=", T, "and H=", H - - END SUBROUTINE SDIRK_ErrorMsg - - -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - SUBROUTINE SDIRK_PrepareMatrix ( H, T, Y, FJAC, & - SkipJac, SkipLU, E, IP, Reject, ISING ) -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -!~~~> Compute the matrix E = 1/(H*GAMMA)*Jac, and its decomposition -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - IMPLICIT NONE - - KPP_REAL, INTENT(INOUT) :: H - KPP_REAL, INTENT(IN) :: T, Y(NVAR) - LOGICAL, INTENT(INOUT) :: SkipJac,SkipLU,Reject - INTEGER, INTENT(OUT) :: ISING, IP(NVAR) -#ifdef FULL_ALGEBRA - KPP_REAL, INTENT(INOUT) :: FJAC(NVAR,NVAR) - KPP_REAL, INTENT(OUT) :: E(NVAR,NVAR) -#else - KPP_REAL, INTENT(INOUT) :: FJAC(LU_NONZERO) - KPP_REAL, INTENT(OUT) :: E(LU_NONZERO) -#endif - KPP_REAL :: HGammaInv - INTEGER :: i, j, ConsecutiveSng - - ConsecutiveSng = 0 - ISING = 1 - -Hloop: DO WHILE (ISING /= 0) - - HGammaInv = ONE/(H*rkGamma) - -!~~~> Compute the Jacobian -! IF (SkipJac) THEN -! SkipJac = .FALSE. -! ELSE - IF (.NOT. SkipJac) THEN - CALL JAC_CHEM( T, Y, FJAC ) - ISTATUS(Njac) = ISTATUS(Njac) + 1 - END IF - -#ifdef FULL_ALGEBRA - DO j=1,NVAR - DO i=1,NVAR - E(i,j) = -FJAC(i,j) - END DO - E(j,j) = E(j,j)+HGammaInv - END DO - CALL DGETRF( NVAR, NVAR, E, NVAR, IP, ISING ) -#else - DO i = 1,LU_NONZERO - E(i) = -FJAC(i) - END DO - DO i = 1,NVAR - j = LU_DIAG(i); E(j) = E(j) + HGammaInv - END DO - CALL KppDecomp ( E, ISING) - IP(1) = 1 -#endif - ISTATUS(Ndec) = ISTATUS(Ndec) + 1 - - IF (ISING /= 0) THEN - WRITE (6,*) ' MATRIX IS SINGULAR, ISING=',ISING,' T=',T,' H=',H - ISTATUS(Nsng) = ISTATUS(Nsng) + 1; ConsecutiveSng = ConsecutiveSng + 1 - IF (ConsecutiveSng >= 6) RETURN ! Failure - H = 0.5d0*H - SkipJac = .TRUE. - SkipLU = .FALSE. - Reject = .TRUE. - END IF - - END DO Hloop - - END SUBROUTINE SDIRK_PrepareMatrix - - -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - SUBROUTINE SDIRK_Solve ( H, N, E, IP, ISING, RHS ) -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -!~~~> Solves the system (H*Gamma-Jac)*x = RHS -! using the LU decomposition of E = I - 1/(H*Gamma)*Jac -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - IMPLICIT NONE - INTEGER, INTENT(IN) :: N, IP(N), ISING - KPP_REAL, INTENT(IN) :: H -#ifdef FULL_ALGEBRA - KPP_REAL, INTENT(IN) :: E(NVAR,NVAR) -#else - KPP_REAL, INTENT(IN) :: E(LU_NONZERO) -#endif - KPP_REAL, INTENT(INOUT) :: RHS(N) - KPP_REAL :: HGammaInv - - HGammaInv = ONE/(H*rkGamma) - CALL WSCAL(N,HGammaInv,RHS,1) -#ifdef FULL_ALGEBRA - CALL DGETRS( 'N', N, 1, E, N, IP, RHS, N, ISING ) -#else - CALL KppSolve(E, RHS) -#endif - ISTATUS(Nsol) = ISTATUS(Nsol) + 1 - - END SUBROUTINE SDIRK_Solve - - - -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -SUBROUTINE BEuler -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -sdMethod = BEL -! Number of stages -rkS = 1 - -! Method coefficients -rkGamma = 1.0_dp -rkA(1,1) = 1.0_dp -rkB(1) = 1.0_dp -rkBhat(1)= 1.0_dp -rkC(1) = 1.0_dp - -! Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i} -rkD(1) = 1.0_dp - -! Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i} -rkE(1) = 1.0_dp - -! Local order of Err estimate -rkElo = -1.0_dp - -! h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j} -rkTheta(1,1) = 0.0_dp - -! Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j} -rkAlpha(1,1) = 0.0_dp - -END SUBROUTINE BEuler - - - -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - SUBROUTINE Sdirk4a -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - sdMethod = S4A -! Number of stages - rkS = 5 - -! Method coefficients - rkGamma = .2666666666666666666666666666666667d0 - - rkA(1,1) = .2666666666666666666666666666666667d0 - rkA(2,1) = .5000000000000000000000000000000000d0 - rkA(2,2) = .2666666666666666666666666666666667d0 - rkA(3,1) = .3541539528432732316227461858529820d0 - rkA(3,2) = -.5415395284327323162274618585298197d-1 - rkA(3,3) = .2666666666666666666666666666666667d0 - rkA(4,1) = .8515494131138652076337791881433756d-1 - rkA(4,2) = -.6484332287891555171683963466229754d-1 - rkA(4,3) = .7915325296404206392428857585141242d-1 - rkA(4,4) = .2666666666666666666666666666666667d0 - rkA(5,1) = 2.100115700566932777970612055999074d0 - rkA(5,2) = -.7677800284445976813343102185062276d0 - rkA(5,3) = 2.399816361080026398094746205273880d0 - rkA(5,4) = -2.998818699869028161397714709433394d0 - rkA(5,5) = .2666666666666666666666666666666667d0 - - rkB(1) = 2.100115700566932777970612055999074d0 - rkB(2) = -.7677800284445976813343102185062276d0 - rkB(3) = 2.399816361080026398094746205273880d0 - rkB(4) = -2.998818699869028161397714709433394d0 - rkB(5) = .2666666666666666666666666666666667d0 - - rkBhat(1)= 2.885264204387193942183851612883390d0 - rkBhat(2)= -.1458793482962771337341223443218041d0 - rkBhat(3)= 2.390008682465139866479830743628554d0 - rkBhat(4)= -4.129393538556056674929560012190140d0 - rkBhat(5)= 0.d0 - - rkC(1) = .2666666666666666666666666666666667d0 - rkC(2) = .7666666666666666666666666666666667d0 - rkC(3) = .5666666666666666666666666666666667d0 - rkC(4) = .3661315380631796996374935266701191d0 - rkC(5) = 1.d0 - -! Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i} - rkD(1) = 0.d0 - rkD(2) = 0.d0 - rkD(3) = 0.d0 - rkD(4) = 0.d0 - rkD(5) = 1.d0 - -! Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i} - rkE(1) = -.6804000050475287124787034884002302d0 - rkE(2) = 1.558961944525217193393931795738823d0 - rkE(3) = -13.55893003128907927748632408763868d0 - rkE(4) = 15.48522576958521253098585004571302d0 - rkE(5) = 1.d0 - -! Local order of Err estimate - rkElo = 4 - -! h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j} - rkTheta(2,1) = 1.875000000000000000000000000000000d0 - rkTheta(3,1) = 1.708847304091539528432732316227462d0 - rkTheta(3,2) = -.2030773231622746185852981969486824d0 - rkTheta(4,1) = .2680325578937783958847157206823118d0 - rkTheta(4,2) = -.1828840955527181631794050728644549d0 - rkTheta(4,3) = .2968246986151577397160821594427966d0 - rkTheta(5,1) = .9096171815241460655379433581446771d0 - rkTheta(5,2) = -3.108254967778352416114774430509465d0 - rkTheta(5,3) = 12.33727431701306195581826123274001d0 - rkTheta(5,4) = -11.24557012450885560524143016037523d0 - -! Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j} - rkAlpha(2,1) = 2.875000000000000000000000000000000d0 - rkAlpha(3,1) = .8500000000000000000000000000000000d0 - rkAlpha(3,2) = .4434782608695652173913043478260870d0 - rkAlpha(4,1) = .7352046091658870564637910527807370d0 - rkAlpha(4,2) = -.9525565003057343527941920657462074d-1 - rkAlpha(4,3) = .4290111305453813852259481840631738d0 - rkAlpha(5,1) = -16.10898993405067684831655675112808d0 - rkAlpha(5,2) = 6.559571569643355712998131800797873d0 - rkAlpha(5,3) = -15.90772144271326504260996815012482d0 - rkAlpha(5,4) = 25.34908987169226073668861694892683d0 - -!~~~> Coefficients for continuous solution -! rkD(1,1)= 24.74416644927758d0 -! rkD(1,2)= -4.325375951824688d0 -! rkD(1,3)= 41.39683763286316d0 -! rkD(1,4)= -61.04144619901784d0 -! rkD(1,5)= -3.391332232917013d0 -! rkD(2,1)= -51.98245719616925d0 -! rkD(2,2)= 10.52501981094525d0 -! rkD(2,3)= -154.2067922191855d0 -! rkD(2,4)= 214.3082125319825d0 -! rkD(2,5)= 14.71166018088679d0 -! rkD(3,1)= 33.14347947522142d0 -! rkD(3,2)= -19.72986789558523d0 -! rkD(3,3)= 230.4878502285804d0 -! rkD(3,4)= -287.6629744338197d0 -! rkD(3,5)= -18.99932366302254d0 -! rkD(4,1)= -5.905188728329743d0 -! rkD(4,2)= 13.53022403646467d0 -! rkD(4,3)= -117.6778956422581d0 -! rkD(4,4)= 134.3962081008550d0 -! rkD(4,5)= 8.678995715052762d0 -! -! DO i=1,4 ! CONTi <-- Sum_j rkD(i,j)*Zj -! CALL Set2zero(N,CONT(1,i)) -! DO j = 1,rkS -! CALL WAXPY(N,rkD(i,j),Z(1,j),1,CONT(1,i),1) -! END DO -! END DO - - rkELO = 4.0d0 - - END SUBROUTINE Sdirk4a - -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - SUBROUTINE Sdirk4b -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - sdMethod = S4B -! Number of stages - rkS = 5 - -! Method coefficients - rkGamma = .25d0 - - rkA(1,1) = 0.25d0 - rkA(2,1) = 0.5d00 - rkA(2,2) = 0.25d0 - rkA(3,1) = 0.34d0 - rkA(3,2) =-0.40d-1 - rkA(3,3) = 0.25d0 - rkA(4,1) = 0.2727941176470588235294117647058824d0 - rkA(4,2) =-0.5036764705882352941176470588235294d-1 - rkA(4,3) = 0.2757352941176470588235294117647059d-1 - rkA(4,4) = 0.25d0 - rkA(5,1) = 1.041666666666666666666666666666667d0 - rkA(5,2) =-1.020833333333333333333333333333333d0 - rkA(5,3) = 7.812500000000000000000000000000000d0 - rkA(5,4) =-7.083333333333333333333333333333333d0 - rkA(5,5) = 0.25d0 - - rkB(1) = 1.041666666666666666666666666666667d0 - rkB(2) = -1.020833333333333333333333333333333d0 - rkB(3) = 7.812500000000000000000000000000000d0 - rkB(4) = -7.083333333333333333333333333333333d0 - rkB(5) = 0.250000000000000000000000000000000d0 - - rkBhat(1)= 1.069791666666666666666666666666667d0 - rkBhat(2)= -0.894270833333333333333333333333333d0 - rkBhat(3)= 7.695312500000000000000000000000000d0 - rkBhat(4)= -7.083333333333333333333333333333333d0 - rkBhat(5)= 0.212500000000000000000000000000000d0 - - rkC(1) = 0.25d0 - rkC(2) = 0.75d0 - rkC(3) = 0.55d0 - rkC(4) = 0.50d0 - rkC(5) = 1.00d0 - -! Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i} - rkD(1) = 0.0d0 - rkD(2) = 0.0d0 - rkD(3) = 0.0d0 - rkD(4) = 0.0d0 - rkD(5) = 1.0d0 - -! Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i} - rkE(1) = 0.5750d0 - rkE(2) = 0.2125d0 - rkE(3) = -4.6875d0 - rkE(4) = 4.2500d0 - rkE(5) = 0.1500d0 - -! Local order of Err estimate - rkElo = 4 - -! h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j} - rkTheta(2,1) = 2.d0 - rkTheta(3,1) = 1.680000000000000000000000000000000d0 - rkTheta(3,2) = -.1600000000000000000000000000000000d0 - rkTheta(4,1) = 1.308823529411764705882352941176471d0 - rkTheta(4,2) = -.1838235294117647058823529411764706d0 - rkTheta(4,3) = 0.1102941176470588235294117647058824d0 - rkTheta(5,1) = -3.083333333333333333333333333333333d0 - rkTheta(5,2) = -4.291666666666666666666666666666667d0 - rkTheta(5,3) = 34.37500000000000000000000000000000d0 - rkTheta(5,4) = -28.33333333333333333333333333333333d0 - -! Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j} - rkAlpha(2,1) = 3. - rkAlpha(3,1) = .8800000000000000000000000000000000d0 - rkAlpha(3,2) = .4400000000000000000000000000000000d0 - rkAlpha(4,1) = .1666666666666666666666666666666667d0 - rkAlpha(4,2) = -.8333333333333333333333333333333333d-1 - rkAlpha(4,3) = .9469696969696969696969696969696970d0 - rkAlpha(5,1) = -6.d0 - rkAlpha(5,2) = 9.d0 - rkAlpha(5,3) = -56.81818181818181818181818181818182d0 - rkAlpha(5,4) = 54.d0 - - END SUBROUTINE Sdirk4b - -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - SUBROUTINE Sdirk2a -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - sdMethod = S2A -! Number of stages - rkS = 2 - -! Method coefficients - rkGamma = .2928932188134524755991556378951510d0 - - rkA(1,1) = .2928932188134524755991556378951510d0 - rkA(2,1) = .7071067811865475244008443621048490d0 - rkA(2,2) = .2928932188134524755991556378951510d0 - - rkB(1) = .7071067811865475244008443621048490d0 - rkB(2) = .2928932188134524755991556378951510d0 - - rkBhat(1)= .6666666666666666666666666666666667d0 - rkBhat(2)= .3333333333333333333333333333333333d0 - - rkC(1) = 0.292893218813452475599155637895151d0 - rkC(2) = 1.0d0 - -! Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i} - rkD(1) = 0.0d0 - rkD(2) = 1.0d0 - -! Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i} - rkE(1) = 0.4714045207910316829338962414032326d0 - rkE(2) = -0.1380711874576983496005629080698993d0 - -! Local order of Err estimate - rkElo = 2 - -! h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j} - rkTheta(2,1) = 2.414213562373095048801688724209698d0 - -! Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j} - rkAlpha(2,1) = 3.414213562373095048801688724209698d0 - - END SUBROUTINE Sdirk2a - -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - SUBROUTINE Sdirk2b -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - sdMethod = S2B -! Number of stages - rkS = 2 - -! Method coefficients - rkGamma = 1.707106781186547524400844362104849d0 - - rkA(1,1) = 1.707106781186547524400844362104849d0 - rkA(2,1) = -.707106781186547524400844362104849d0 - rkA(2,2) = 1.707106781186547524400844362104849d0 - - rkB(1) = -.707106781186547524400844362104849d0 - rkB(2) = 1.707106781186547524400844362104849d0 - - rkBhat(1)= .6666666666666666666666666666666667d0 - rkBhat(2)= .3333333333333333333333333333333333d0 - - rkC(1) = 1.707106781186547524400844362104849d0 - rkC(2) = 1.0d0 - -! Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i} - rkD(1) = 0.0d0 - rkD(2) = 1.0d0 - -! Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i} - rkE(1) = -.4714045207910316829338962414032326d0 - rkE(2) = .8047378541243650162672295747365659d0 - -! Local order of Err estimate - rkElo = 2 - -! h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j} - rkTheta(2,1) = -.414213562373095048801688724209698d0 - -! Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j} - rkAlpha(2,1) = .5857864376269049511983112757903019d0 - - END SUBROUTINE Sdirk2b - - -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - SUBROUTINE Sdirk3a -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - sdMethod = S3A -! Number of stages - rkS = 3 - -! Method coefficients - rkGamma = .2113248654051871177454256097490213d0 - - rkA(1,1) = .2113248654051871177454256097490213d0 - rkA(2,1) = .2113248654051871177454256097490213d0 - rkA(2,2) = .2113248654051871177454256097490213d0 - rkA(3,1) = .2113248654051871177454256097490213d0 - rkA(3,2) = .5773502691896257645091487805019573d0 - rkA(3,3) = .2113248654051871177454256097490213d0 - - rkB(1) = .2113248654051871177454256097490213d0 - rkB(2) = .5773502691896257645091487805019573d0 - rkB(3) = .2113248654051871177454256097490213d0 - - rkBhat(1)= .2113248654051871177454256097490213d0 - rkBhat(2)= .6477918909913548037576239837516312d0 - rkBhat(3)= .1408832436034580784969504064993475d0 - - rkC(1) = .2113248654051871177454256097490213d0 - rkC(2) = .4226497308103742354908512194980427d0 - rkC(3) = 1.d0 - -! Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i} - rkD(1) = 0.d0 - rkD(2) = 0.d0 - rkD(3) = 1.d0 - -! Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i} - rkE(1) = 0.9106836025229590978424821138352906d0 - rkE(2) = -1.244016935856292431175815447168624d0 - rkE(3) = 0.3333333333333333333333333333333333d0 - -! Local order of Err estimate - rkElo = 2 - -! h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j} - rkTheta(2,1) = 1.0d0 - rkTheta(3,1) = -1.732050807568877293527446341505872d0 - rkTheta(3,2) = 2.732050807568877293527446341505872d0 - -! Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j} - rkAlpha(2,1) = 2.0d0 - rkAlpha(3,1) = -12.92820323027550917410978536602349d0 - rkAlpha(3,2) = 8.83012701892219323381861585376468d0 - - END SUBROUTINE Sdirk3a - -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - END SUBROUTINE SDIRK ! AND ALL ITS INTERNAL PROCEDURES -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - SUBROUTINE FUN_CHEM( T, Y, P ) -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - USE KPP_ROOT_Parameters, ONLY: NVAR, LU_NONZERO - USE KPP_ROOT_Global, ONLY: FIX, RCONST, TIME - USE KPP_ROOT_Function, ONLY: Fun - USE KPP_ROOT_Rates, ONLY: Update_SUN, Update_RCONST, Update_PHOTO - IMPLICIT NONE - - KPP_REAL :: T, Told - KPP_REAL :: Y(NVAR), P(NVAR) - - Told = TIME - TIME = T - IF ( Do_Update_SUN ) CALL Update_SUN() - IF ( Do_Update_RCONST ) CALL Update_RCONST() - - CALL Fun( Y, FIX, RCONST, P ) - - TIME = Told - - END SUBROUTINE FUN_CHEM - - -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - SUBROUTINE JAC_CHEM( T, Y, JV ) -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - USE KPP_ROOT_Parameters, ONLY: NVAR, LU_NONZERO - USE KPP_ROOT_Global, ONLY: FIX, RCONST, TIME - USE KPP_ROOT_Jacobian - USE KPP_ROOT_Jacobian, ONLY: Jac_SP - USE KPP_ROOT_Rates, ONLY: Update_SUN, Update_RCONST, Update_PHOTO - IMPLICIT NONE - - KPP_REAL :: T, Told - KPP_REAL :: Y(NVAR) -#ifdef FULL_ALGEBRA - KPP_REAL :: JS(LU_NONZERO), JV(NVAR,NVAR) - INTEGER :: i, j -#else - KPP_REAL :: JV(LU_NONZERO) -#endif - - Told = TIME - TIME = T - IF ( Do_Update_SUN ) CALL Update_SUN() - IF ( Do_Update_RCONST ) CALL Update_RCONST() - -#ifdef FULL_ALGEBRA - CALL Jac_SP(Y, FIX, RCONST, JS) - DO j=1,NVAR - DO i=1,NVAR - JV(i,j) = 0.0D0 - END DO - END DO - DO i=1,LU_NONZERO - JV(LU_IROW(i),LU_ICOL(i)) = JS(i) - END DO -#else - CALL Jac_SP(Y, FIX, RCONST, JV) -#endif - TIME = Told - - END SUBROUTINE JAC_CHEM - -END MODULE KPP_ROOT_Integrator - - diff --git a/int/sdirk.f90 b/int/sdirk.f90 index f50d3d3..a43af72 100644 --- a/int/sdirk.f90 +++ b/int/sdirk.f90 @@ -3,6 +3,8 @@ ! * Sdirk 2a, 2b: L-stable, 2 stages, order 2 ! ! * Sdirk 3a: L-stable, 3 stages, order 2, adj-invariant ! ! * Sdirk 4a, 4b: L-stable, 5 stages, order 4 ! +! * Backward Euler: L-stable, 1 stage, order 1 ! +! Fixed step size = RCNTRL(3) = Hstart ! ! By default the code employs the KPP sparse linear algebra routines ! ! Compile with -DFULL_ALGEBRA to use full linear algebra (LAPACK) ! ! ! @@ -10,6 +12,7 @@ ! Virginia Polytechnic Institute and State University ! ! Contact: sandu@cs.vt.edu ! ! Revised by Philipp Miehe and Adrian Sandu, May 2006 ! +! Backward Euler added by Adrian Sandu, December 2012 ! ! This implementation is part of KPP - the Kinetic PreProcessor ! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! @@ -18,9 +21,10 @@ MODULE KPP_ROOT_Integrator USE KPP_ROOT_Precision USE KPP_ROOT_Global USE KPP_ROOT_Parameters - USE KPP_ROOT_JacobianSP , ONLY: LU_DIAG - USE KPP_ROOT_LinearAlgebra, ONLY: KppDecomp, & - KppSolve, Set2zero, WLAMCH, WCOPY, WAXPY, WSCAL, WADD + USE KPP_ROOT_JacobianSP, ONLY : LU_DIAG + USE KPP_ROOT_LinearAlgebra, ONLY : KppDecomp, KppSolve, Set2zero, & + WLAMCH, WCOPY, WAXPY, & + WSCAL, WADD IMPLICIT NONE PUBLIC @@ -38,6 +42,14 @@ MODULE KPP_ROOT_Integrator Nrej=5, Ndec=6, Nsol=7, Nsng=8, & Ntexit=1, Nhexit=2, Nhnew=3 + ! mz_rs_20171120+ + ! description of the error numbers IERR + CHARACTER(LEN=50), PARAMETER, DIMENSION(-1:1) :: IERR_NAMES = (/ & + 'DUMMY ERROR MESSAGE ', & ! -1 + ' ', & ! 0 (not used) + 'Success ' /) ! 1 + ! mz_rs_20171120- + CONTAINS SUBROUTINE INTEGRATE( TIN, TOUT, ICNTRL_U, RCNTRL_U, & @@ -49,7 +61,7 @@ SUBROUTINE INTEGRATE( TIN, TOUT, ICNTRL_U, RCNTRL_U, & KPP_REAL, INTENT(IN) :: TIN ! Start Time KPP_REAL, INTENT(IN) :: TOUT ! End Time - !~~~> Optional input parameters and statistics + ! Optional input parameters and statistics INTEGER, INTENT(IN), OPTIONAL :: ICNTRL_U(20) KPP_REAL, INTENT(IN), OPTIONAL :: RCNTRL_U(20) INTEGER, INTENT(OUT), OPTIONAL :: ISTATUS_U(20) @@ -73,14 +85,14 @@ SUBROUTINE INTEGRATE( TIN, TOUT, ICNTRL_U, RCNTRL_U, & ICNTRL(15) = 5 ! Call Update_SUN and Update_RCONST from w/in the int. !~~~> if optional parameters are given, and if they are /= 0, - !~~~> then use them to overwrite default settings + ! then use them to overwrite default settings IF ( PRESENT( ICNTRL_U ) ) THEN WHERE( ICNTRL_U /= 0 ) ICNTRL = ICNTRL_U ENDIF IF ( PRESENT( RCNTRL_U ) ) THEN WHERE( RCNTRL_U > 0 ) RCNTRL = RCNTRL_U ENDIF - + !~~~> Determine the settings of the Do_Update_* flags, which determine !~~~> whether or not we need to call Update_* routines in the integrator !~~~> (or not, if we are calling them from a higher-level) @@ -97,7 +109,7 @@ SUBROUTINE INTEGRATE( TIN, TOUT, ICNTRL_U, RCNTRL_U, & Do_Update_RCONST, & Do_Update_PHOTO, & Do_Update_Sun ) - + !~~~> In order to remove the prior EQUIVALENCE statements (which !~~~> are not thread-safe), we now have declared VAR and FIX as !~~~> threadprivate pointer variables that can point to C. @@ -106,8 +118,8 @@ SUBROUTINE INTEGRATE( TIN, TOUT, ICNTRL_U, RCNTRL_U, & !~~~> Call the integrator T1 = TIN; T2 = TOUT - CALL SDIRK( NVAR, T1, T2, VAR, RTOL, ATOL, & - RCNTRL, ICNTRL, RSTATUS,ISTATUS, Ierr ) + CALL SDIRK( NVAR,T1,T2,VAR,RTOL,ATOL, & + RCNTRL,ICNTRL,RSTATUS,ISTATUS,Ierr ) !~~~> Free pointers VAR => NULL() @@ -128,7 +140,7 @@ SUBROUTINE INTEGRATE( TIN, TOUT, ICNTRL_U, RCNTRL_U, & IF ( PRESENT( RSTATUS_U ) ) RSTATUS_U = RSTATUS IF ( PRESENT( IERR_U ) ) IERR_U = IERR -END SUBROUTINE INTEGRATE + END SUBROUTINE INTEGRATE !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -150,11 +162,14 @@ SUBROUTINE SDIRK(N, Tinitial, Tfinal, Y, RelTol, AbsTol, & ! * Sdirk 2a, 2b: L-stable, 2 stages, order 2 ! * Sdirk 3a: L-stable, 3 stages, order 2, adjoint-invariant ! * Sdirk 4a, 4b: L-stable, 5 stages, order 4 +! * Backward Euler: L-stable, 1 stage, order 1 +! Fixed step size = RCNTRL(3) = Hstart ! ! (C) Adrian Sandu, July 2005 ! Virginia Polytechnic Institute and State University ! Contact: sandu@cs.vt.edu -! Revised by Philipp Miehe and Adrian Sandu, May 2006 +! Revised by Philipp Miehe and Adrian Sandu, May 2006 +! Backward Euler added by Adrian Sandu, Dec. 2012 ! This implementation is part of KPP - the Kinetic PreProcessor !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! @@ -216,7 +231,7 @@ SUBROUTINE SDIRK(N, Tinitial, Tfinal, Y, RelTol, AbsTol, & ! = 4 : Call Update_SUN from within the integrator ! = 5 : Call Update_SUN and Update_RCONST from within the int. ! = 6 : Call Update_SUN and Update_PHOTO from within the int. -! = 7 : Call Update_SUN, Update_PHOTO and Update_RCONST w/in the int. +! = 7 : Call Update_SUN, Update_PHOTO, Update_RCONST from w/in the int. ! !~~~> Real parameters ! @@ -278,7 +293,7 @@ SUBROUTINE SDIRK(N, Tinitial, Tfinal, Y, RelTol, AbsTol, & !~~~> SDIRK method coefficients, up to 5 stages INTEGER, PARAMETER :: Smax = 5 - INTEGER, PARAMETER :: S2A=1, S2B=2, S3A=3, S4A=4, S4B=5 + INTEGER, PARAMETER :: S2A=1, S2B=2, S3A=3, S4A=4, S4B=5, BEL=6 KPP_REAL :: rkGamma, rkA(Smax,Smax), rkB(Smax), rkC(Smax), & rkD(Smax), rkE(Smax), rkBhat(Smax), rkELO, & rkAlpha(Smax,Smax), rkTheta(Smax,Smax) @@ -307,17 +322,19 @@ SUBROUTINE SDIRK(N, Tinitial, Tfinal, Y, RelTol, AbsTol, & !~~~> ICNTRL(3) - method selection SELECT CASE (ICNTRL(3)) CASE (0,1) - CALL Sdirk2a + CALL Sdirk2a() CASE (2) - CALL Sdirk2b + CALL Sdirk2b() CASE (3) - CALL Sdirk3a + CALL Sdirk3a() CASE (4) - CALL Sdirk4a + CALL Sdirk4a() CASE (5) - CALL Sdirk4b + CALL Sdirk4b() + CASE (6) + CALL BEuler() CASE DEFAULT - CALL Sdirk2a + CALL Sdirk2a() END SELECT !~~~> The maximum number of time steps admitted @@ -373,7 +390,7 @@ SUBROUTINE SDIRK(N, Tinitial, Tfinal, Y, RelTol, AbsTol, & !~~~> Starting step size: (positive value) IF (RCNTRL(3) == ZERO) THEN - Hstart = MAX(Hmin,Roundoff) + Hstart = MAX(Hmin,101.*Roundoff) ! factor 101 added mz_rs_20171120 ELSEIF (RCNTRL(3) > ZERO) THEN Hstart = MIN(ABS(RCNTRL(3)),ABS(Tfinal-Tinitial)) ELSE @@ -460,6 +477,24 @@ SUBROUTINE SDIRK(N, Tinitial, Tfinal, Y, RelTol, AbsTol, & END DO END IF +!~~~> Special treatment of backward Euler + IF (sdMethod == BEL) THEN + StartNewton = .FALSE. + FacMin = 1.0_dp + FacMax = 1.0_dp + FacSafe= 1.0_dp + FacRej = 1.0_dp + Qmin = 1.0_dp + Qmax = 1.0_dp + Hmax = Hstart + IF ( Hstart <= 100_dp*Roundoff ) THEN + PRINT * , ' Backward Euler calling error:' + PRINT * , ' it requires Hstart > 100*Roundoff' + STOP + END IF + END IF + + IF (Ierr < 0) RETURN CALL SDIRK_Integrator( N,Tinitial,Tfinal,Y,Ierr ) @@ -620,7 +655,7 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Ierr ) CYCLE Tloop END IF -!~~~> End of implified Newton iterations +!~~~> End of simplified Newton iterations END DO stages @@ -630,18 +665,24 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Ierr ) !~~~> Error estimation !~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ISTATUS(Nstp) = ISTATUS(Nstp) + 1 - CALL Set2zero(N,TMP) - DO i = 1,rkS - IF (rkE(i)/=ZERO) CALL WAXPY(N,rkE(i),Z(1,i),1,TMP,1) - END DO - CALL SDIRK_Solve( H, N, E, IP, IER, TMP ) - CALL SDIRK_ErrorNorm(N, TMP, SCAL, Err) + IF (sdMethod /= BEL) THEN ! All methods but Backward Euler + CALL Set2zero(N,TMP) + DO i = 1,rkS + IF (rkE(i)/=ZERO) CALL WAXPY(N,rkE(i),Z(1,i),1,TMP,1) + END DO + + CALL SDIRK_Solve( H, N, E, IP, IER, TMP ) + CALL SDIRK_ErrorNorm(N, TMP, SCAL, Err) !~~~> Computation of new step size Hnew - Fac = FacSafe*(Err)**(-ONE/rkELO) - Fac = MAX(FacMin,MIN(FacMax,Fac)) - Hnew = H*Fac + Fac = FacSafe*(Err)**(-ONE/rkELO) + Fac = MAX(FacMin,MIN(FacMax,Fac)) + Hnew = H*Fac + ELSE ! Backward Euler Method + Err = ZERO + Hnew = MIN( Hstart, 1.2_dp*H ) + END IF !~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> Accept/Reject step @@ -883,8 +924,41 @@ SUBROUTINE SDIRK_Solve ( H, N, E, IP, ISING, RHS ) END SUBROUTINE SDIRK_Solve + +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + SUBROUTINE BEuler() +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + sdMethod = BEL + ! Number of stages + rkS = 1 + + ! Method coefficients + rkGamma = 1.0_dp + rkA(1,1) = 1.0_dp + rkB(1) = 1.0_dp + rkBhat(1)= 1.0_dp + rkC(1) = 1.0_dp + + ! Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i} + rkD(1) = 1.0_dp + + ! Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i} + rkE(1) = 1.0_dp + + ! Local order of Err estimate + rkElo = -1.0_dp + + ! h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j} + rkTheta(1,1) = 0.0_dp + + ! Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j} + rkAlpha(1,1) = 0.0_dp + + END SUBROUTINE BEuler + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - SUBROUTINE Sdirk4a + SUBROUTINE Sdirk4a() !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ sdMethod = S4A ! Number of stages @@ -1002,7 +1076,7 @@ SUBROUTINE Sdirk4a END SUBROUTINE Sdirk4a !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - SUBROUTINE Sdirk4b + SUBROUTINE Sdirk4b() !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ sdMethod = S4B @@ -1090,7 +1164,7 @@ SUBROUTINE Sdirk4b END SUBROUTINE Sdirk4b !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - SUBROUTINE Sdirk2a + SUBROUTINE Sdirk2a() !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ sdMethod = S2A @@ -1133,7 +1207,7 @@ SUBROUTINE Sdirk2a END SUBROUTINE Sdirk2a !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - SUBROUTINE Sdirk2b + SUBROUTINE Sdirk2b() !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ sdMethod = S2B @@ -1177,7 +1251,7 @@ END SUBROUTINE Sdirk2b !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - SUBROUTINE Sdirk3a + SUBROUTINE Sdirk3a() !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ sdMethod = S3A @@ -1240,11 +1314,10 @@ END SUBROUTINE SDIRK ! AND ALL ITS INTERNAL PROCEDURES SUBROUTINE FUN_CHEM( T, Y, P ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - USE KPP_ROOT_Parameters, ONLY : NVAR, LU_NONZERO - USE KPP_ROOT_Global, ONLY : FIX, RCONST, TIME - USE KPP_ROOT_Function, ONLY : Fun - USE KPP_ROOT_Rates, ONLY : Update_SUN, Update_RCONST - + USE KPP_ROOT_Parameters, ONLY: NVAR, LU_NONZERO + USE KPP_ROOT_Global, ONLY: FIX, RCONST, TIME + USE KPP_ROOT_Function, ONLY: Fun + USE KPP_ROOT_Rates, ONLY: Update_SUN, Update_RCONST, Update_PHOTO IMPLICIT NONE KPP_REAL :: T, Told @@ -1266,12 +1339,11 @@ END SUBROUTINE FUN_CHEM SUBROUTINE JAC_CHEM( T, Y, JV ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - USE KPP_ROOT_Parameters, ONLY : NVAR, LU_NONZERO - USE KPP_ROOT_Global, ONLY : FIX, RCONST, TIME + USE KPP_ROOT_Parameters, ONLY: NVAR, LU_NONZERO + USE KPP_ROOT_Global, ONLY: FIX, RCONST, TIME USE KPP_ROOT_Jacobian - USE KPP_ROOT_Jacobian, ONLY : Jac_SP - USE KPP_ROOT_Rates, ONLY : Update_SUN, Update_RCONST - + USE KPP_ROOT_Jacobian, ONLY: Jac_SP + USE KPP_ROOT_Rates, ONLY: Update_SUN, Update_RCONST, Update_PHOTO IMPLICIT NONE KPP_REAL :: T, Told From a28e7e8465046d0108e8a666b6efebc60c78f962 Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Tue, 18 Feb 2025 18:08:30 -0500 Subject: [PATCH 23/28] Only use 100*Roundoff for Backward Euler when computing Hstart int/sdirk.f90 - Indented & lined up code for readability - Added () to subroutines taking no arguments - Test if the method is Backward Euler (BEL) before multiplying the RoundOff variable by 101 when computing Hstart - Trimmed trailing whitespace Signed-off-by: Bob Yantosca --- int/sdirk.f90 | 256 ++++++++++++++++++++++++++------------------------ 1 file changed, 131 insertions(+), 125 deletions(-) diff --git a/int/sdirk.f90 b/int/sdirk.f90 index a43af72..abaa0ed 100644 --- a/int/sdirk.f90 +++ b/int/sdirk.f90 @@ -25,7 +25,7 @@ MODULE KPP_ROOT_Integrator USE KPP_ROOT_LinearAlgebra, ONLY : KppDecomp, KppSolve, Set2zero, & WLAMCH, WCOPY, WAXPY, & WSCAL, WADD - + IMPLICIT NONE PUBLIC SAVE @@ -41,7 +41,7 @@ MODULE KPP_ROOT_Integrator INTEGER, PARAMETER :: Nfun=1, Njac=2, Nstp=3, Nacc=4, & Nrej=5, Ndec=6, Nsol=7, Nsng=8, & Ntexit=1, Nhexit=2, Nhnew=3 - + ! mz_rs_20171120+ ! description of the error numbers IERR CHARACTER(LEN=50), PARAMETER, DIMENSION(-1:1) :: IERR_NAMES = (/ & @@ -80,9 +80,9 @@ SUBROUTINE INTEGRATE( TIN, TOUT, ICNTRL_U, RCNTRL_U, & !~~~> fine-tune the integrator: ICNTRL(2) = 0 ! 0 - vector tolerances, 1 - scalar tolerances - ICNTRL(6) = 0 ! starting values of Newton iterations: + ICNTRL(6) = 0 ! starting values of Newton iterations: ! interpolated (0), zero (1) - ICNTRL(15) = 5 ! Call Update_SUN and Update_RCONST from w/in the int. + ICNTRL(15) = 5 ! Call Update_SUN and Update_RCONST from w/in the int. !~~~> if optional parameters are given, and if they are /= 0, ! then use them to overwrite default settings @@ -102,14 +102,14 @@ SUBROUTINE INTEGRATE( TIN, TOUT, ICNTRL_U, RCNTRL_U, & ! = 2 ! Call Update_PHOTO from within the integrator ! = 3 ! Call Update_RCONST and Update_PHOTO from w/in the int. ! = 4 ! Call Update_SUN from within the integrator - ! = 5 ! Call Update_SUN and Update_RCONST from within the int. + ! = 5 ! Call Update_SUN and Update_RCONST from within the int. ! = 6 ! Call Update_SUN and Update_PHOTO from within the int. ! = 7 ! Call Update_SUN, Update_PHOTO, Update_RCONST w/in int. CALL Integrator_Update_Options( ICNTRL(15), & Do_Update_RCONST, & Do_Update_PHOTO, & Do_Update_Sun ) - + !~~~> In order to remove the prior EQUIVALENCE statements (which !~~~> are not thread-safe), we now have declared VAR and FIX as !~~~> threadprivate pointer variables that can point to C. @@ -133,7 +133,7 @@ SUBROUTINE INTEGRATE( TIN, TOUT, ICNTRL_U, RCNTRL_U, & IF (Ierr < 0) THEN PRINT *,'SDIRK: Unsuccessful exit at T=',TIN,' (Ierr=',Ierr,')' ENDIF - + !~~~> if optional parameters are given for output !~~~> use them to store information in them IF ( PRESENT( ISTATUS_U ) ) ISTATUS_U = ISTATUS @@ -159,17 +159,17 @@ SUBROUTINE SDIRK(N, Tinitial, Tfinal, Y, RelTol, AbsTol, & ! This code is based on the SDIRK4 routine in the above book. ! ! Methods: -! * Sdirk 2a, 2b: L-stable, 2 stages, order 2 -! * Sdirk 3a: L-stable, 3 stages, order 2, adjoint-invariant -! * Sdirk 4a, 4b: L-stable, 5 stages, order 4 -! * Backward Euler: L-stable, 1 stage, order 1 -! Fixed step size = RCNTRL(3) = Hstart +! * Sdirk 2a, 2b: L-stable, 2 stages, order 2 +! * Sdirk 3a: L-stable, 3 stages, order 2, adjoint-invariant +! * Sdirk 4a, 4b: L-stable, 5 stages, order 4 +! * Backward Euler: L-stable, 1 stage, order 1 +! Fixed step size = RCNTRL(3) = Hstart ! ! (C) Adrian Sandu, July 2005 ! Virginia Polytechnic Institute and State University ! Contact: sandu@cs.vt.edu -! Revised by Philipp Miehe and Adrian Sandu, May 2006 -! Backward Euler added by Adrian Sandu, Dec. 2012 +! Revised by Philipp Miehe and Adrian Sandu, May 2006 +! Backward Euler added by Adrian Sandu, Dec. 2012 ! This implementation is part of KPP - the Kinetic PreProcessor !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! @@ -204,7 +204,7 @@ SUBROUTINE SDIRK(N, Tinitial, Tfinal, Y, RelTol, AbsTol, & ! ! Note: For input parameters equal to zero the default values of the ! corresponding variables are used. -!~~~> +!~~~> ! ICNTRL(1) = not used ! ! ICNTRL(2) = 0: AbsTol, RelTol are NVAR-dimensional vectors @@ -282,15 +282,15 @@ SUBROUTINE SDIRK(N, Tinitial, Tfinal, Y, RelTol, AbsTol, & !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IMPLICIT NONE -! Arguments +! Arguments INTEGER, INTENT(IN) :: N, ICNTRL(20) KPP_REAL, INTENT(IN) :: Tinitial, Tfinal, & RelTol(NVAR), AbsTol(NVAR), RCNTRL(20) KPP_REAL, INTENT(INOUT) :: Y(NVAR) INTEGER, INTENT(OUT) :: Ierr - INTEGER, INTENT(INOUT) :: ISTATUS(20) + INTEGER, INTENT(INOUT) :: ISTATUS(20) KPP_REAL, INTENT(OUT) :: RSTATUS(20) - + !~~~> SDIRK method coefficients, up to 5 stages INTEGER, PARAMETER :: Smax = 5 INTEGER, PARAMETER :: S2A=1, S2B=2, S3A=3, S4A=4, S4B=5, BEL=6 @@ -298,7 +298,7 @@ SUBROUTINE SDIRK(N, Tinitial, Tfinal, Y, RelTol, AbsTol, & rkD(Smax), rkE(Smax), rkBhat(Smax), rkELO, & rkAlpha(Smax,Smax), rkTheta(Smax,Smax) INTEGER :: sdMethod, rkS ! The number of stages -! Local variables +! Local variables LOGICAL :: StartNewton KPP_REAL :: Hmin, Hmax, Hstart, Roundoff, & FacMin, Facmax, FacSafe, FacRej, & @@ -306,7 +306,7 @@ SUBROUTINE SDIRK(N, Tinitial, Tfinal, Y, RelTol, AbsTol, & INTEGER :: ITOL, NewtonMaxit, Max_no_steps, i KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0 - + ISTATUS(1:20) = 0 RSTATUS(1:20) = ZERO Ierr = 0 @@ -319,24 +319,24 @@ SUBROUTINE SDIRK(N, Tinitial, Tfinal, Y, RelTol, AbsTol, & ITOL = 0 END IF -!~~~> ICNTRL(3) - method selection +!~~~> ICNTRL(3) - method selection SELECT CASE (ICNTRL(3)) - CASE (0,1) - CALL Sdirk2a() - CASE (2) - CALL Sdirk2b() - CASE (3) - CALL Sdirk3a() - CASE (4) - CALL Sdirk4a() - CASE (5) - CALL Sdirk4b() - CASE (6) - CALL BEuler() - CASE DEFAULT - CALL Sdirk2a() + CASE (0,1) + CALL Sdirk2a() + CASE (2) + CALL Sdirk2b() + CASE (3) + CALL Sdirk3a() + CASE (4) + CALL Sdirk4a() + CASE (5) + CALL Sdirk4b() + CASE (6) + CALL BEuler() + CASE DEFAULT + CALL Sdirk2a() END SELECT - + !~~~> The maximum number of time steps admitted IF (ICNTRL(4) == 0) THEN Max_no_steps = 200000 @@ -363,7 +363,7 @@ SUBROUTINE SDIRK(N, Tinitial, Tfinal, Y, RelTol, AbsTol, & StartNewton = .TRUE. ELSE StartNewton = .FALSE. - END IF + END IF !~~~> Unit roundoff (1+Roundoff>1) Roundoff = WLAMCH('E') @@ -377,7 +377,7 @@ SUBROUTINE SDIRK(N, Tinitial, Tfinal, Y, RelTol, AbsTol, & PRINT * , 'User-selected RCNTRL(1)=', RCNTRL(1) CALL SDIRK_ErrorMsg(-3,Tinitial,ZERO,Ierr) END IF - + !~~~> Upper bound on the step size: (positive value) IF (RCNTRL(2) == ZERO) THEN Hmax = ABS(Tfinal-Tinitial) @@ -387,17 +387,25 @@ SUBROUTINE SDIRK(N, Tinitial, Tfinal, Y, RelTol, AbsTol, & PRINT * , 'User-selected RCNTRL(2)=', RCNTRL(2) CALL SDIRK_ErrorMsg(-3,Tinitial,ZERO,Ierr) END IF - + !~~~> Starting step size: (positive value) IF (RCNTRL(3) == ZERO) THEN - Hstart = MAX(Hmin,101.*Roundoff) ! factor 101 added mz_rs_20171120 + ! If the starting timestep is not supplied, then take the + ! smaller of Hmin and the machine roundoff. Backward Euler + ! needs a starting timestep of at least 100*Roundoff. + ! -- Bob Yantosca (18 Feb 2025) + IF (sdMethod == BEL) THEN + Hstart = MAX(Hmin,101.*Roundoff) ! factor 101 added mz_rs_20171120 + ELSE + Hstart = MAX(Hmin,Roundoff) + ENDIF ELSEIF (RCNTRL(3) > ZERO) THEN Hstart = MIN(ABS(RCNTRL(3)),ABS(Tfinal-Tinitial)) ELSE PRINT * , 'User-selected Hstart: RCNTRL(3)=', RCNTRL(3) CALL SDIRK_ErrorMsg(-3,Tinitial,ZERO,Ierr) END IF - + !~~~> Step size can be changed s.t. FacMin < Hnew/Hexit < FacMax IF (RCNTRL(4) == ZERO) THEN FacMin = 0.2_dp @@ -476,22 +484,22 @@ SUBROUTINE SDIRK(N, Tinitial, Tfinal, Y, RelTol, AbsTol, & END IF END DO END IF - + !~~~> Special treatment of backward Euler IF (sdMethod == BEL) THEN - StartNewton = .FALSE. - FacMin = 1.0_dp - FacMax = 1.0_dp - FacSafe= 1.0_dp - FacRej = 1.0_dp - Qmin = 1.0_dp - Qmax = 1.0_dp - Hmax = Hstart - IF ( Hstart <= 100_dp*Roundoff ) THEN - PRINT * , ' Backward Euler calling error:' - PRINT * , ' it requires Hstart > 100*Roundoff' - STOP - END IF + StartNewton = .FALSE. + FacMin = 1.0_dp + FacMax = 1.0_dp + FacSafe = 1.0_dp + FacRej = 1.0_dp + Qmin = 1.0_dp + Qmax = 1.0_dp + Hmax = Hstart + IF ( Hstart <= 100_dp*Roundoff ) THEN + PRINT* , ' Backward Euler calling error:' + PRINT* , ' it requires Hstart > 100*Roundoff' + STOP + END IF END IF @@ -514,26 +522,26 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Ierr ) USE KPP_ROOT_Parameters IMPLICIT NONE -!~~~> Arguments: +!~~~> Arguments: INTEGER, INTENT(IN) :: N KPP_REAL, INTENT(INOUT) :: Y(NVAR) KPP_REAL, INTENT(IN) :: Tinitial, Tfinal INTEGER, INTENT(OUT) :: Ierr - -!~~~> Local variables: - KPP_REAL :: Z(NVAR,Smax), G(NVAR), TMP(NVAR), & + +!~~~> Local variables: + KPP_REAL :: Z(NVAR,Smax), G(NVAR), TMP(NVAR), & NewtonRate, SCAL(NVAR), RHS(NVAR), & T, H, Theta, Hratio, NewtonPredictedErr, & Qnewton, Err, Fac, Hnew,Tdirection, & NewtonIncrement, NewtonIncrementOld INTEGER :: j, IER, istage, NewtonIter, IP(NVAR) LOGICAL :: Reject, FirstStep, SkipJac, SkipLU, NewtonDone - -#ifdef FULL_ALGEBRA + +#ifdef FULL_ALGEBRA KPP_REAL FJAC(NVAR,NVAR), E(NVAR,NVAR) -#else +#else KPP_REAL FJAC(LU_NONZERO), E(LU_NONZERO) -#endif +#endif KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0 @@ -567,14 +575,14 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Ierr ) IF (IER /= 0) THEN CALL SDIRK_ErrorMsg(-8,T,H,Ierr); RETURN END IF - END IF + END IF IF (ISTATUS(Nstp) > Max_no_steps) THEN CALL SDIRK_ErrorMsg(-6,T,H,Ierr); RETURN - END IF + END IF IF ( (T+0.1d0*H == T) .OR. (ABS(H) <= Roundoff) ) THEN CALL SDIRK_ErrorMsg(-7,T,H,Ierr); RETURN - END IF + END IF stages:DO istage = 1, rkS @@ -584,7 +592,7 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Ierr ) !~~~> Starting values for Newton iterations CALL Set2zero(N,Z(1,istage)) - + !~~~> Prepare the loop-independent part of the right-hand side CALL Set2zero(N,G) IF (istage > 1) THEN @@ -601,7 +609,7 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Ierr ) !~~~> Initializations for Newton iteration NewtonDone = .FALSE. Fac = 0.5d0 ! Step reduction factor if too many iterations - + NewtonLoop:DO NewtonIter = 1, NewtonMaxit !~~~> Prepare the loop-dependent part of the right-hand side @@ -615,17 +623,17 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Ierr ) !~~~> Solve the linear system CALL SDIRK_Solve ( H, N, E, IP, IER, RHS ) - + !~~~> Check convergence of Newton iterations CALL SDIRK_ErrorNorm(N, RHS, SCAL, NewtonIncrement) IF ( NewtonIter == 1 ) THEN Theta = ABS(ThetaMin) - NewtonRate = 2.0d0 + NewtonRate = 2.0d0 ELSE Theta = NewtonIncrement/NewtonIncrementOld IF (Theta < 0.99d0) THEN NewtonRate = Theta/(ONE-Theta) - ! Predict error at the end of Newton process + ! Predict error at the end of Newton process NewtonPredictedErr = NewtonIncrement & *Theta**(NewtonMaxit-NewtonIter)/(ONE-Theta) IF (NewtonPredictedErr >= NewtonTol) THEN @@ -640,14 +648,14 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Ierr ) END IF NewtonIncrementOld = NewtonIncrement ! Update solution: Z(:) <-- Z(:)+RHS(:) - CALL WAXPY(N,ONE,RHS,1,Z(1,istage),1) - + CALL WAXPY(N,ONE,RHS,1,Z(1,istage),1) + ! Check error in Newton iterations NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol) IF (NewtonDone) EXIT NewtonLoop - + END DO NewtonLoop - + IF (.NOT.NewtonDone) THEN !CALL RK_ErrorMsg(-12,T,H,Ierr); H = Fac*H; Reject=.TRUE. @@ -657,7 +665,7 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Ierr ) !~~~> End of simplified Newton iterations - + END DO stages @@ -695,10 +703,10 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Ierr ) !~~~> Update time and solution T = T + H ! Y(:) <-- Y(:) + Sum_j rkD(j)*Z_j(:) - DO i = 1,rkS + DO i = 1,rkS IF (rkD(i)/=ZERO) CALL WAXPY(N,rkD(i),Z(1,i),1,Y,1) - END DO - + END DO + !~~~> Update scaling coefficients CALL SDIRK_ErrorScale(N, ITOL, AbsTol, RelTol, Y, SCAL) @@ -723,7 +731,7 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Ierr ) ! If convergence is fast enough, do not update Jacobian ! SkipJac = (Theta <= ThetaMin) SkipJac = .FALSE. - + ELSE accept !~~~> Step is rejected IF (FirstStep .OR. Reject) THEN @@ -733,16 +741,16 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Ierr ) END IF Reject = .TRUE. SkipJac = .TRUE. - SkipLU = .FALSE. + SkipLU = .FALSE. IF (ISTATUS(Nacc) >= 1) ISTATUS(Nrej) = ISTATUS(Nrej) + 1 - + END IF accept - + END DO Tloop ! Successful return Ierr = 1 - + END SUBROUTINE SDIRK_Integrator @@ -763,14 +771,14 @@ SUBROUTINE SDIRK_ErrorScale(N, ITOL, AbsTol, RelTol, Y, SCAL) END DO END IF END SUBROUTINE SDIRK_ErrorScale - - + + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE SDIRK_ErrorNorm(N, Y, SCAL, Err) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -! +! INTEGER :: i, N - KPP_REAL :: Y(N), SCAL(N), Err + KPP_REAL :: Y(N), SCAL(N), Err Err = ZERO DO i=1,N Err = Err+(Y(i)*SCAL(i))**2 @@ -819,17 +827,17 @@ SUBROUTINE SDIRK_ErrorMsg(Code,T,H,Ierr) PRINT *, "T=", T, "and H=", H END SUBROUTINE SDIRK_ErrorMsg - - + + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE SDIRK_PrepareMatrix ( H, T, Y, FJAC, & SkipJac, SkipLU, E, IP, Reject, ISING ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> Compute the matrix E = 1/(H*GAMMA)*Jac, and its decomposition !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - + IMPLICIT NONE - + KPP_REAL, INTENT(INOUT) :: H KPP_REAL, INTENT(IN) :: T, Y(NVAR) LOGICAL, INTENT(INOUT) :: SkipJac,SkipLU,Reject @@ -846,9 +854,9 @@ SUBROUTINE SDIRK_PrepareMatrix ( H, T, Y, FJAC, & ConsecutiveSng = 0 ISING = 1 - + Hloop: DO WHILE (ISING /= 0) - + HGammaInv = ONE/(H*rkGamma) !~~~> Compute the Jacobian @@ -858,8 +866,8 @@ SUBROUTINE SDIRK_PrepareMatrix ( H, T, Y, FJAC, & IF (.NOT. SkipJac) THEN CALL JAC_CHEM( T, Y, FJAC ) ISTATUS(Njac) = ISTATUS(Njac) + 1 - END IF - + END IF + #ifdef FULL_ALGEBRA DO j=1,NVAR DO i=1,NVAR @@ -889,7 +897,7 @@ SUBROUTINE SDIRK_PrepareMatrix ( H, T, Y, FJAC, & SkipLU = .FALSE. Reject = .TRUE. END IF - + END DO Hloop END SUBROUTINE SDIRK_PrepareMatrix @@ -911,19 +919,17 @@ SUBROUTINE SDIRK_Solve ( H, N, E, IP, ISING, RHS ) #endif KPP_REAL, INTENT(INOUT) :: RHS(N) KPP_REAL :: HGammaInv - + HGammaInv = ONE/(H*rkGamma) CALL WSCAL(N,HGammaInv,RHS,1) -#ifdef FULL_ALGEBRA +#ifdef FULL_ALGEBRA CALL DGETRS( 'N', N, 1, E, N, IP, RHS, N, ISING ) #else CALL KppSolve(E, RHS) #endif ISTATUS(Nsol) = ISTATUS(Nsol) + 1 - - END SUBROUTINE SDIRK_Solve - + END SUBROUTINE SDIRK_Solve !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE BEuler() @@ -934,7 +940,7 @@ SUBROUTINE BEuler() rkS = 1 ! Method coefficients - rkGamma = 1.0_dp + rkGamma = 1.0_dp rkA(1,1) = 1.0_dp rkB(1) = 1.0_dp rkBhat(1)= 1.0_dp @@ -1041,7 +1047,7 @@ SUBROUTINE Sdirk4a() rkAlpha(5,2) = 6.559571569643355712998131800797873d0 rkAlpha(5,3) = -15.90772144271326504260996815012482d0 rkAlpha(5,4) = 25.34908987169226073668861694892683d0 - + !~~~> Coefficients for continuous solution ! rkD(1,1)= 24.74416644927758d0 ! rkD(1,2)= -4.325375951824688d0 @@ -1070,9 +1076,9 @@ SUBROUTINE Sdirk4a() ! CALL WAXPY(N,rkD(i,j),Z(1,j),1,CONT(1,i),1) ! END DO ! END DO - + rkELO = 4.0d0 - + END SUBROUTINE Sdirk4a !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -1160,7 +1166,7 @@ SUBROUTINE Sdirk4b() rkAlpha(5,2) = 9.d0 rkAlpha(5,3) = -56.81818181818181818181818181818182d0 rkAlpha(5,4) = 54.d0 - + END SUBROUTINE Sdirk4b !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -1203,7 +1209,7 @@ SUBROUTINE Sdirk2a() ! Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j} rkAlpha(2,1) = 3.414213562373095048801688724209698d0 - + END SUBROUTINE Sdirk2a !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -1246,7 +1252,7 @@ SUBROUTINE Sdirk2b() ! Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j} rkAlpha(2,1) = .5857864376269049511983112757903019d0 - + END SUBROUTINE Sdirk2b @@ -1314,24 +1320,25 @@ END SUBROUTINE SDIRK ! AND ALL ITS INTERNAL PROCEDURES SUBROUTINE FUN_CHEM( T, Y, P ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - USE KPP_ROOT_Parameters, ONLY: NVAR, LU_NONZERO - USE KPP_ROOT_Global, ONLY: FIX, RCONST, TIME - USE KPP_ROOT_Function, ONLY: Fun - USE KPP_ROOT_Rates, ONLY: Update_SUN, Update_RCONST, Update_PHOTO + USE KPP_ROOT_Parameters, ONLY : NVAR, LU_NONZERO + USE KPP_ROOT_Global, ONLY : FIX, RCONST, TIME + USE KPP_ROOT_Function, ONLY : Fun + USE KPP_ROOT_Rates, ONLY : Update_SUN, Update_RCONST + IMPLICIT NONE KPP_REAL :: T, Told KPP_REAL :: Y(NVAR), P(NVAR) - + Told = TIME TIME = T IF ( Do_Update_SUN ) CALL Update_SUN() IF ( Do_Update_RCONST ) CALL Update_RCONST() - + CALL Fun( Y, FIX, RCONST, P ) - + TIME = Told - + END SUBROUTINE FUN_CHEM @@ -1339,13 +1346,14 @@ END SUBROUTINE FUN_CHEM SUBROUTINE JAC_CHEM( T, Y, JV ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - USE KPP_ROOT_Parameters, ONLY: NVAR, LU_NONZERO - USE KPP_ROOT_Global, ONLY: FIX, RCONST, TIME + USE KPP_ROOT_Parameters, ONLY : NVAR, LU_NONZERO + USE KPP_ROOT_Global, ONLY : FIX, RCONST, TIME USE KPP_ROOT_Jacobian - USE KPP_ROOT_Jacobian, ONLY: Jac_SP - USE KPP_ROOT_Rates, ONLY: Update_SUN, Update_RCONST, Update_PHOTO + USE KPP_ROOT_Jacobian, ONLY : Jac_SP + USE KPP_ROOT_Rates, ONLY : Update_SUN, Update_RCONST + IMPLICIT NONE - + KPP_REAL :: T, Told KPP_REAL :: Y(NVAR) #ifdef FULL_ALGEBRA @@ -1354,7 +1362,7 @@ SUBROUTINE JAC_CHEM( T, Y, JV ) #else KPP_REAL :: JV(LU_NONZERO) #endif - + Told = TIME TIME = T IF ( Do_Update_SUN ) CALL Update_SUN() @@ -1378,5 +1386,3 @@ SUBROUTINE JAC_CHEM( T, Y, JV ) END SUBROUTINE JAC_CHEM END MODULE KPP_ROOT_Integrator - - From 2ca322da6cb48b902d3a54cbc07d125010a3ee0e Mon Sep 17 00:00:00 2001 From: Rolf Sander Date: Wed, 19 Feb 2025 08:40:20 +0100 Subject: [PATCH 24/28] comment corrected --- int/sdirk.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/int/sdirk.f90 b/int/sdirk.f90 index abaa0ed..490586f 100644 --- a/int/sdirk.f90 +++ b/int/sdirk.f90 @@ -391,7 +391,7 @@ SUBROUTINE SDIRK(N, Tinitial, Tfinal, Y, RelTol, AbsTol, & !~~~> Starting step size: (positive value) IF (RCNTRL(3) == ZERO) THEN ! If the starting timestep is not supplied, then take the - ! smaller of Hmin and the machine roundoff. Backward Euler + ! larger of Hmin and the machine roundoff. Backward Euler ! needs a starting timestep of at least 100*Roundoff. ! -- Bob Yantosca (18 Feb 2025) IF (sdMethod == BEL) THEN From 24f5efa4fbe5425b3d0a89a2e09b0b1a3bf516c1 Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Wed, 19 Feb 2025 10:29:50 -0500 Subject: [PATCH 25/28] Updated "KPP revision history" with info about updated sdirk docs/source/getting_started/00_revision_history.rst - Added notes about beuler integrator being removed, and that users can use sdirk with ICNTRL(3)=6 to request Backward Euler Signed-off-by: Bob Yantosca --- docs/source/getting_started/00_revision_history.rst | 3 +++ 1 file changed, 3 insertions(+) diff --git a/docs/source/getting_started/00_revision_history.rst b/docs/source/getting_started/00_revision_history.rst index f4c342e..3465aaa 100644 --- a/docs/source/getting_started/00_revision_history.rst +++ b/docs/source/getting_started/00_revision_history.rst @@ -20,6 +20,9 @@ KPP 3.2.0 - Updated code in :code:`src/gen.c` to generate the :code:`UPDATE_RCONST` routine with an optional argument :code:`Y` - Updated C-I tests to print the compiler versions that are used +- Updated :literal:`int/sdirk.f90` to a newer version +- Removed :literal:`int/beuler.f90`; Users can select Backward Euler + with :literal:`sdirk` integrator and :literal:`ICNTRL(3)=6` .. _kpp311: From 1a3949536fd2ee6d688bc4068b4826088310cae1 Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Wed, 19 Feb 2025 12:19:15 -0500 Subject: [PATCH 26/28] Updated "KPP revision history" for PR #127 docs/source/getting_started/00_revision_history.rst - Added a line to KPP 3.2.0 about MacOS architecture specific build compiler flags Signed-off-by: Bob Yantosca --- docs/source/getting_started/00_revision_history.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/source/getting_started/00_revision_history.rst b/docs/source/getting_started/00_revision_history.rst index 3465aaa..978b758 100644 --- a/docs/source/getting_started/00_revision_history.rst +++ b/docs/source/getting_started/00_revision_history.rst @@ -23,6 +23,7 @@ KPP 3.2.0 - Updated :literal:`int/sdirk.f90` to a newer version - Removed :literal:`int/beuler.f90`; Users can select Backward Euler with :literal:`sdirk` integrator and :literal:`ICNTRL(3)=6` +- Added MacOS architecture-specific compilation flags to the build sequence .. _kpp311: From a97b8e3bea716854bdd8e6bb8b1c9756b548b4c9 Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Thu, 27 Feb 2025 11:18:15 -0500 Subject: [PATCH 27/28] Update version numbers to 3.2.0 in the usual places Updated version numbers in - CHANGELOG.md - docs/source/conf.py (for RTD documentation) - src/gdata.h Signed-off-by: Bob Yantosca --- CHANGELOG.md | 2 +- docs/source/conf.py | 2 +- src/gdata.h | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 6e57ca6..8839071 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,7 +8,7 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). -## [Unreleased] - TBD +## [3.2.0] - 2025-02-27 ### Added - Added new inline key `F90_RCONST_USE` in `src/gdata.h` and `src/scanner.c` - Added documentation about `F90_RCONST_USE` for ReadTheDocs diff --git a/docs/source/conf.py b/docs/source/conf.py index 6272ce4..fdac8ec 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -25,7 +25,7 @@ # The full version, including alpha/beta/rc tags # (version numbers must be synchronized in CHANGELOG.md, src/gdata.h, # and docs/source/conf.py) -release = "3.1.1" +release = "3.2.0" # -- General configuration --------------------------------------------------- diff --git a/src/gdata.h b/src/gdata.h index 435337b..7cd638e 100644 --- a/src/gdata.h +++ b/src/gdata.h @@ -31,7 +31,7 @@ // Version numbers must be synchronized in CHANGELOG.md, src/gdata.h, // and docs/source/conf.py -#define KPP_VERSION "3.1.1" +#define KPP_VERSION "3.2.0" #ifndef _GDATA_H_ #define _GDATA_H_ From 00f18781683cd5d439716f4642addebc38b99430 Mon Sep 17 00:00:00 2001 From: Rolf Sander Date: Thu, 27 Feb 2025 19:35:09 +0100 Subject: [PATCH 28/28] minor comment changes --- .ci-pipelines/ci-common-defs.sh | 2 +- src/Makefile.defs | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/.ci-pipelines/ci-common-defs.sh b/.ci-pipelines/ci-common-defs.sh index 62340d3..ce35bf8 100755 --- a/.ci-pipelines/ci-common-defs.sh +++ b/.ci-pipelines/ci-common-defs.sh @@ -46,7 +46,7 @@ function get_ci_test_path() { return } -# Prints a headeer with the compiler versions +# Prints a header with the compiler versions function print_compiler_versions() { echo \ "###########################################################################" diff --git a/src/Makefile.defs b/src/Makefile.defs index 2199256..a6c9244 100644 --- a/src/Makefile.defs +++ b/src/Makefile.defs @@ -73,9 +73,9 @@ FLEX := flex FLEX_LIB_DIR := /usr/lib BISON := bison INCLUDE_DIR := /usr/include -SYSTEM := $(shell uname) # e.g. "Linux", "Darwin" -SYSTEM_M := $(shell uname -m) # e.g. "x86_64", "arm64" -HOST := $(shell hostname) # name of machine +SYSTEM := $(shell uname) # kernel name, e.g., "Linux", "Darwin" +SYSTEM_M := $(shell uname -m) # machine hardware name, e.g., "x86_64", "arm64" +HOST := $(shell hostname) # name of machine ############################################################ ### Keep looking for the flex library until found ###