Unexpected (output) behaviour when using dependencies on aggregated variables

[FabianGrosze]

During model development, I recently noticed the following behaviour when using (user-defined) standard variables with aggregate_variable = .true. and dependencies on such variable.

In my specific case I have three algae functional types and a standard variable:
type(type_qsim_variables) :: qsim_variables

using:

type :: type_qsim_variables
   type(type_interior_standard_variable) ::                                                                              &
      algae_biomass = type_interior_standard_variable("algae_biomass", "mg/l", aggregate_variable = .true.)
end type type_qsim_variables

For each algae FT, the model does:
call self%add_to_aggregate_variable(qsim_variables%algae_biomass, self%id_algae_biomass),
where self%if_algae_biomass points to the biomass of the individual group.

The model has one zooplankton group whose grazing depends on the total algal biomass (i.e. qsim_variables%algae_biomass). Due to grazing, the contributing terms (i.e. the individual algal groups' biomasses) are changed via _ADD_SOURCE_(...) by the zooplankton model.
If the zooplankton model uses register_(state_)dependency(self%id_algae_biomass, qsim_variables%algae_biomass), i.e. it uses the aggregated variable, the values of qsim_variables%algae_biomass written to the output have a lag of one time step compared to the individual contributing terms, i.e. at least for the output:
qsim_variables%algae_biomass(t+1) = algae_biomass_1(t) + algae_biomass_2(t) + algae_biomass_3(t).

Note that register_dependency vs. register_state_dependency makes no difference.

If I use the same dependency in a model that does not change any of the contributing variables, the lag does not occur, i.e. qsim_variables%algae_biomass(t) = algae_biomass_1(t) + algae_biomass_2(t) + algae_biomass_3(t).

I am curious to know how this behaviour can be explained and whether this is intended or not?

Cheers,
Fabian

P.S. For our zooplankton model, this behaviour can be avoided as the zooplankton model also requires the individual biomasses, so the dependency on the aggregated variable can be removed.

[jornbr]

Hi Fabian,

That sounds surprising. What is the version of FABM you are using, and what is the host model? And is self%if_algae_biomass a state variable in the algae FT? Or a diagnostic?

You can get detailed information about the order in which FABM calls the different model components by adding log: true at the top of fabm.yaml. If you then run, you should get a file fabm_jobs.log, among others, which describes the data loads/stores and calls made by FABM, and their ordering. It is typically a very verbose log file, but it could help understand when the aggregate total is updated.

Cheers,

Jorn

[FabianGrosze]

Hi Jorn,

Thanks for the hint re. log: true, I didn't know about that option (and it made me notice some other bugs).
I did a small test, one with the dependency, one without. Attached are the fabm_jobs.log files of both tests.

With register_state_dependency: fabm_jobs_with_rsd.log
Without register_state_dependency: fabm_jobs_without_rsd.log

Note the different location of /algae_biomass_calculator: do in the sequence.

I am using FABM 3.0. We are currently building the host model, which is essentially an offline host for the hydrodynamic (river) model we are using, i.e. it reads hydrodynamics from that model's output and calculates transport etc. on that basis.

self%id_algae_biomass is a state variable for the different algal groups.

Cheers,
Fabian

[jornbr]

Hi Fabian,

Thanks! From what I can tell, everything is behaving as I would expect: if the rotifer module requests the value of the total algal biomass, FABM will ensure that its value is computed before the rotifer module is called (fabm_jobs_with_rsd.log). If it is not requested by any other FABM module, then the only reason to calculate it is if the user wants it as part of output. Then, the calculation of total algal biomass is scheduled last, as part of the host call to finalize_outputs (fabm_jobs_without_rsd.log - note that FABM has determined that no other module is interested in the value of algae_biomass_calculator/result calculated by algae_biomass_calculator: do). Nevertheless, the total biomass is updated, just later in the call sequence. If the updated value is not seen in the output, then the likely reason is that the host does not call FABM's finalize_outputs before output is sent/written. Can that be the case?

Cheers,

Jorn

[FabianGrosze]

Hi Jorn,

Our BGC do_timestep routine has this sequence (in shortened pseudocode form):

subroutine do_timestep(self, current_time, time_settings)
   class(type_bgc_model_data), intent(inout) :: self 

   ### get and apply surface fluxes
   ### get and apply bottom fluxes
   ### get and apply interior fluxes

   call self%bgc_model%finalize_outputs()

   ### Retrieve interior diagnostics
   ### Retrieve horizontal diagnostics
end subroutine do_timestep

This is the only location where finalize_outputs is/was called. self%bgc_model is declared class(type_fabm_model), pointer.

In the main routine, the call sequence during the simulation is this (again in shortened pseudocode form):

### update weather
### update boundaries
### update hydrodnamics
call self%bgc_model%prepare_inputs()
### write output (initial state)

do step = 1, n_timesteps
   call self%bgc_model_data%do_timestep()     ! <-- this is where finalize_outputs is called (see above)
   call self%bgc_model_data%validate_state()

   ### do transport
   ### write output (state after time step)
         
   ! initialize next loop iteration 
   if (step /= n_steps) then
      ### update weather
      ### update boundaries
      ### update hydrodynamics
      call self%bgc_model%prepare_inputs()
   endif
enddo

For testing, I added another call self%bgc_model%finalize_ouputs() between ### do transport and ### write output but it doesn't change anything.

In the previously uploaded log files, finalize_outputs occurs in the same place, with the only difference that /algae_biomass_calculator: do is part of that call for the without_rsd case while it's not for the with_rsd case. This doesn't change for the with_rsd case when adding the second finalize_outputs.

Provided I interpret the logs correctly, the model only does the aggregation once for eacht aggregation variable and it does it at the first required necessary location. However, perhaps something in the call sequence above is not as it should be to achieve the desired behaviour?

[jornbr]

It's puzzling... Regarding the code snippets above, the only thing that looks unusual is that you have the call to prepare_inputs at the end of the timestep. That seems unlikely to explain your issue, but it could cause stale values for FABM variables that depend directly on time (and then only if there is logic such as a time update in between the prepare_inputs and the do_timestep, not shown in the pseudocode). Is there any particular reason not to call prepare_inputs at the very start of do_timestep?

[FabianGrosze]

Hi Jorn,

There is a time update between do transport and write_output in the above pseudocode, which - based on your explanation - shouldn't matter, though, as prepare_inputs is the very last action in the loop and bgc_model_data%do_timestep is the very first one.

Re. prepare_inputs at the end of the loop: I can't give you the original reasoning for it as the colleague who started the host development has left. I can only assume, it was implemented this way because the initial state is written before entering the loop and so that it isn't done twice (for the first iteration).
I tested a few different sequences, among others, moving the calls from that step /= n_steps section (incl. ``prepare_inputs) to the start of the loop, so directly before do_timestep, but nothing changes.

Is there any way to enforce aggregation of all aggregation_variables before output? Also, one thing I noticed: for the initial state written to output, diagnostic variables and aggregation_variables still are zero. Is there a simple way to enforce their calculation before the first actual time step?

[jornbr]

The strange thing is that aggregate variables should be updated: they are by default marked for output and - if no other FABM module requests them - are therefore scheduled to be calculated by finalize_outputs. The log indicates that this happens as expected. What output variable (full name) shows the lag, and how do you access its values?

Regarding the vey first output: diagnostics are calculated somewhere in the sequence prepare_inputs -> get_xx_sources -> finalize_outputs. FABM decides where based on inter-module dependencies and computational efficiency. Thus, the typical approach is to call that full sequence (e.g. the full do_timestep) at the very start to ensure diagnostics are up to date - but without updating (time-integrating) the model state. A common approach is to have an equivalent of your do_timestep computing diagnostics and sources, and to later (after output) time-integrate the state based on those sources. That ensures that state and diagnostics are always fully in sync in output.

[FabianGrosze]

The output variable is a user-defined standard variable:

module qsim_standard_variables
   use fabm_types 
   implicit none
   private
   
   public :: qsim_variables

   type :: type_qsim_variables
      type(type_interior_standard_variable) ::                                                                              &
         algae_biomass     = type_interior_standard_variable("algae_biomass"    , "mg/l"    , aggregate_variable = .true.)
   end type type_qsim_variables

   type(type_qsim_variables), parameter :: qsim_variables = type_qsim_variables()

end module qsim_standard_variables

I am adding it to the output via algae_biomass in our settings.yaml containing all simulation settings. In the code I access it via (standard_variable = )qsim_variables%algae_biomass. This generally works, and it also works if I have a (state_)dependency in a model that doesn't change any of its contributors (i.e. any of the algae FTs' biomasses). It appears to me that, only because the rotifers modify the contributors via _ADD_SOURCE_, the mechanism somehow fails. This confuses me the most.

In fact, I will do another test with algae_biomass being used in another model without modifying its contributors and check the fabm_jobs.log.

[jornbr]

My guess is that it is the rotifers (or any other module) depending on algae biomass that makes the difference, because that forces the calculation of the algae biomass at an earlier (before-rotifer) stage - whether or not a module adds sources to one of the constituents of the aggregate variable should not matter, as that impact is essentially dealt with outside FABM by the host perfoming time-integration.

[jornbr]

One more thing to note: in the sequence

   ### get and apply surface fluxes
   ### get and apply bottom fluxes
   ### get and apply interior fluxes

I'm not sure what "apply" means, but I'd note it is not safe to change the state (e.g. time-integrate) there, as it will mean that some parts of prepare_inputs, get_xx_sources, finalize_outputs will see a different state than others. That will lead to inconsistent results since FABM is free to schedule calls anywhere in that sequence, and can change the schedule of any particular call when the model configuration (fabm.yaml) changes, as we've seen. Thus, in between prepare_inputs and finalize_outputs, no changes to state (or environment) should be made. The source terms would typically be saved there and applied to the state only after finalize_outputs. FABM's diagnostics (which include aggregate totals) will reflect the state that was set during the prepare_inputs ... finalize_outputs sequence.

[FabianGrosze]

I did some more tests, and indeed I had missed something in my test model (that - supposedly - depended on algae_biomass but didn't change any of its constituents). That said, I can confirm that it really is the register_(state_)dependency that moves up algae_biomass_calculator do in the calculation sequence.

I noticed one other thing that might be a clue for what's going wrong. The temporal lag also occurs for diagnostics, which, in fact, I suspect to relate to you addition above, because the code for ### get and apply interior fluxes currently looks like this:

      !------------------------------------------------------------------------
      ! 9.9.4 Interior source terms
      !------------------------------------------------------------------------
      allocate(int_source_terms(n_horizontal, n_interior_state), stat = status, errmsg = error_text)
      if (status /= 0) call fatal_error(proc_name, error_text)
      do i_vertical = 1, n_vertical
         int_source_terms = 0._rk
         call self%bgc_model%get_interior_sources(1, n_horizontal, i_vertical, int_source_terms)
         
         do ivar = 1, n_interior_state
            self%current%interior_state(:,i_vertical,ivar) = self%current%interior_state(:,i_vertical,ivar) &
                                                           + int_source_terms(:,ivar) * delta_t
         enddo 
      enddo

This seems a likely suspect for the entire behaviour - would you agree?

[jornbr]

Yes that is very likely the cause. I'd suggest declaring int_source_terms for the full domain (adding the i_vertical dimension), feeding int_source_terms(:,i_vertical,:) into get_interior_sources, and doing the interior_state update based on that after finalize_outputs with something like

self%current%interior_state = self%current%interior_state + int_source_terms * delta_t

If you further do the output in between finalize_outputs and the above state update, you'd get the nice result that state and diagnostics are in sync in output.

[FabianGrosze]

Jorn, thank you so much for letting me pick your brain!
After some code splitting and rearrangement following your suggestions I got it to work as intended!
All state, diagnostic and aggregation variables are finally in sync also in the output.

I have one last question (re. this work complex, possibly more to come in the further process). My entire run routine now looks like this:

subroutine run(self)
      class(type_simulation), intent(inout) :: self 

      integer             :: step, n_steps
      type(type_datetime) :: current_time
      
      call log_header("Main Loop", "=", centered = .true.)
      
      call log_message("Starting loop")
      
      n_steps = (self%time_settings%get_endtime() - self%time_settings%get_starttime()) / self%time_settings%get_timestep()
      do step = 1, n_steps

         current_time = self%time_settings%get_starttime() + self%time_settings%get_timestep() * (step - 1)
         call log_progress(step - 1, n_steps, current_time)
         
         call self%weather%update(current_time)
         call self%boundaries%update(current_time)
         call self%hydro%hydrodynamics%update(current_time)
         call self%bgc_model%prepare_inputs()

         call self%bgc_model_data%get_sources_sinks(current_time, self%time_settings)   ! <-- this includes finalize_outputs and get_diagnostics
         call self%bgc_model_data%obtain_conserved_quantities()
         call self%bgc_model_data%validate_state()
         
         call self%output%write_data(current_time, self%time_settings)

         call self%bgc_model_data%integrate_time(current_time, self%time_settings)
         call self%hydro%hydrodynamics%do_timestep(current_time, self%time_settings)
         call self%transport%do_timestep(current_time, self%time_settings)
      enddo

      ! get and write state at end of last time step
      current_time = self%time_settings%get_starttime() + self%time_settings%get_timestep() * n_steps
      call log_progress(n_steps, n_steps, current_time)
      call self%weather%update(current_time)
      call self%boundaries%update(current_time)
      call self%hydro%hydrodynamics%update(current_time)
      call self%bgc_model%prepare_inputs()

      call self%bgc_model_data%get_sources_sinks(current_time, self%time_settings)
      call self%bgc_model_data%obtain_conserved_quantities()
      call self%bgc_model_data%validate_state()
         
      call self%output%write_data(current_time, self%time_settings)

      call log_message("Main loop completed.")
   end subroutine run

So, the output written during the first iteration is my initial state with all diagnostic and aggregation variables being properly set/initialised.
To get the final state at the end of the simulation, I added that block of code after the loop. I am wondering if there is a more elegant/standard way to get both initial and final state (and obviously all intermediate ones) into the output?

[jornbr]

Great! Glad to hear it's resolved. That logic looks fine to me. If you'd want to compact it, I guess you could loop up to n_steps + 1 and exit the loop just after write_data if step == nsteps + 1 (or loop from 0 to nsteps and exit if step == nsteps). Alternatively, you could put the current_time ... write_data block in a dedicated subroutine.