umd-aosc / umd-sst Goto Github PK

add var.x to the ctests

• use BUMP covariance in the yaml
• use the sample observations that are already in this repository

@travissluka I reproduced your initial JEDI test, it's pretty fast on Cheyenne, taking about 5 minutes with 1 node (36 processors) to finished 1 month run.

Now I'll try to test for 0.25x0.25 grid, I've uploaded the initial background file "20080101_regridded_sst.nc" and the landmask.nc files at my feature/test branch. Please help check if they are ok when you get a minute.

use the land mask when initializing BUMP horizontal convolution.

Implement a simple "dumb" identity background error covariance in our Covariance class,

• similar to what is done here https://github.com/JCSDA-internal/fv3-jedi/blob/develop/src/fv3jedi/ErrorCovariance/ErrorCovariance.cc
• enable the the ErrorCovariance test, it should pass.

@travissluka I posted the following issue to at JCSDA/JEDI, I thought it would be good to post it here as well.

I am trying to use JEDI (and ECMWF/atlas) to implement the 3D-Var algorithm for SST data assimilation. The State and Observation will involve only one variable of SST. I have implemented my own "Geometry" class and now I am implementing the "State" class, which has a data member of "Geoemtry" object "boost::shared_ptr geom_;". I got confused here about where to store the actual "State" (i.e. background) SST variable.

My "Geometry" class has usual members of
"std::unique_ptratlas::functionspace::StructuredColumns atlasFunctionSpace_;" and "std::unique_ptratlas::FieldSet atlasFieldSet_;"
My understanding is that, atlasFieldSet_ is where we actually store the "State" or Field (in my case SST) variables, these fields are created by "atlasFunctionSpace_->createField(...)" function. So when I try to implement those member functions of class "State", e.g.
"void zero();"
"double norm() const;"
"void accumul(const double &, const State &);"
"void read(const eckit::Configuration &);"
I actually need to operate on the "geom_->atlasFieldSet_". Am I right? What is the standard structure to do this? I could not find much helpful examples or documents.

What's confusing me is that for the models "oops/qg" and "oops/l95", whose "State" classes both have a member that is similar to "FieldSet" (e.g. in qg model, it has "std::unique_ptr fields_;"), in which it contains an object of "Geometry" class. Should I also implement such FieldsQG class?

@travissluka , Please take a look. Thanks,

we have SST observation files already for use with SOCA tests.
I will copy these into the oisst project

A recent oops update (https://github.com/JCSDA/oops/pull/932) has added new interfaces required for increment

@loganchen39 OISST will not compile currently, I'll issue a PR shortly that fixes this

• fill in implementation for State::print()

First real application!
hofx_nomodel.x is responsible for applying the observation operator to a background state given a set of observation locations.

• requires sample observations (#16 )
• requires GetValues implemented (#18 )

I'll starting working on the Covariance class.

implement the staticbinit.x application. This will use BUMP to calculate the files needed for a basic horizontal correlation. operator

• enable a ctest for staticbinit
  • instead of umdsstCovar as the listed covariance in the yaml file, use BUMP
  • specify a constant horizontal correlation length (let's say, 300km ?)
• change the dirac ctest to instead uses the BUMP static B generated from the above mentioned staticbinit ( change covariance model: BUMP)

at this point, the above dirac ctest should produce output fields with increments that are spread out horizontally

#51 added code to calculate the gridcell area, but this only works correctly for a regular, global, latlon grid. This should be generalized to support any arbitrary grid that ATLAS generates.

.... not a high priority, here simply as a reminder if we think we might need it in the future.

The Field::read() and Field::write() are converting to/from Kelvin so that the value inside jedi is in Celsius, however this means the saved/read Increment files are wrong (as exposed by PR #43). This isn't a huge problem now, but will make interpreting dirac output awkward in the future.

The easiest solution is probably to add a kelvin: true line to the appropriate state I/O sections of the yaml files and check inside Field::read() and Field::write() if that parameter exists, and only do the conversion if that parameter exists and is true.

SST observations within JEDI are currently in Celsius, but our model state here is in Kelvin. We need to either

1. convert the state between K and C when reading/writing the state or
2. Convert between K and C within GetValues.

1 might be more convenient, especially if we plan on using the CoolSkin observation operator, which I think currently expects the get the GeoVaLs as Celsius.

Description

Obs formats are changing. Prepare this repo to run with the IODAv2 branches.

• https://github.com/JCSDA-internal/ufo/pull/1135
• https://github.com/JCSDA-internal/ioda/pull/246
• JCSDA-internal/ioda-converters#534

• Replace landmask file with a better one (mainly fixing Central America).
• make sure there are observations in the dirac test that are near Panama so that the network check on the land mask can be properly tested.

I realize the input file is probably using sst, but once inside JEDI we try to use the CF naming conventions whenever possible (http://cfconventions.org/Data/cf-standard-names/76/build/cf-standard-name-table.html). I also realize it doesn't make much of a difference now, but might be helpful further along when hooking up to UFO

Description

Make changes needed to deal with https://github.com/JCSDA-internal/oops/pull/1199

Interfaces with JEDI typically all assume double instead of float. I've been warned that there might be problems if we continue using float. (as you've seen the tolerances need to be increased quite a bit when using float).
Some people on the JEDI team are looking into whether parts of the interface can use float, but for safety reasons now we should use double.

@travissluka Can we follow the suggestion in post "ecbuild 3.5/3.6 migration" from https://github.com/orgs/JCSDA-internal/teams/jedi , to include fckit, eckit, atlas etc. in the containers and HPC modules instead of building them with the bundles? It needs us to modify the file "UMD-SST/bundle/CMakeLists.txt".

I've been having trouble recompiling JEDI/UMD-SST on Cheyenne, mainly due to their migration to ecbuild35. Below is the link error. The reason is previously we used eckit version "jcsda-1.11.6.jcsda2" (under /glade/work/miesch/modules/gnu-9.1.0/openmpi-4.0.3/eckit), now after upgrade all of these oops/atlas/fckit need to use upgraded eckit version "ecmwf-1.16.0".
I tried several times providing the path to only the upgraded eckit lib, which failed. I also tried to compile everything from the beginning, which also failed due to the upgraded ecbuild. I figured it would be more effective if we just follow JEDI's suggestion.

[ 96%] Linking CXX executable ../../../bin/umdsst_convertstate.x
/usr/bin/ld: ../../../lib/liboops.so: undefined reference to eckit::Main::createMetricsLogTarget() const' /usr/bin/ld: ../../../lib/libatlas.so.0.20: undefined reference to eckit::Buffer::Buffer(void const*, unsigned long)'
/usr/bin/ld: ../../../lib/libatlas_f.so.0.20: undefined reference to `eckit_git_sha1'
......

@travissluka The longitude of our observational SST data (as well as the initial background SST) is within [-180, 180], do you have experience to create such RegularLonLatGrid, most likely with an YAML configuration file? The default longitude for Atlas RegularLonLatGrid is within [0, 360]. If you do not, don't bother, I'll test myself, or at least I can convert the longitude from [-180, 180] to [0, 360].

• missing value from netcdf should be read/written
• within the state field, missing values should be stored using the values given in this file https://github.com/JCSDA-internal/oops/blob/develop/src/oops/util/missingValues.h
• math operators should skip over missing values.

for example, the output from the state test looks like this because missing values are not given

Using the basic "identity" umdsstCovar implemented by #34, get the dirac.x application working.

• check Increment::dirac() to make sure it working correctly. It should use the x,y coordinates in the yaml configuration to set that grid box of the increment field to 1.0. Since we have an MPI decomposed grid, you'll need to convert the global x,y to a local x,y and only do something with it if it is indeed on the PE's local section of the grid. (i'm sure atlas provides the infromation somewhere for converting between global/local indexes)
• enable the ctest for dirac.x

at this point dirac test should produce output files that show increments at the exact x,y locations specified in the yaml. The increment will not yet be spread out horizontally.

pass the land mask to the unstructured grid interpolation when setting up the obs interpolation in GetValues

Currently a fixed global correlation length is used by bump. Add the code needed to pass it the cor_rh field needed to use horizontally varying correlation lengths.

We can decide in a subsequent PR how we want to specify those lengths (two options would be 1) use simple latitude based rossby radius formula, or 2) read in a file with predetermined lengths )

This class does the interpolation of State variables to observations locations, and passes it to UFO to run the forward operator.
I need to check and see what the best way to do the interpolation is now (there are generic interpolation routines in OOPS that use SABER/BUMP... I'll try to track these down for you)

TestGetValues should pass afterward.

Similar to GetValues except it operates on Increment and so performs the tangent linear / adjoint for interpolation.
TestLinearGetValues should pass after this

GetValues is probably easier, so that one should be done first.

To reduce redundant code, and improve maintainability in the future, create a base class that State and Increment inherit from, moving the common methods into that class.

@loganchen39 since this is just reorganizing existing code, and not adding anything new, I can do this while you work on implementing the GetValues class

Implement the Increment class in a similar way to what was done with State.

I would recommend first creating a common base class, say Fields that State and Increment will both inherit from. Move the atlasFieldSet_ variable into this class. Then, as you're implementing the members of Increment, if the logic is identical to what is in State, move the entire member function from State to Fields... that way you wont' have to write as much duplicate code in Increment

Currently, the background error variance is effectively 1C. A horizontally varying background error variance needs to be used by providing a LinearVariableChange

• I'm not sure if the saber::StdDevVariableChange can be used for this, probably not, it might only be used for the standard deviation when SABER is used to calculate the full covariance from a sample ensemble.
• otherwise, we need our own variable change class that reads in the background error standard deviation and implements the multiply and multiplyAD methods. Conceptually it shouldn't be very difficult. See https://github.com/JCSDA-internal/soca/blob/develop/src/soca/Transforms/BkgErr/BkgErr.h for an example with SOCA (though we of course we will but using all C++ instead)

For actual standard deviation values to use: we can just create a file with some made up values for now, and then later work on using scientifically meaningful values (perhaps use the data files from the current OISST?)

on creation Geometry should read in a land mask from a file (separate from what is being done in State now), specified by the yaml configuration and save the mask in an atlas field contained in Geometry

BUMP is no longer working for us after a recent update to SABER. It appears that it now expects the area field to be provided (actually it seems that it was expecting area before, but was somehow working even though we weren't providing it.)

To confirm this, I tested by adding the following quick hack to the geometry constructor

atlas::Field area = atlasFunctionSpace_->createField<double>(
                   atlas::option::levels(1) |
                   atlas::option::name("area"));
make_view<double, 2>(area).assign(110000.0 * 110000.0);
atlasFieldSet_->add(area);

And the staticbinit.x is able to run again, though with slightly different answers.

We obviously will want to fill in the area with the actual gridcell area. I'm assuming that ATLAS should have a convenient way of calculating and returning the grid cell area, but I wasn't able to find such a method at first glance, perhaps you'll have better luck @loganchen39 ?

This will have to be fixed before the other pending PR

We need QC filters to properly remove SST observations that are on, or too close, to land (The current method of relying on bad large values will stop working once I add masking to the interpolation)

1. Add sea_area_fraction to GetValues. Similar to how sea_surface_temperature works there, except pull the data from the land mask that was just added to Geometry. This does not need to be added to the LinearGetValues
2. add a QC filter in the yaml for the var:

  filter: Domain Check
  where:
  - variable: {name: sea_area_fraction@GeoVaLs}
    minvalue: 1.0

You should be able to disable the current QC bound check filter, use this land mask check instead, and have to Var produce valid answers still