Skip to content

User guide

Jump to bottom
Jorn Bruggeman edited this page Mar 18, 2025 · 32 revisions

For now, this page lists answers to common questions. It will evolve into a proper user guide over time.

Input

Dependencies

Biogeochemical models often need additional inputs that fabmos does not provide. For instance, variables describing physical, chemical or biological properties of the environment. This shows up as error messages about unfulfilled dependencies when you execute the run script. For instance:

INFO:FABM:UNFULFILLED DEPENDENCY: absorption_of_silt
INFO:FABM:  This is an interior field.
INFO:FABM:  It is needed by the following model instances:
INFO:FABM:    light
INFO:FABM:UNFULFILLED DEPENDENCY: bottom_stress
INFO:FABM:  This is a horizontal field.
INFO:FABM:  It has units Pa
INFO:FABM:  It is needed by the following model instances:
INFO:FABM:    R4
INFO:FABM:    R6
INFO:FABM:    R8
INFO:FABM:    P1
INFO:FABM:    P2
INFO:FABM:    P3
INFO:FABM:    P4
INFO:FABM:    L2
INFO:FABM:    erosion
INFO:FABM:UNFULFILLED DEPENDENCY: mole_fraction_of_carbon_dioxide_in_air
INFO:FABM:  This is a horizontal field.
INFO:FABM:  It has units 1e-6
INFO:FABM:  It is needed by the following model instances:
INFO:FABM:    O3

Such variables must be provided manually in the run script. Typically, this is done by retrieving the dependency with sim.fabm.get_dependency, and then assigning it a value. For instance:

sim.fabm.get_dependency("mole_fraction_of_carbon_dioxide_in_air").set(280.0)
sim.fabm.get_dependency("absorption_of_silt").set(0.02)
sim.fabm.get_dependency("bottom_stress").set(0.0)

Here, all three dependencies are set to a value that is constant in both time and space. This is not a requirement, however. You can provide time- and/or space dependent values, for instance, by reading values from a NetCDF file like this:

sim.fabm.get_dependency("absorption_of_silt").set(
  fabmos.input.from_nc("<NETCDF_FILE>", "<NETCDF_VARNAME>")
)

The result of fabmos.input.from_nc is a xarray.DataArray that supports lazy arithmetic operations. That allows you to perform simple transformations of NetCDF values. For instance, to scale with a factor 1000 and then add an offset of 10, you could use:

sim.fabm.get_dependency("absorption_of_silt").set(
  fabmos.input.from_nc("<NETCDF_FILE>", "<NETCDF_VARNAME>") * 1000 + 10.0
)

These arithmetic operations are performed lazily: only when and where fabmos actually needs the data. This is very efficient and does not require all data from the NetCDF variable to be read into memory at once.

If you want to provide time-varying climatologies, e.g., representative values for each month, you need to ensure that the input variable describes a single year. You can then use that variable throughout a multi-annual simulation by calling set with the additional argument climatology=True. This works provided the input variable has a valid time coordinate. It is not uncommon to find climatologies in NetCDF files that use a deprecated convention to store the time coordinate, involving reference year 0. This causes error messages such as "Failed to decode variable 'TIMEVAR': unable to decode time units 'hour since 0000-01-01 00:00:00'". To work around this, you need to call fabmos.input.from_nc with extra argument decode_times=False and then provide your own time coordinate by appending something like .assign_coords(TIMEVAR=[cftime.datetime(2000, month, 16) for month in range(1, 13)]) to the result of from_nc.

In general, the set method accepts any xarray.DataArray. For instance, to use a time-varying atmospheric pCO2, you could use:

import pandas as pd
mauna_loa = pd.read_csv(
    "https://gml.noaa.gov/webdata/ccgg/trends/co2/co2_annmean_mlo.txt",
    sep="\s+",
    names=["date", "co2", "unc"],
    index_col="date",
    comment="#",
    converters={"date": lambda y: cftime.datetime(int(y), 7, 1, calendar=calendar)}
)
sim.fabm.get_dependency("mole_fraction_of_carbon_dioxide_in_air").set(mauna_loa.co2.to_xarray())

Note that this downloads pCO2 from the NOAA web server every time you execute the run script. In practice, you will likely want to download this dataset yourself, once, and the replace the URL argument to pd.read_csv by a local file path.

Initial conditions

By default, state variables in the biogeochemical model are initialized to the space-independent value set in your FABM configuration (typically, fabm.yaml). In many cases, you will want to replace this with more realistic, space-dependent values. This can be done in a way that is very similar to the treatment of dependencies (previous section). For instance, to read initial values from a NetCDF file:

sim["VARIABLE_NAME"].set(fabmos.input.from_nc("<NETCDF_FILE>", "<NETCDF_VARNAME>"))

It is worth noting that fabmos has functionality built in to access some datasets that commonly used for model initialization:

  1. World Ocean Atlas (WOA). To download this, use python -m pygetm.input.woa <VARID>, where <VARID> is one of t (temperature), s (salinity), n (nitrate), p (phosphate), i (silicate), o (oxygen). This will download the WOA data for the chosen variable and save it to a new NetCDF file named <VARID>.nc. To see all supported options, use python -m pygetm.input.woa -h.

    You can use the downloaded dataset as usual, for instance:

    sim["N1_p"].set(fabmos.input.from_nc("p.nc", "p"))

  2. Global Ocean Data Analysis Project (GLODAP). To download the this, use python -m pygetm.input.glodap. This will download the mapped GLODAP data product and save it in a new NetCDF file named glodap.nc. To see all supported options, use python -m pygetm.input.glodap -h.

    You can use the downloaded dataset as usual, for instance:

    sim["N1_p"].set(fabmos.input.from_nc("glodap.nc", "PO4"))

    Note that GLODAP provides nutrients as amount per seawater mass (typically, micromol kg-1), whereas FABM-based biogeochemical models typically need amount per volume (e.g., mmol m-3). The former can be converted into the latter by multiplying with seawater density. To make this conversion easier, the pygetm.input.glodap download tool calculates density in kg m-3 and adds it to the NetCDF file. You can use the combine the original nutrient and density values like this:

    density = fabmos.input.from_nc("glodap.nc", "density")
sim["N1_p"].set(fabmos.input.from_nc("glodap.nc", "PO4") * density)

    Note that the final unit of the provided phosphate values here is micromol kg-1 * kg m-3 = micromol m-3. This will not be appropriate for all biogeochemical models. For instance, a model that needs mmol m-3 as units for phosphate would use:

    sim["N1_p"].set(fabmos.input.from_nc("glodap.nc", "PO4") * density * 0.001)

    and a model that needs mol L-1 would use:

    sim["N1_p"].set(fabmos.input.from_nc("glodap.nc", "PO4") * density * 1e-9)

Output

What variables are available?

The best way to find out which variables can be added to output is to (temporarily) customize your run script to print a list of all available variables to screen. The following line does this for all state and diagnostic variables that your chosen biogeochemical model makes available by default:

print([v.name for v in sim.fabm.default_outputs])

Or slightly prettier:

for v in sorted(sim.fabm.default_outputs, key=lambda v: v.name):
    print(f"{v.name}: {v.long_name} ({v.units})")

The following shows all available variables including internally used ones that are typically hidden from output:

print(sim.output_manager.fields.keys())

but this typically shows so many variables to become overwhelming.

Either of these lines can be inserted directly after you create your Simulation (typically called sim). You may want to follow your print line with sys.exit(0) to just get the available variables and not start a simulation. You'll probably wan to run in serial for this specific purpose, i.e., with python <RUNSCRIPT> instead of mpiexec -n <NCPUS> python <RUNSCRIPT>.

For completeness: if you provide an unknown variable name to <OUTPUTFILE>.request, it will complain and show all available variables itself as part of the error message. And you can potentially also diagnose which variables are available by running a short period while saving everything, i.e., with <OUTPUTFILE>.request(*sim.fabm.default_outputs)

Restart files

Restart files are used to store the entire model state produced by one simulation and then use that exact same state as starting point for another simulation. Creating a restart file is very similar to creating a regular NetCDF output file:

out = sim.output_manager.add_restart("res.nc")

Here, the restart file is named "res.nc". The line above can be added anywhere after the simulator (sim) is created and before the simulation is started with sim.start(...).

By default, the call to add_restart:

  • adds all variables that are part of the model state
  • configures the restart file to be written at the very end of the simulation

This is sufficient for most purposes. However, you can add any additional variables you like by calling out.request, just as you would for normal output files. You could use this to add key diagnostics that you consider to particularly informative. You can also save the model state multiple times over the course of the simulation by explicitly setting the output interval. For instance:

out = sim.output_manager.add_restart("res.nc",
    interval=1,
    interval_units=fabmos.TimeUnit.YEARS,
)

This allows you to perform a multi-annual simulation while writing the model state at the end of every year. You can then start a new simulation from any of previously simulated years.

To use restart file to initialize a new simulation, call the following just before sim.start:

out = sim.load_restart("res.nc")

By default, this loads the final state stored in the restart file. If the restart file contains multiple time points and you want to start from an earlier one, you do this by explicitly specifying the time you want to load the state for. For instance:

out = sim.load_restart("res.nc", time=cftime.datetime(2005, 1, 1))

Running in parallel

You can speed up simulations by running fabmos on multiple processing cores, either on your local workstation, or on a High Performance Computing (HPC) cluster. That's done in the same way in both cases: you prefix your call to the run script with mpiexec -n <NCORES>, like this:

mpiexec -n <NCORES> python <RUN_SCRIPT>

Here, <NCORES> is the number of CPU cores you want to use.

How many cores should I use?

In general, you might expect more cores to always increase performance. It is therefore tempting to assign the maximum number of cores available. However, adding more cores does not always produce the expected increase in performance. The problem is that different tasks are needed to solve the model equations. Some parallelize well because each CPU core can handle them completely independently. That’s the case for the biogeochemical source terms. Others parallelize poorly because they require global information, causing each core to request data from all others in order to update the grid points it is responsible for. That quickly becomes a bottleneck, as inter-core communication is relatively slow. This applies to the transport equation of the Transport Matrux Method: each CPU core needs the global model state, essentially because horizontal transports allow the water columns on one core to influence those on other cores. Similarly, output to NetCDF also requires gathering data from all cores (not normally a bottleneck – but if you’d request output at very high frequency for many variables, it could be).

As a result, the model (i.e., fabmos) does not scale linearly (“double the cores halves the runtime”). As you add more cores, the time spent on biogeochemical sources will drop quickly, but the time spent on transport not so much. In fact, at some point, the added cost for inter-core communication can actually increase runtime as you add more cores…

You can test this by running your run script like this:

pygetm-test-scaling --plot <RUN_SCRIPT>

That will perform the same simulation with different numbers of cores and plot the run times for each. Before doing this, you’ll likely want to set the simulation duration to something short (1 month – 1 year) – otherwise you’ll wait a long time for this analysis to complete. NB this is not so easy to run on a cluster – we generally only use it locally.

Typically, adding the first few cores decreases runtime noticeably, but after that gains are limited. And for high core numbers, runtime can even increase. For instance, if you allocate all cores on your workstation, background tasks will typically interfere with the simulation. These tasks will compete with fabmos on some cores, which then holds up the entire simulation.

Your mileage will vary – it’ll depend on your system, MPI implementation, biogeochemical model, time steps for transport and biogeochemistry, output frequency, and more… For instance, more complex biogeochemcial models typically make better use of additonal cores, as they spend comparatively more time on source terms, and less on transport.

logo

Clone this wiki locally