Dev doeclim #206

cahartin · 2017-11-22T02:43:13Z

This PR includes the a new temperature component (DOECLIM) to better represent ocean heat uptake and global mean temperature change.

This document contains more detailed information on the changes.

Project file to Xcode 6.; C++ runtime used is now compiler default (clang) instead of GNU C++; cleaned up header search paths; added linker flags for boost `filesystem` and `system`. This allows building on Mac OS X 10.11, Boost 1.60, and GSL 2.1.

linking boost and gsl libraries in both hector and hector-lib project files

Fix VS project

See issue #136. This commit changes `main.cpp` so it returns nonzero exit codes (1=program exception, 2=standard exception, 3=other exception).

- 'SOx', 'NH3', 'C6F14', 'HALON1202'

Xcode project file updates

Remove input variables which are not being used

Rc1.1.3

boost (header only). This opens the door to release hector under different licenses that GPLv2.

Replace GSL algorithms including ODE integration with various libs from

information instead of including the license in each file. This will make less version control headache to release under different licenses.

Have header and source files reference LICENSE.md for licensing

instead just use the logger with cout echo turned on. This closes #123.

Remove duplicative cout statements in the carbon-cycle-solver

…windows for that matter) the Xcode and xcode directories get merged as one causing extra confusion.

Remove extra visual studio and Xcode project files.

bpbond

This looks like a great improvement to Hector, and nice job on the authorea document! Need more clarity in some spots, and more defensive programming in others (e.g. use of const).

bpbond · 2017-11-22T16:50:48Z

headers/components/component_names.hpp

@@ -26,6 +26,7 @@
 #define OH_COMPONENT_NAME   "OH"
 #define N2O_COMPONENT_NAME "N2O"
 #define TEMPERATURE_COMPONENT_NAME "temperature"
+#define TEMP_DOECLIM_COMPONENT_NAME "temp_doeclim"


Should we use a more meaningful name?

Ah, the "DOE" here doesn't refer to the department of energy. OK then.

@bpbond and @rplzzz yea it's diffusive ocean energy climate model. are there any licensing issues we need to be aware of of linking this one model with Hector? Can we still call it Hector?
https://github.com/scrim-network/BRICK

We are using DOEClim with the permission of its authors, right? If so, then it's fine. We should, however, acknowledge the code they donated. I have a note about that elsewhere.

I have not done any investigations into permissions. DOECLIM has been used several times by Klaus Keller's group and I adapted their C++ implementation here. I'll check in with Ryan and Klaus about it.

It looks like one or more versions of DOECLIM have been released under GPLv3. Just get verbal confirmation from one of the original authors that they're ok with us releasing it under GPLv2. (The differences aren't really relevant to what we're doing; I like v2 because it's shorter and written in plain language.) Then we can just indicate that the code is used with permission.

renamed component. Still need to double check licensing.
(CH - tracking her changes)

A short update on the permissions: we now have emailed permission from Elmar Kriegler (original creator of DOECLIM) to use DOECLIM in Hector under GPLv2. I'm waiting for a response from Klaus Keller (coauthor of the paper for which the GPLv3 version was written).

We now also have emailed permission from the author of the DOECLIM version released under GPLv3

bpbond · 2017-11-22T16:51:16Z

headers/components/ocean_component.hpp

-	unitval k_max;           //!< maximum ocean heat uptake efficiency, W/m2/K
-    unitval k_min;          //!< minimum ocean heat uptake efficiency, W/m2/K
-    unitval slope;           //!< slope of the curve, 1/K
-//	unitval cp;             //!< Specific heat


Why is this commented out? Remove?

thanks, i missed that one.

bpbond · 2017-11-22T16:53:20Z

headers/components/temp_doeclim_component.hpp

+    double ak;               // slope in climate feedback - land-sea heat exchange linear relationship
+    double bk;               // offset in climate feedback - land-sea heat exchange linear relationship, W/m2/K
+    double csw;              // specific heat capacity of seawater W*yr/m3/K
+    double earth_area;       // m2


Seems like many of these should be const

instead of doubles?

No, I mean, something like const double earth_area = x.xxx;

Based on @rplzzz 's comment below it looks like we should move these initializations to the header (some are called in more than one function so if we want to keep them together they need to be there). Those with well-known values will have const static, while those that are potentially tweakable will have const. Does that sound like what we want?

bpbond · 2017-11-22T16:54:06Z

headers/components/temp_doeclim_component.hpp

+    double A[4];
+    double IB[4];
+
+    //Dependent DOECLIM parameters


Space after //. Dependent on what? ~~What's the difference between ocean_area below and the same variable name above?~~ misread

@bpbond If these are just consts then we don't actually need them in the header?

Mmm. Feeling fuzzy on my C++. I think actually we'd want static const x = .... As for header or not, it depends; is it referred to in multiple component functions? If not, right, just move into the run method or wherever.

const means that once it's initialized, it can't be changed. static means that there is one copy for the entire class, instead of one for each object created. A member can be either const, static, neither, or both, depending on what you are trying to accomplish. For physical constants like the area of the earth or the specific heat of seawater, static const is the way to go. For something that might be a settable parameter (ak, perhaps?), you would generally do that as a normal (i.e., not static) member.

A static member is initialized just like a global variable, irrespective of whether or not it happens to be const. In fact, static members are basically just namespaced global variables, except that C++ does respect their private or protected visibility.

A member that is const, but not static must be initialized in the constructor. This is one of the rare cases in which the distinction between initialization and assignment matters. In this case you have to use the : notation to make it an initialization instead of an assignment. If you can't do that for some reason, such as because you have to construct the object before the initialization data is available, then you can't make the members const. Your best bet in this case is to make them ordinary members and to add a comment to the effect that once they are set they should be treated as const.

A constant that is only used in one method could be moved into that method, as Ben suggests, but sometimes it's best to make them const or static const. The reason why is that if someone wants to know, what are they using for the xyz parameter, it's a lot easier to find it if it's with all the other constants, rather than tucked away at the beginning of some particular method they might not know the name of. Also, if you end up refactoring so that more than one method uses the constant, then you'll have to make the constant a member anyhow. Starting it off that way future proofs your design a bit.

bpbond · 2017-11-22T16:58:55Z

source/components/temp_doeclim_component.cpp

+ *  \param[in] y        Inverted 1-d matrix
+ *  \returns            void, inverse is stored in y
+ */
+void TempDOECLIMComponent::invert_1d_2x2_matrix(double * x, double * y) {


Since we depend on Boost anyway, why not use its ublas library?

Any Boost inversion implementation I saw in my cursory google search looked a good deal more complex than this simple algorithm for 2x2 inversion. But I could be wrong!

That's correct. If you know the matrix is 2x2, you can invert it very simply using analytic formulae. It's similar to how you could solve an equation using the Newton-Raphson method, but if you happen to know it's quadratic, it's simpler to use the quadratic formula.

bpbond · 2017-11-22T17:01:02Z

source/components/temp_doeclim_component.cpp

+    tauksl  = (1.0 - flnd) * cas / kls;
+    taukls  = flnd * cal / kls;
+
+    // First order


Some more documentation needed here--what's going on?

@bvegawe can you add a few comments here? thank you!

My guess is that this has all been explained in a paper somewhere. We should include a pointer to the paper and document any mappings between names used in the paper and names used here that aren't obvious.

bpbond · 2017-11-22T17:01:49Z

source/components/temp_doeclim_component.cpp

+    //// to match a temperature record), need to track the Ftot that *would* have
+    //// produced the observed temperature record. This way there's a smooth
+    //// transition when we exit the constraint period, after which internal_Ftot
+    //// will rise in parallel with the value reported by ForcingComponent.


Nice commenting and logic.

@bpbond i think that was you back in the day!

bpbond · 2017-11-22T17:04:47Z

Note also my comment about using Boost–I doubt we want to roll our own matrix-inversion functions, but will defer to @rplzzz here.

rplzzz

Lots of cleanup needed in various places, though most of it shouldn't be too much work. I did have a question about ocean heat flux that I would like to see some discussion on, but everything else is minor.

rplzzz · 2017-11-27T20:59:25Z

input/hector_rcp85_doeclim.ini

@@ -0,0 +1,401 @@
+; Config file for hector model: RCP85


What is the story with this file? We have already added DOEclim to the rcp85.ini file above, so why is this one necessary?

good point. this was our initial test case file.

rplzzz · 2017-11-27T21:03:08Z

headers/components/temperature_component.hpp

@@ -1,87 +0,0 @@
-/* Hector -- A Simple Climate Model


I really don't love that we are deleting the original temperature component entirely. The whole idea behind Hector's design is that we are supposed to be able to add and drop components as needed. I get that we don't want components to proliferate because then we have to maintain them all, but so far maintaining components hasn't been a huge amount of effort. Is the original temperature component really unable to coexist with the DOEClim temperature (in the code base, I mean; obviously only one of them can be active in any single run)?

In my experience, it works fine with both components and just specifying which version you want to use in the .ini file. It does lead to a couple warnings like below:
WARNING:registerCapability: temperature is declaring capability Tgav previously registered

They can exist together and do currently. My concern is that the old way was a hack with a bunch of bugs that were found. I don't think I want people having an easy option to get back incorrect results by turning on that old component.

That sounds reasonable. Since that leaves the doeclim temperature component as the only temperature component, maybe we should just rename it to "temperature_component". (We would still explain its doeclim origins in comments in the source file, of course.)

rplzzz · 2017-11-27T21:06:46Z

source/components/ocean_component.cpp


 	lastflux_annualized.set( 0.0, U_PGC );

    // Register the data we can provide
    core->registerCapability( D_OCEAN_CFLUX, getComponentName() );
 	core->registerCapability( D_OCEAN_C, getComponentName() );
-    core->registerCapability( D_HEAT_FLUX, getComponentName() );


I take it these heat flux capabilities have been moved to the DOEClim component?

Edit: Actually, it doesn't look like they have. Are we just not supporting retrieving the ocean heat flux anymore, or is it going to known by a different name going forward? Apart from being a serious backward-incompatible change, it would be kind of unfortunate if we could no longer get this information from the model.

With DOECLIM, we calculate a flux into the mixed layer and a flux into the interior. For backward compatibility, I'll make a note to also do a total heat flux.

Heat flux is calculated in the DOECLIM component now. @bvegawe that's a good idea to have total heat flux into the ocean.

rplzzz · 2017-11-27T21:08:25Z

source/components/temp_doeclim_component.cpp

+/*! \brief              Calculates inverse of x and stores in y
+ *  \param[in] x        Assume x is setup like x = [a,b,c,d] -> x = |a, b|
+ *                                                                  |c, d|
+ *  \param[in] y        Inverted 1-d matrix


This should be documented as param[out] or param[inout], depending on whether the function expects y to be initialized.

rplzzz · 2017-11-27T21:10:25Z

source/components/temp_doeclim_component.cpp

+/*! \brief              Calculates element-wise sum of matrix x and y, stores in z
+ *  \param[in] x,y      Assume x,y setup like x = [a,b,c,d] -> x = |a, b|
+ *                                                                 |c, d|
+ *  \param[in] z        Summed 1-d matrix


This also is param[out].

rplzzz · 2017-11-27T21:40:23Z

source/components/temp_doeclim_component.cpp

+    //
+    // If the user has supplied temperature data, use that (except if past its end)
+    if( tgav_constrain.size() && runToDate <= tgav_constrain.lastdate() ) {
+    //    H_LOG( logger, Logger::NOTICE ) << "** Using user-supplied temperature" << std::endl;


Delete this line. Saving an old version of an error message in a commented line is perhaps the most egregious example of code hoarding I can think of.

@bpbond @rplzzz we should probably get this constraint working again in the new component. for now i'll leave the comments in and we can test this later.

rplzzz · 2017-11-27T21:40:56Z

source/components/temp_doeclim_component.cpp

+    //    H_LOG( logger, Logger::NOTICE ) << "** Using user-supplied temperature" << std::endl;
+        H_LOG( logger, Logger::NOTICE ) << "** ERROR - DOECLIM can't currently handle user-supplied temperature" << std::endl;
+    }
+    //  tgav = tgav_constrain.get( runToDate );


Lots of commented-out code from here to line 385. Clean it up.

rplzzz · 2017-11-27T21:44:37Z

source/components/temp_doeclim_component.cpp

+    double aero_forcing = double(core->sendMessage( M_GETDATA, D_RF_BC ).value( U_W_M2 )) + double(core->sendMessage( M_GETDATA, D_RF_OC ).value( U_W_M2 )) + double(core->sendMessage( M_GETDATA, D_RF_SO2d ).value( U_W_M2 )) + double(core->sendMessage( M_GETDATA, D_RF_SO2i ).value( U_W_M2 ));
+    forcing[tstep] = double(core->sendMessage( M_GETDATA, D_RF_TOTAL ).value( U_W_M2 )) - ( 1.0 - alpha ) * aero_forcing;
+
+    // Initialize variables for time-stepping through the model


Are the next 80 lines or so copied/translated from the original DOEClim code base? If so, we should maybe include a pointer to that code, both so that anyone who wants to look at it can, and also to document that it was copied with the knowledge and permission of its authors and maintainers.

rplzzz · 2017-11-27T21:51:24Z

headers/components/component_names.hpp

@@ -26,6 +26,7 @@
 #define OH_COMPONENT_NAME   "OH"
 #define N2O_COMPONENT_NAME "N2O"
 #define TEMPERATURE_COMPONENT_NAME "temperature"


If we are removing the old temperature component, then shouldn't we also remove this #define?

rplzzz · 2017-11-27T21:55:11Z

headers/components/component_names.hpp

@@ -26,6 +26,7 @@
 #define OH_COMPONENT_NAME   "OH"
 #define N2O_COMPONENT_NAME "N2O"
 #define TEMPERATURE_COMPONENT_NAME "temperature"
+#define TEMP_DOECLIM_COMPONENT_NAME "temp_doeclim"


We are using DOEClim with the permission of its authors, right? If so, then it's fine. We should, however, acknowledge the code they donated. I have a note about that elsewhere.

rplzzz

In addition to what I wrote before (which seems to be interspersed with the line comments for some reason, so you have to hunt for it a bit), here are two more things.

Can we please get rid of this junk:

Wed Nov 22 02:46:24 2017:NOTICE:run: Attempting ODE solver 1878->1879 (1878->1879)
Wed Nov 22 02:46:24 2017:NOTICE:run: Success: we have reached 1879
Wed Nov 22 02:46:24 2017:NOTICE:run: ODE solver success at t= 1879  last dt= 0.5

As far as I can tell, it's useful only when something goes wrong, so we should downgrade it to debug level and make notice level output the default.

The code is failing its tests, presumably because we've changed the model, so naturally the outputs no longer match. Once we are satisfied with the results using the new model, we need to update the comparison data for the tests.

additional descriptions of model parameters and calculations. added a total ocean heat flux output parameter to match the parameter output by the original temperature component

Dev doeclim: improved DOECLIM descriptions, total heat flux output

…v_doeclim

renamed temp_doeclim to temperature in code and ini files updates to header to include a few more parameters not specified

rplzzz

My biggest comment has been addressed. There were a couple other things, but the world won't come to an end if they're not addressed, and we should move this thing along.

Forgot to add the factor for the slightly smaller surface area of the interior ocean

includes full citations and the fact that we have permission.

set hard-coded DOECLIM parameters to const; moved their initializations to the header. if these initializations are static I get a warning, so just const for now.

dev_doeclim: a bug fix, const parameters, updated references

bpbond and others added 30 commits January 19, 2016 07:35

Xcode project file updates

bd1ae89

Project file to Xcode 6.; C++ runtime used is now compiler default (clang) instead of GNU C++; cleaned up header search paths; added linker flags for boost `filesystem` and `system`. This allows building on Mac OS X 10.11, Boost 1.60, and GSL 2.1.

VS project files

99cc46c

linking boost and gsl libraries in both hector and hector-lib project files

fixed compile error in ini_to_core.cpp

2c2e142

Merge pull request #135 from JGCRI/fix-VC-project

ae397ad

Fix VS project

Return nonzero exit codes on error

af19e2f

See issue #136. This commit changes `main.cpp` so it returns nonzero exit codes (1=program exception, 2=standard exception, 3=other exception).

Remove input variables which are not being used

c19f847

- 'SOx', 'NH3', 'C6F14', 'HALON1202'

Merge pull request #134 from JGCRI/xcode-MacOS10.11

6b627f5

Xcode project file updates

Merge pull request #143 from rgieseke/remove-unused-variables

bf6499f

Remove input variables which are not being used

Merge pull request #145 from JGCRI/rc1.1.3

7f87172

Rc1.1.3

Replace GSL algorithms including ODE integration with various libs from

d98accf

boost (header only). This opens the door to release hector under different licenses that GPLv2.

Clean up commented out code in carbon-cycle-solver.

523874a

Update wiki submodule tag to reflect latest documentation updates.

e1954b9

Merge pull request #149 from JGCRI/replace-gsl

caea4fe

Replace GSL algorithms including ODE integration with various libs from

Have header and source files reference LICENSE.md for licensing

b88a5d2

information instead of including the license in each file. This will make less version control headache to release under different licenses.

Merge pull request #154 from JGCRI/consolidate_license

9630601

Have header and source files reference LICENSE.md for licensing

Remove duplicative cout statements in the carbon-cycle-solver and

c715b7a

instead just use the logger with cout echo turned on. This closes #123.

Merge pull request #155 from JGCRI/cc_solver_logging

5a3fbaf

Remove duplicative cout statements in the carbon-cycle-solver

Remove extra visual studio and Xcode project files. Note on OSX (and …

5459c00

…windows for that matter) the Xcode and xcode directories get merged as one causing extra confusion.

Merge pull request #156 from JGCRI/remove_extra_proj_files

46e79cd

Remove extra visual studio and Xcode project files.

Added test travis-ci script

58b9ee2

Test commit

f784a2e

Don't include BOOSTROOT in includes

855fb45

Explicitly download and built Boost (not apt-get)

0608330

Back BOOSTROOT back into Makefile...

b606fb8

Added b2 Boost build command

804a697

Added build status badge

c84f5f5

Comment tweaks to .travis.yml

d44e98a

Might as well test against latest Boost

7c816b3

Run RCP 4.5 after buliding

76edd94

Build using both gcc and clang

ea901d3

cahartin added the enhancement Changes to the Hector software not affecting the representation of the Earth system label Nov 22, 2017

cahartin added this to the v2.0 milestone Nov 22, 2017

cahartin requested review from bpbond and rplzzz November 22, 2017 02:43

bpbond requested changes Nov 22, 2017

View reviewed changes

Merge remote-tracking branch 'JGCRI/dev_doeclim' into dev_doeclim

14eeb24

rplzzz suggested changes Nov 27, 2017

View reviewed changes

delete redundant ini file

b478e76

bpbond mentioned this pull request Nov 29, 2017

Automated calibration of Hector against time series #208

Closed

cahartin changed the base branch from master to rc2.0.0 November 29, 2017 18:30

rgieseke mentioned this pull request Nov 30, 2017

Error: bad Tk value #207

Closed

bvegawe and others added 8 commits November 30, 2017 15:35

DOECLIM descriptions and heatflux output

f8a9ac3

additional descriptions of model parameters and calculations. added a total ocean heat flux output parameter to match the parameter output by the original temperature component

Merge pull request #209 from bvegawe/dev_doeclim

7c24918

Dev doeclim: improved DOECLIM descriptions, total heat flux output

Merge branch 'dev_doeclim' of https://github.com/JGCRI/hector into de…

bdcfac8

…v_doeclim

Renamed temp_doeclim

c92f9e5

renamed temp_doeclim to temperature in code and ini files updates to header to include a few more parameters not specified

renaming temp doeclim to simply temperature

379c9b2

renaming temp doeclim to temperature component

f558b90

VS project file - rename temp doeclim to temperature

97e0e8a

one more VS proj file

c424684

rplzzz approved these changes Jan 3, 2018

View reviewed changes

bvegawe added 3 commits January 3, 2018 12:18

correct bug in total heatflux calc

73a8bb4

Forgot to add the factor for the slightly smaller surface area of the interior ocean

improved DOECLIM references

848ef26

includes full citations and the fact that we have permission.

const DOECLIM params

bd4d80d

set hard-coded DOECLIM parameters to const; moved their initializations to the header. if these initializations are static I get a warning, so just const for now.

cahartin added this to To Do in checklist for v2.0 Jan 5, 2018

Merge pull request #211 from bvegawe/dev_doeclim

448e65a

dev_doeclim: a bug fix, const parameters, updated references

cahartin merged commit b5e7c3b into rc2.0.0 Jan 5, 2018

cahartin moved this from To Do to Done in checklist for v2.0 Jan 8, 2018

rplzzz deleted the dev_doeclim branch June 20, 2018 13:44

Dev doeclim #206

Dev doeclim #206

Conversation

cahartin commented Nov 22, 2017 • edited

bpbond left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

bvegawe Dec 14, 2017 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

bpbond Nov 22, 2017 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

bpbond commented Nov 22, 2017 • edited

rplzzz left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

rplzzz left a comment • edited

Choose a reason for hiding this comment

rplzzz left a comment

Choose a reason for hiding this comment

cahartin commented Nov 22, 2017 •

edited

bvegawe Dec 14, 2017 •

edited

bpbond Nov 22, 2017 •

edited

bpbond commented Nov 22, 2017 •

edited

rplzzz left a comment •

edited