Create variables to track non-passive species, so that they can be excluded from the Rosenbrock error norm (in int/rosenbrock*.f90 integrators only)

[yantosca]

Overview

This is the companion PR to #66 and #144. In this PR we allow optional arguments to be passed to Rosenbrock F90 integrators so that we can exclude dummy species (i.e. those added to chemical reactions for the purpose of obtaining diagnostic information) from the Rosenbrock error norm computation. This should prevent the issue where integration results are changed depending on the number of dummy species in the mechanism (cf geoschem/geos-chem#3175).

Technical description

Two optional arguments are passed to the top-level INTEGRATE functions in the int/rosenbrock.f90, int/rosenbrock_autoreduce.f90, int/rosenbrock_adj.f90 and int/rosenbrock_tlm.f90 integrators. These are:

1. NonDummySpc_Count: The number of non-dummy species
2. NonDummySpc_Indices: The KPP species indices for non-dummy species

These optional variables can be defined outside of KPP by using a simple test of absolute tolerance, such as:

! For example, in GEOS-Chem
INTEGER :: NonDummySpc_Count
INTEGER :: NonDummySpc_Indices(NVAR)
. . .
NonDummySpcCount = 0
DO I = 1, NVAR
   IF ( ATOL(I) < 1.0d25 ) THEN
      NonDummySpc_Count = NonDummySpc_Count + 1
      NonDummySpc_Indices(NonDummySpc_Count) = I
   ENDIF
ENDDO

This only has to be done once, e.g. at model initialization. Then these can be passed to the Rosenbrock integrator:

       CALL Integrate( TIN, TOUT, ICNTRL, RCNTRL, ISTATUS, RSTATE, IERR, &
                       NonDummySpc_Count, NonDummySpc_Indices )

We have also updated the algorithm in the Ros_ErrorNorm function (in each of the int/rosenbrock*.f90 files) to compute the error norm on the basis of non-dummy species when the optional arguments are passed. If they are not passed, it will fall back to the default behavior (using all species to compute the error norm). Also we have optimized the function to facilitate vectorization and to remove unnecessary ELSE blocks, which can be a bottleneck:

   Err = ZERO

   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   ! If optional args NonDummySpc_Count and NonDummySpc_Indices are passed,
   ! then use them to exclude dummy species from the error norm computation.
   !
   ! Use separate loops for scalar and vector tolerances, as this
   ! will facilitate loop vectorization for better performance.
   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   IF ( PRESENT( NonDummySpc_Count   )  .and. &
        PRESENT( NonDummySpc_Indices ) ) THEN

      !~~~> Vector Tolerance (per-species AbsTol & RelTol)
      IF ( VectorTol ) THEN

         DO IDX = 1, NonDummySpc_Count
            I     = NonDummySpc_Indices(IDX)
            Ymax  = MAX( ABS( Y(I) ), ABS( Ynew(I) ) )
            Scale = AbsTol(I) + RelTol(I) * Ymax
            Err   = Err + ( Yerr(I) / Scale )**2
         ENDDO

         ! Normalize the error norm by the number of non-dummy species
         ! and prevent it from getting smaller than 1e-10
         Err           = SQRT( Err / NonDummySpc_Count )
         ros_ErrorNorm = MAX( Err, 1.0d-10 )
         RETURN

      ENDIF

      !~~~> Scalar Tolerance (same AbsTol & RelTol for all species)
      DO IDX = 1, NonDummySpc_Count
         I     = NonDummySpc_Indices(IDX)
         Ymax  = MAX( ABS( Y(I) ), ABS( Ynew(I) ) )
         Scale = AbsTol(1) + RelTol(1) * Ymax
         Err   = Err + ( Yerr(I) / Scale )**2
      ENDDO

      ! Normalize the error norm by the number of non-dummy species
      ! and prevent it from getting smaller than 1e-10
      Err           = SQRT( Err / NonDummySpc_Count )
      ros_ErrorNorm = MAX( Err, 1.0d-10 )
      RETURN

   ENDIF

   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   ! If optional arguments are not passed, fall back to the default,
   ! which is to include all species in the error norm computation.
   !
   ! Use separate loops for scalar and vector tolerances, as this
   ! will facilitate loop vectorization for better performance
   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

   !~~~> Vector Tolerance (per-species AbsTol & RelTol)
   IF ( VectorTol ) THEN

      DO I = 1, N
         Ymax  = MAX( ABS( Y(I) ), ABS( Ynew(I) ) )
         Scale = AbsTol(I) + RelTol(I) * Ymax
         Err   = Err + ( Yerr(I) / Scale )**2
      ENDDO

      ! Normalize the error norm by the number of non-dummy species
      ! and prevent it from getting smaller than 1e-10
      Err           = SQRT( Err / N )
      ros_ErrorNorm = MAX( Err, 1.0d-10 )
      RETURN

   ENDIF

   !~~~> Scalar Tolerance (same AbsTol & RelTol for all species)
   DO I = 1, N
      Ymax  = MAX( ABS( Y(I) ), ABS( Ynew(I) ) )
      Scale = AbsTol(1) + RelTol(1) * Ymax
      Err   = Err + ( Yerr(I) / Scale )**2
   ENDDO

   ! Normalize the error norm by the number of non-dummy species
   ! and prevent it from getting smaller than 1e-10
   Err           = SQRT( Err / N )
   ros_ErrorNorm = MAX( Err, 1.0d-10 )

Related GitHub issues

• Closes [BUG/ISSUE] Prod/loss species are contributing to the Error Norm when they should be ignored #66
• Closes Adding dummy species into reactions (to track P/L rates) changes results (in KPP 3.3.0) #144
• See Adding inert reaction tags causes changes to chemistry in 14.6.3 geoschem/geos-chem#3175

Tagging @RolfSander @msl3v @jimmielin @obin1 @fsolmon @ujongebloed @mje1398

[RolfSander]

It's good to see that there is now a possibility to exclude dummy
species from the Rosenbrock error calculation! I do have a suggestion
though:

With the current implementation, you can use the optional argument only
for the four selected integrators. When you switch to another integrator
and forget to remove the optional argument again, the compiler will
produce an error. Of course, we could add the new optional arguments to
all integrators but that's a bit tedious.

My suggestion is that we put NonDummySpc_Count into ICNTRL(19). The
calculation of NonDummySpc_Indices could be done by KPP automatically at
the end of SUBROUTINE Initialize (which is generated by
GenerateInitialize() in gen.c).

Would that work?

[yantosca]

Thanks @RolfSander, I think that's a good idea. I'll try it.

[yantosca]

@RolfSander: I just noticed that there is the #DUMMYINDEX option that lets you set the indices of species that aren't in the mechanism to zero. I was thinking I should rename e.g. NonDumSpc_Count etc. so as not to confuse our nomenclature. Maybe I should refer to these as "diagnostic species" (e.g. NonDiagSpc_Count, etc)? Or maybe "actual" species?

[RolfSander]

NonDiagSpc_Count

Yes, we should avoid confusion. Unfortunately, diag occurs quite often
in the KPP code, meaning diagonal, not diagnostic.

ActualSpec_Count would be ok.

Here is yet another idea: During the discussion, @obin1 used the name
"P/L species". So maybe ProdLossSpc_Count?

[obin1]

I like ProdLossSpc_Count! This is a good name if we only use them for the production/loss diagnostics.

If we are looking for something more general, would PassiveSpc_Count work? I would imagine we want all passive species (not necessarily inert, but passive diagnostic species) to be excluded from the adaptive timestepping strategy.

[yantosca]

Thanks @RolfSander. I think NonPassive_SpcCount etc is more general. I'll try that.