Cleanup of integrator files

[RolfSander]

One of the features I like most about KPP is that you can choose the
most suitable integrator for your problem and even fine-tune it.
However, for me as a non-mathematician, it is also important that the
integrator works out-of-the-box. Therefore I suggest that we check all
int/*.f90 files and either decide to maintain them, or move them into
the int/user_contributed/ directory. For those that we decide to
maintain, I suggest that we look at a couple of issues:

• A consistent naming scheme would be good. For example, why do some
  files have the prefix kpp_? Can we remove that prefix?

• I'd like to check if [IR]CNTRL and [IR]STATUS are used consistently
  between different integrators, as far as possible.

• I'd like to check if the error numbers and messages are used
  consistently between different integrators, as far as possible.

• Do we need the empty file int/none.f90 ?

• We have several rosenbrock variants as individual files. Would it be
  possible to merge them into one file? Their different features could
  be switched on/off via ICNTRL. Well, this is probably more of a
  long-term project...

@yantosca, @msl3v: Any comments?

[yantosca]

Thanks @RolfSander, I think it is a good idea to check for consistency between the integrator options. Just the other day I ran into an issue because the GEOS-Chem interface was expecting the Nhnew parameter to be found in the integrator file (which is true if you are using rosenbrock but not if you are using feuler). So that led to a compile time error. There may be other things like that to sort out. If you can take the point on it, that would be great.

You may be right, we might not need int/none.f90. As long as #DRIVER none isn't trying to inline that we can probably skip it.

I think some of the integrators that use kpp_ in the file name are just historical baggage, and we can probably remove that.

It'd be great to merge the rosenbrock options. I'd have to take a look to see if they are the same file but with just different ICNTRL settings. If so we can merge that in.

[RolfSander]

I think it is a good idea to check for consistency between the integrator options.
If you can take the point on it, that would be great.

Yes, I'll work on this.

You may be right, we might not need int/none.f90. As long as #DRIVER
none isn't trying to inline that we can probably skip it.

I think that's okay because in gen.c we have:

if( strcmp( driver, "none" ) != 0 )
    GenerateDriver();

I think some of the integrators that use kpp_ in the file name are just
historical baggage, and we can probably remove that.

I have renamed them.

It'd be great to merge the rosenbrock options. I'd have to take a look
to see if they are the same file but with just different ICNTRL
settings. If so we can merge that in.

Here's a quick overview:

• rosenbrock_adj.f90 and rosenbrock_tlm.f90 are the adjoint and the
  tangent linear model, respectively. I guess these should remain as
  separate files.

• rosenbrock_mz.f90 is based on an older version (from KPP-2.1) and
  specific to our global model. We may want to move it into
  user_contributed.

• rosenbrock.f90 is probably the original file on which the following
  files are based.

• rosenbrock_posdef.f90 discards negative values which can avoid
  numerical problems. We could activate this as a new KPP command
  "#POSDEF".

• rosenbrock_posdef_h211b_qssa.f90 has an adaptive time-stepping
  according to Söderlind (2002, doi:10.1145/641876.641877)

• rosenbrock_split.f90 separates dc/dt into production and
  destruction. Maybe we could activate this feature with a new KPP
  command "#SPLIT". @msl3v what do you think?

Further differences between the files are:

• old BLAS functions replaced by simple Fortran90 commands

• cosmetic changes (white space, separating lines, etc.)

• renaming of variables (NVAR->N)

Merging the Rosenbrock files is certainly worth the effort, but I think
there's quite a bit of work ahead for us...

[RolfSander]

I removed the trivial differences (whitespace etc.) between
rosenbrock.f90, rosenbrock_posdef.f90 and rosenbrock_split.f90.
Now it's easier to see what the real differences are.

I think it should be straightforward to merge rosenbrock_posdef.f90
into rosenbrock.f90 with a new KPP command "#POSDEF". Or, even better,
use ICNTRL to switch between the options.

Regarding rosenbrock_split.f90, I'm not so sure. Maybe we can return
Aout, P_VARout and D_VARout as OPTIONAL parameters from
SUBROUTINE Fun_SPLIT in the same way as we return Vdotout from
SUBROUTINE Fun? Of course, only if useAggregate=0.
@msl3v and @yantosca, what do you think?

[msl3v]

@RolfSander , are you not sure about merging rosenblock_split.f90 into rosenbrock.f90? It doesn't seem straightforward. Computationally, I think it's cleaner to keep it separate, since it is driven by different preprocessor output, and I've been religiously avoiding branching statements (if else).

[RolfSander]

merging rosenblock_split.f90 into rosenbrock.f90?

I think it's cleaner to keep it separate, since it is driven by
different preprocessor output, and I've been religiously avoiding
branching statements (if else).

I see your point @msl3v. #FUNCTION AGGREGATE is evaluated already when KPP is
run. Do you think that such a branching would slow down the code
cosiderably?

IF (split) THEN
  CALL Fun_Split
ELSE
  CALL Fun
ENDIF

An alternative could be that we insert a placeholder like
CALL Fun_OR_Fun_Split into rosenbrock.f90 and let KPP modify it to
either CALL Fun or CALL Fun_Split. This would be similar to the way
that KPP replaces KPP_REAL by REAL(kind=sp) or REAL(kind=dp).

The reason why I prefer a unified file is that I fear that
rosenblock_split.f90 and rosenbrock.f90 would diverge soon if we
keep them separate.

[msl3v]

Re. branching: I don't know how much of an impact it would have. Though my understanding is that, using

IF (split) THEN
CALL Fun_Split
ELSE
CALL Fun
ENDIF

the compiler would try to preload Fun_Split in every iteration, which would definitely slow things down if (.not. split).

Having KPP overwrite something like _KPPFUNCTION is a great idea! Indeed this approach in general would allow really fine control of the code structure w/o introduction of CPP tags or branching statements.

[RolfSander]

OK, then let's use this approach for Fun_Split and Fun!

this approach in general would allow really fine control of the code
structure w/o introduction of CPP tags or branching statements.

Yes, we can use this whenever we need to modify something in the
integrator file, as long as it is something that will not change during
runtime of the model. Then we have to use IF blocks and ICNTRL.

[yantosca]

@RolfSander: I have run into an issue with the radau5.f90 integrator. The function k_3rd_jpl_activation causes the radau90 C-I test to fail

radau90_Rates.f90:133:30:`
  133 |     k_3rd_jpl_activation(ASSOC)  = k_f
      |                              1
Error: Symbol 'assoc' at (1) has no IMPLICIT type
radau90_Rates.f90:134:31:

  134 |     k_3rd_jpl_activation(DISSOC) = k_fCA
      |                               1
Error: Symbol 'dissoc' at (1) has no IMPLICIT type

Where are the ASSOC and DISSOC parameters defined?

[RolfSander]

Sorry, somehow the CI-tests don't start when I push into this branch.

I hope the fix in gen.c works.

[yantosca]

The C-I tests now work again. Had to change #INTEGRATOR kpp_radau5 to #INTEGRATOR radau5 in the radau90.kpp file. All good!

[RolfSander]

I have performed several tests with different integrators and a simple
ozone+NOx chemistry. Here is a summary of the results:

• We already have the integrators LSODE, RADAU5, ROSENBROCK,
  RUNGE_KUTTA, and SDIRK in our CI-tests. It's no surprise that they
  also passed all of my tests.

• The integrators ROSENBROCK_POSDEF and SDIRK4 are not included in our
  CI-tests but they passed as well.

• The integrators EXPONENTIAL and ROSENBROCK_SPLIT complain that
  "fun_split is already defined". It seems that there is a duplicate
  generation of the SPLIT routines. This needs to be fixed.

• The integrators FEULER and BEULER run into numerical problems. I guess
  that even my simple ozone and NOx chemistry was already too stiff for
  these solvers. Maybe we can create a chemical system for which these
  solvers do work?

• The integrators DVODE and SEULEX produce compiler errors. My
  suggestion: Repair or remove.

• The integrators GILLESPIE and TAU_LEAP have the #STOCHASTIC option
  switched on. The compilation crashed because the Makefile was not
  compatible (apparently you need the command "make stochastic" instead
  of "make"). My suggestion: Repair or remove.

• I did not check the _adj and _tlm variants.

The decision "repair or remove" may be tough but I prefer not to have
any dysfunctional integrators in our official release 3.0.0.

@yantosca, @msl3v: What do you think?

[RolfSander]

More questions for you, @yantosca and @msl3v. For 4 of our integrators,
we only have FORTRAN77 versions:

atm_lsodes.f
atm_odessa_ddm.f
atm_radau5.f
odessa_ddm.f

I suggest to delete them. If anyone is interested, we can always recover
the old code from KPP-2.5.0 and convert it to Fortran 90 later.

[yantosca]

I have no problem with this. We typically use Rosenbrock and Forward Euler for our simulations.

I guess this is a wider question -- should we also go thru the KPP code and remove code_f77.c and all other references to F77 output?

[RolfSander]

Removing all references to F77 might be a lot of work. I think it would
be sufficient to state in our documentation that F77 is deprecated and
not maintained anymore.

[yantosca]

@RolfSander that's fine w/ me. I can add that to the docs.

[RolfSander]

OK, thanks!

What is your opinion on DVODE, SEULEX, GILLESPIE, and TAU_LEAP (see above)?

[yantosca]

Hi @RolfSander, let me scope out the issues with DVODE & SEULEX. If it's a quick fix I'll push an update. Otherwise we'll remove them.

[RolfSander]

I guess as long as there are F90 compilers, I'd leave them. We do
reference Gillespie in the docs.

Yes, I'd also like to keep GILLESPIE and TAU_LEAP. If the Makefile is
the only issue with them, we could simply mention in the docs that you
compile with make stochastic instead of make.

let me scope out the issues with DVODE & SEULEX. If it's a quick fix
I'll push an update. Otherwise we'll remove them.

Sounds good. Thanks for checking!

[yantosca]

SEULEX Is OK, it just needed a 'USE KPP_Root_Global` statement at the top of the module.

I've managed to fix all the warnings and errors except for:

small_strato_Integrator.f90:361:78:

  361 | F90(FUN_CHEM,NVAR,VAR,T1,T2,ITASK,ISTATE,OPTIONS,J_FCN=JAC_CHEM)
      |                                                                1

Error: There is no specific subroutine for the generic ‘dvode_f90’ at (1
make: *** [Makefile_small_strato:249: small_strato_Integrator.o] Error 1

I think it has to do with the manual INTERFACE declaration for J_FCN, and that not lining up w/ the actual arguments of FUN_CHEM. I need to move on to something else for now but I'll try to fix that later.

[RolfSander]

Thanks, it's great that you were able to repair SEULEX!

Hmm, the INTERFACE DVODE_F90 with MODULE PROCEDURE statements looks
pretty complex.

Next week, I'll be attending the EGU conference (virtually) and after
that I'll take a few days off. I'm not sure if I can find much time to
work on KPP until the beginning of June...

[yantosca]

No worries @RolfSander. I will be going to IGC10 in St. Louis in early June so will have to prep for that. At least we got KPP 2.5.0 out which will be something for us to trumptet at IGC10.

[RolfSander]

I'm using ICNTRL(16)=1 now in rosenbrock.f90 to set negative
concentrations to zero. The separate files rosenbrock_posdef.* are not necessary anymore.

[RolfSander]

I'm working on FUN and FUN_SPLIT now.

@yantosca : In commit b95b4ea, you added this code to
GenerateFun() in gen.c:

 for (i = VarNr; i < VarNr; i++) {
      sum = Const(0);
      for (j = 0; j < EqnNr; j++)
        sum = Add( sum, Mul( Const( Stoich[i][j] ), Elm( A, j ) ) );
      Assign( Elm( Vdot, i ), sum );
    }

I'm trying to understand what it means. The loop starts with i = VarNr
and ends at i < VarNr. What is the purpose of this for-loop?

[RolfSander]

There is no separate rosenblock_split integrator anymore. Instead, you
now just use the command #FUNCTION SPLIT and use the normal
rosenbrock.f90 integrator. It should work fine but please let me know
if you discover any problems. I also adapted the exponential.f90
solver but that hasn't been tested yet. There is also a new CI test
(ros_split) but for unknown reasons the CI tests don't execute when I
push to the cleanup_int branch, and ci-manual-testing-script.sh
doesn't work for me either.

[msl3v]

I have to applaud the work you guys are doing. I wish I was more available to help with this directly. If you need more specific things, please let me know.

…

On 5/22/22 07:01, Rolf Sander wrote: There is no separate |rosenblock_split| integrator anymore. Instead, you now just use the command |#FUNCTION SPLIT| and use the normal |rosenbrock.f90| integrator. It should work fine but please let me know if you discover any problems. I also adapted the |exponential.f90| solver but that hasn't been tested yet. There is also a new CI test (|ros_split|) but for unknown reasons the CI tests don't execute when I push to the cleanup_int branch, and |ci-manual-testing-script.sh| doesn't work for me either. — Reply to this email directly, view it on GitHub <#40 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABVE2334CK5N7BBTRCH5RKDVLIHZJANCNFSM5V5ZVMXQ>. You are receiving this because you were mentioned.Message ID: ***@***.***>

[msl3v]

@RolfSander @yantosca I think the AR-boxmodel can be ingested into the $KPP_HOME/models dir. Would you guys like me to do this?

[RolfSander]

Thanks, @msl3v!

Maybe you could just try to run Rosenbrock with the SPLIT option on your computer and let me know if it works? If not, I can take care of fixing the problem...

[RolfSander]

I don't know the AR-boxmodel. If you think you can put it into a *.def, a *.spc and an *.eqn file in models/, it would certainly be a nice addition to KPP!

[RolfSander]

I'm working on ICNTRL consistency amongst the solvers now. As a start, here is a summary table for the f90 solvers:

ICNTRL	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16
beuler	-	+	+	+	+	+	-	-	-	-	-	-	-	-	+	-
dvode	-	-	-	-	-	-	-	-	-	-	-	-	-	-	+	-
exponential	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
feuler	?	?	-	-	-	-	-	-	-	-	-	-	-	-	+	-
gillespie	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
lsode	-	+	-	+	?	-	-	-	-	-	-	-	-	-	+	-
radau5	-	+	-	+	+	+	-	-	-	-	+	-	-	-	+	-
rosenbrock_adj	+	+	+	+	-	?	?	?	-	-	-	-	-	-	+	-
rosenbrock	+	+	+	+	-	-	-	-	-	-	-	-	-	-	+	+
rosenbrock_tlm	+	+	+	+	-	-	-	-	-	-	-	?	-	-	+	-
runge_kutta_adj	-	+	+	+	+	?	?	?	?	?	+	-	-	-	+	-
runge_kutta	-	+	+	+	+	+	-	-	-	?	+	-	-	-	+	-
runge_kutta_tlm	-	+	+	-	+	+	?	-	?	?	+	?	-	-	+	-
sdirk4	-	+	-	+	-	-	-	-	-	-	-	-	-	-	+	-
sdirk_adj	-	+	+	+	+	+	?	?	-	-	-	-	-	-	+	-
sdirk	-	+	+	+	+	+	-	-	-	-	-	-	-	-	+	-
sdirk_tlm	-	+	+	+	+	+	?	-	?	-	-	?	-	-	+	-
seulex	+	+	-	+	-	-	-	-	-	?	?	?	?	?	+	-
tau_leap	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-

- = not used
+ = used
? = inconsistent (solver-specific) usage

[RolfSander]

Here is a list of tasks related to ICNTRL:

1. In our documentation, ICNTRL is currently in the "Information for KPP
   developers" section (07_numerical_methods). I don't think you
   have to be a developer to use ICNTRL. I suggest to move the
   description to the "Input for KPP" section (04_input_for_kpp).
   Likewise, the description of I/RSTATUS could go to "Output from KPP"
   (05_output_from_kpp).

2. The solver feuler.f90 uses ICNTRL(1) for verbose error output and
   ICNTRL(2) to stop upon negative results. This is incompatible with the
   usage in other solvers. I suggest:
   a) Use ICNTRL(16) to deal with negative results in a similar way as
   rosenbrock.f90 does it now, e.g.:
   0 -> don't do anything
   1 -> set negative values to zero
   2 -> stop.
   b) Use ICNTRL(17) to control error output. This is a nice feature that
   we could also use in other integrators. Let's use a new index in the
   ICNTRL array for this.

3. Using ICNTRL(5) for maximal order in lsode.f90 collides with many
   other solvers which use it for the maximum number of Newton
   iterations. We could use ICNTRL(10) instead, which is already
   compiler-specific.

4. There is a bug in the usage of ICNTRL(11) in the solvers radau5.f90,
   runge_kutta_adj.f90, and runge_kutta_tlm.f90. It seems to me that the
   code in runge_kutta.f90 is correct and could be adopted.

5. The solvers gillespie.f90 and tau_leap.f90 don't even read ICNTRL.
   Maybe this is intentional, I don't know.

If you agree, @yantosca and @msl3v, I could work on the first 3 items. Regarding 4 and 5, I'm not sure what to do. Maybe ask Adrian.

[yantosca]

Thanks @RolfSander. Please go ahead with the #1, #2, and #3. I think those are all good updates.

I remember looking at the tau_leap and gillespie integrators. Not sure if many of the ICNTRL options apply to them (maybe only ICNTRL(15) would, or ICNTRL(17) if we add it in). We could just add those features in as time allows.

[msl3v]

OK. This is good guidance. I will make the needed changes to feuler.f90

[RolfSander]

Thanks, @msl3v, for updating feuler.f90!

I will take care of 1. and 3.

I don't know anyone who is actually using tau_leap and gillespie. For
now, I suggest that we simply mention in the docs that they don't make
use of ICNTRL.

[yantosca]

Hi @RolfSander and @msl3v. I was away for a few days, I won't be able to test the rosenbrock_split code in GEOS-Chem until I get back from IGC10 in mid-June.

About the weird for loop, let's comment it out and then see if it doesn't break anything.

[RolfSander]

No problem. Enjoy the GEOS-Chem meeting and let's continue when you're back...

[yantosca]

We can close out this issue as PR #44 and PR #47 have been merged.

I will open a new issue where we can discuss further cleanup of the rate law functions.

[msl3v]

Modified this: 1) IF( ICNTRL(16) .eq. 0 ) Don't even perform the search for negatives. 2) Three ICNTRL(16) options instead of two     0 -> Don't do anything     1 -> Set negative values to zero     2 -> Return     3 -> Stop Pushed to forwardeuler branch. Pull requesting in a moment. Bob this will require changing the ICNTRL settings in the carbon cycle code. ccyclechem_mod.F90. You could probably do this easily in GCv14.0.0, or I could do it in the CCycle branch and create a pull request. Your call.

…

On 5/23/22 10:41, Rolf Sander wrote: The solver feuler.f90 uses ICNTRL(1) for verbose error output and ICNTRL(2) to stop upon negative results. This is incompatible with the usage in other solvers. I suggest: a) Use ICNTRL(16) to deal with negative results in a similar way as rosenbrock.f90 does it now, e.g.: 0 -> don't do anything 1 -> set negative values to zero 2 -> stop. b) Use ICNTRL(17) to control error output. This is a nice feature that we could also use in other integrators. Let's use a new index in the ICNTRL array for this.