walkeranthonyp
commented
5 years ago Nice one @rgknox , this is great! A few initial comments:
-
From talking to Nate, the solver in FATES is a fixed point iteration. The difference between the current guess and the old guess is evaluated against the convergence criteria. Once this difference is smaller than the convergence criteria the current guess is deemed the solution. A small note: apparently this is not a residual, residual in solver terminology refers to a function constructed for other numerical methods where the residual = f(x) - x or x - f(x) and x is what you are trying to solve.
-
The normalisation Ci / p * 1e6 is converting Ci in Pa to umol mol-1 (so solved values should be of the rough order 100-1000). The second normalisation Ci / Ca is the Ci:Ca ratio (so will give solved values of the rough order 0.5 - 0.9). So the convergence criteria of 1e-6 is extremely strict, especially when the difference is in units of umol mol-1.
-
Given that extremely strict criteria and the distribution of the difference I would say that the fact that the solution is not within the tolerance is probably not important. At least in this specific BCI test case. The size of this difference might change in other locations with different parameters or environmental variables.
-
Also I wonder what level of precision we need in the photosynthesis solve so that the canopy temperature solver will converge? Seems like that converges in 5 - 10ish loops. But again we have to be careful that this is at BCI, for this PFT, etc.
-
I'm working on a semi-analytical solver that doesn't need any iteration. It finds only an approximate solution but it's pretty damn close - see point above about what level of precision is needed. ELM and CLM4.5 and up use Brent's method, coded up by Jinyun I believe. We could use my method to come up with pretty tight brackets for the Brent method if people are not comfortable with an approximate solution (working on looking at how approximate right now). That would speed things up if Brent's method takes a substantial number of iterations to converge. Also Nate's Newton-Raphson method might work well. Or using previous solutions as the starting guess for the fixed-point iteration.
-
It would be interesting to know how many iterations Brent's method takes to converge in ELM or CLM and how many iterations the canopy temperature solve takes. Presumably about the same number as you show in your plots.
rgknox
dlawrenncar
alistairrogers
@walkeranthonyp, Nathan Collier and I met up at the DOE Modeling PI meeting and had a fun discussion about improving the photosynthesis code. @walkeranthonyp has some more fundamental ways to improve how we go about this process regarding which schemes we use and how the equations themselves are cast, and Nate was interested in how we actual perform out solve for C_i (interstitial CO2), which is an iterative loop converging on a solution to non-linear algebraic equations.
See here where we perform the iteration: https://github.com/NGEET/fates/blob/master/biogeophys/FatesPlantRespPhotosynthMod.F90#L992
An important note here, is that the photosynthesis scheme calculates stomatal conductance (among other things), which in-turn effects latent-heat transfer from the surface, which effects land-atmospheric stability, which should tightly feed-back on these fluxes again. Thus, photosynthesis is called inside this broader iteration loop (CLM/ELM side of code) that seeks to converge on a stability parameter and land-fluxes with a low residual (change since last iteration). Note: These coupled iteration loops make the search for C_i (and the quadratic root finders) the most often called subroutines in the model. We cap this loop iteration at 5 steps currently.
I wanted to diagnose how efficient our code is currently performing the solve for C_I, and see if there is any low-hanging fruit that can be picked to help speed things up and make it more efficient.
The first thing I did was do some counting statistics to see the distributions of how many iterations are performed. The following plot shows results from a 1 year simulation at BCI, using a forest inventory initialization.
Top Left Panel: Histogram of counts needed to converge on the outer loop (stability parameter) Top Right Panel: Histogram of the final residuals in the C_i loop, using the existing formulation. The red-line is the residual error criteria for success Bottom Left Panel: Histogram of an ALTERNATIVE way to formulate the residual in the C_i loop, assuming that the residual is the change in C_i from previous time-point, normalized by the canopy CO2 partial pressure. Bottom Right Panel: Histogram of the number of inner loop counts for C_i (which is capped at 5. Spoiler alert, almost every single iteration attempt used up 5).
So all iterations used up the max 5 attempts. Also, it seemed that the residual is strangely formed. Why is it the difference in C_i normalized by air-pressure and not CO2 partial pressure which should have the same units? I can't rectify the units on the current formulation, it is pascals of partial pressure divided by atmospheric pressure * 1 million? (maybe they thought that C_i was in units of PPM?) But it seems incredibly strict, the residuals at completion (due to step-count maxed) are all way above the error limit.
I was curious if there would be a significant shift in how many iterations were required, if I changed both the way the residual is formed, and gave it a lower success threshold of 1e-5. The residual that is evaluated here is the change in C_i divided by air CO2 partial pressure (which are the same units Pa, and similar magnitudes).
My take home from this exercise is that by adjusting the form of the residual there is a very modest effect on how quickly it converges. Keep in mind here that the physics was not changed, nor was the algorithm changed that searches for a solution.
Nate Collier did some calculations, and thinks that with Newton-Raphson, a solution should be found within 2-3 steps (I think), so I think that would be a worthwhile undertaking, as this would speed up the code. I have to look at the code more thoroughly, to see what method we use now (bracketing?) but it is not Newton-Raphson.