additional history output for Newton-Krylov Spinup

[wwieder]

This is related to #1451, but we're discussing using the Newton-Krylov spinup used for ocean BGC tracers in MIMICS too.
This method, however, could also be useful for spinning up soil C pools with the CENTURY soilBGC code and may be easier to maintain than the matrix approach in the long term (mainly because it doens't really touch the existing CTSM code and uses a seperate python package that @klindsay28 has developed for the ocean (see links below).

It sounds like two things are needed for this

1. .yml file with names of pools to be read off history files (SOILC1, SOILC2, etc) This is simple.
2. CTSM restart files with pools of interest (this is already done)
3. Additional history file output for the preconditional, which is where the magic occurs in helping to accelerate spinup.

This issue deals with the last requirements, # 3, for the preconditional, which is supposed to be the tendency of the first derivative for each pool. It seems like this can be obtained by writing out some additional fields to history files that will be used in spinup. These include the:

• decomp_k matrix, which is already dimensioned by pool (e.g. decomp_k_lit1, decomp_k_som1, decomp_k_cwd)
• pathfrac_decomp_cascade this is currently dimensioned by transitions, but will have to be dimensioned for each flux (pathfrac_l1_s1, pathfrac_s1s2). This will be pretty straight forward for century, but more involved for mimics with three fluxes going from m1_s* pools.
• Average states for each pool, this is already done.

One nuance here is that these history fields should be the time average for time period being iterated over for spinup (e.g. 20 years in a typical CLM spinup, but maybe only over a year or two for testing. This approach will also omit any information about advection and diffusion terms in the preconditional, but we'll see how many Newton-Krylov iterations we need for steady-state without this information, as the land model is inexpensive enough to run that these terms may not speed up our solution much.

@slevisconsulting since we're talking about doing this first with the century-bgc code, it seems like this should go into a different PR off the latest CTSM tag? Does this make sense @ekluzek?

@klindsay28, please suggest additional though, corrections, or ideas here to correct, clarify, or improve these notes.

---

Links to the Newton-Krylov tools @klindsay28 has developed.
The Newton-Krylov based spinup tool that I've put together is on github at
https://github.com/klindsay28/Newton-Krylov_OOC
with documentation on Read the Docs at
https://newton-krylov-ooc.readthedocs.io/en/latest/index.html

For our conversation, I think it'd be useful if you could read the background and terminology pages from the Description portion of the documentation:
https://newton-krylov-ooc.readthedocs.io/en/latest/description/background.html
https://newton-krylov-ooc.readthedocs.io/en/latest/description/terminology.html

Keith gathered NK-related documents here

[wwieder]

@klindsay28 and @slevisconsulting how do we handle the parallel N pools and tendencies for the N-K solver. This is maybe less of an issue with the BGC code, which has fixed stoichiometry. But I wondered how you handle the other tracers in POP?

[klindsay28]

There are 6 dissolved organic matter tracers in POP that are somewhat analogous to soil tracers in CLM. There are 2 tracers each for C, N, and P. Remineralization timescales are different for each element, and for each of the 2 tracers for that element. The N-K solver is applied to all 6 tracers individually. There is a preconditioner matrix for each of the 6 tracers.

If stoichiometry of SOM tracers is fully fixed, including source and sink terms, then it doesn't seem necessary to apply the solver to each element. If the SOM source terms differ, but the derivatives d(SOM_remin)/d(SOM) are the same for each element, then the solver can be applied to each pool. The solver could use the same preconditioner matrix for each element in that scenario. If these derivatives are different for each element, then you would use element specific preconditioners.

[wwieder]

I think we're fine to ignore N pools in CENTURY output, as it passively follows C through the decomposition cascade.

I don't think the stoichiometry changes much over time in MIMICS either, so assuming we can just update the soil N states based on the updated soil C states and initial C:N ratios that go into the solver we can likely ignore solving for nitrogen in the N-K solver?

For our meeting tomorrow would it be helpful to have dummy data to start looking at @klindsay28, or can you advise us on next steps with information from @slevisconsulting's work on #1457?

[wwieder]

@klindsay28 the second 20 year run with restart and tendency files is here /glade/scratch/wwieder/archive/RandomBoreal_1.
The first 20 year simulation is here /glade/scratch/wwieder/archive/RandomBoreal_0

The case directory is here /glade/work/wwieder/ctsm/ctsm5.1_N-K_test/cime/scripts/RandomBoreal_1
I've also started an AD and postAD spinup for this point that we can use to compare to results generated with N-K.

Is it helpful to discuss what's in these files before you dive into the N-K code?

[klindsay28]

Having a chat about the output would be helpful to me. I'm looking to ensure that I understand the relationship between the output variables and equations 21.1 and 21.2 in the Decomposition chapter of the CLM5 Tech Note. How about we chat after the CLM Meeting today?

[wwieder]

@slevisconsulting and @klindsay28 have started working on this, and I wondered if we need a longer-term vision for maintaining spinup capabilities for different soil BGC models? This likely warrants a longer conversation, but I thought I'd introduce it here.

The topic is relevant in the short term to see if we should just the N-K working in MIMICS, or if it should be tested with CENTURY too? Assuming the N-K approach works efficiently with MIMICS, do we also see this as a replacement for AD mode or the MATRIX in CENTURY simulations (which would facilitate code maintainability), or is it redundant to have 3 different ways to spin-up the CENTURY model?

[slevis-lmwg]

My thoughts:

1. A potentially helpful advantage of getting N-K to work for CENTURY (rather than just for MIMICS) is the ability to evaluate the N-K spin-up against the other methods. Ideally this should give us the confidence of a good starting point before adapting N-K for MIMICS.
2. I doubt that it adds much effort (and hope that it doesn't) if @klindsay28 and I do this work for CENTURY first, rather than doing only for MIMICS.

[ekluzek]

@wwieder I agree with @slevisconsulting. Technically you are right that it is redundant to have 3 different ways to spin-up the CENTURY model, but having the three in place allows you to compare and contrast them and decide on which is the best to keep. And it's also possible that when they are in place -- you'll find ways of combining them together that you hadn't thought of now. I suspect that it will take enough development time that we will want to have all three in place for some time before we may get rid of any one of them. The MATRIX approach does still require AD mode, so you can't get rid of if.

The MATRIX and AD mode also is a good example of using the two spinup modes together. Originally it was thought that MATRIX spinup would replace AD mode, but we get the shortest spinup times using the two together. It's possible there will be something like that with the N-K mode, that you get the shortest spinup times in some combination with MATRIX and/or AD spinup mode. So having all three in place at the same time will facilitate this type of activity.

[wwieder]

OK, let's go ahead with the plan to implement N-K with CENTURY first and then adapt it to MIMICS after we learn a bit more, @slevisconsulting .