Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Simple Slab Land #8

Merged
merged 20 commits into from
Jun 25, 2021
Merged

Simple Slab Land #8

merged 20 commits into from
Jun 25, 2021

Conversation

LenkaNovak
Copy link
Collaborator

@LenkaNovak LenkaNovak commented Jun 7, 2021

Dry heat advection-diffusion equation in an atmospheric single stack coupled to a prognostic surface temperature equation reflecting a slab land. Dry setup. This case is mainly for software design before implementing more complex cases.

Design: https://github.com/CliMA/CouplerMachine/blob/ln/slabland/experiments/SlabLand/README.md
Conservation tests: under the 16 June section

Project.toml Outdated
@@ -13,6 +13,7 @@ KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
Oceananigans = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We should keep our dependencies as trim as possible - so no Plots in the CouplerMachine dependencies (Pkg.rm("Plots") in the CouplerMachine environment). If you add Plots from the default Julia environment (no --project) then it will still be accessible when you want to use it in other projects without adding it to the dependency list:

julia -e 'using Pkg; Pkg.add("Plots")'


(constant conductances and ConstantViscosity( FT(0), FT(0), FT(1e-5), u = 0 )

## **Pipeline** (probably delete)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this pipeline could be worth holding onto for a little while - it's a clear walkthrough of what is happening during the coupling.

@LenkaNovak LenkaNovak marked this pull request as ready for review June 19, 2021 20:02
@LenkaNovak LenkaNovak requested a review from kmdeck June 19, 2021 20:13
- Upward-pointing net long-wave flux
$$R_{LW} = F_{sfc}^u - F_{sfc}^d = \epsilon \sigma \theta_{sfc}^4 - \epsilon F_{a}$$
- Soil storage
$$G = \kappa_s (\theta_{sfc}-\theta_h)/h_s$$
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

\theta_{sfc} -> T_{sfc}

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good catch! :)


3. Flux / state storage and communication
- the coupler stores 2 fields:
- `EnergyFluxAtmos` which collects `F_ρθ_accum / Δt` calculated and accumulated in Atmos (with $Δt$ denoting the period of the coupling cycle), and is called by land for its $T_{sfc}$ equation
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

maybe rename, e.g. BoundaryEnergyFlux?

@@ -0,0 +1,165 @@
# Boundary Conditions (BCs)

This section is an attempt to bridge mathematical considerations of BCs and their implementations in climate models.
Copy link
Contributor

@jb-mackay jb-mackay Jun 21, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this doc is great and very helpful, but may perhaps belong with the RQA work in ClimateMachine. What do you think @LenkaNovak ?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

RQA said they will add it once they've revamped their BCs, so perhaps we can leave it here until then?

@@ -0,0 +1,98 @@
# Land put and get calls
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I like the coupler directory set up and am glad we are separating out the pre- and post-step function definition from the balance laws; so this is great! There will probably be some other coupler-specific configuration/definitions that we'll need to do in the future so it's good to have a separate place for it.

- Lower (coupled) boundary:
- ρ: `Impenetrable()`
- u: `FreeSlip()`
- θ: `CoupledPrimary()`: $\frac{\partial θ}{\partial t} = - ∇ \cdot F_{tot}$ where $F_{tot}$ is the surface flux calculated by `calculate_land_sfc_fluxes(_...)` using Atmos and Land properties (see below)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is the boundary condition $$\partial θ / \partial n = F_{tot}$$? (not the divergence of the flux)

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not quite, the way we set the BCs in our case is ultimately via the normal flux, i.e. fluxᵀn.ρθ += - F_tot / scaling_params (referred to as local_flux in different parts of the code), and then the model calculates the div of that.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ok! that is what I thought. I was just commenting on the way it is written. (the math expression does not reflect the BC, which is a Neumann BC on flux like you wrote - the math instead reflects the entire RHS at the boundary).

5) Coupler converts Atmos field `state.F_ρθ_accum` to the coupler field `EnergyFluxAtmos`
6) Coupler sends its field `EnergyFluxAtmos` and transforms it into Land field `auxiliary.F_ρθ_prescribed`
5) Land performs its own initialization and physics (`G`), which is then added to the `auxiliary.F_ρθ_prescribed` (corresponding to $F_{tot}$ above) in the `source!` function to solve:
$$T_s^{n+1} = T_s^n +\frac{F_{tot}^{n...}}{\rho_s c_s h_s}\Delta t$$
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i also found this pipeline really helpful.

the use of F_tot is inconsistent here - in line 87 it refers to the non G fluxes, here it refers to G + non G fluxes.

Copy link
Contributor

@kmdeck kmdeck Jun 22, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

edited, reread more carefully and understand the ^n notation now!

If one wanted to use a different explicit stepper on the land side (not forward Euler), is that possible? Would the "coupling cycle" dt be e.g. the different dts for the RK stages?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@kmdeck unless I'm mistaken, yes, the choice of the land timestepper will determine how $T_s^{n+1}$ is computed, and this does not have to be forward Euler, of course.

The "coupling cycle" refers to the synchronization timescale of the two solvers. At the beginning/end of the coupling cycle, coupled field information is exchanged. Within that coupling cycle, the individual components can take as many substeps as they need, but they don't exchange information.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

thanks! @jb-mackay. I'll think about it more! May still follow up with ?s, but that sounds good.

@codecov
Copy link

codecov bot commented Jun 22, 2021

Codecov Report

Merging #8 (28f775a) into main (eb2b6ac) will not change coverage.
The diff coverage is n/a.

❗ Current head 28f775a differs from pull request most recent head e82eab9. Consider uploading reports for the commit e82eab9 to get more accurate results
Impacted file tree graph

@@           Coverage Diff           @@
##             main       #8   +/-   ##
=======================================
  Coverage   40.98%   40.98%           
=======================================
  Files           4        4           
  Lines          61       61           
=======================================
  Hits           25       25           
  Misses         36       36           
Impacted Files Coverage Δ
src/CouplerMachine.jl 100.00% <ø> (ø)
src/CplModel.jl 0.00% <ø> (ø)
src/CplSolver.jl 0.00% <ø> (ø)
src/CplState.jl 77.41% <ø> (ø)

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update eb2b6ac...e82eab9. Read the comment docs.

- it is called by Atmos, which is assumed to be the model with the shortest timestep. This means that the fluxes are only calculated once and are consistent with Atmos's time stepping.

2. Flux accumulation
- In order to integrate fluxes in time consistently with the Atmos time stepper, we define a new prognostic variable `F_ρθ_accum` whose source function accumulates the fluxes, so that $F\_ρθ\_accum = \int F_{tot} dt$
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

maybe this language was causing confusion? the source function doesnt accumulate the fluxes, it is the RHS term of the ODE d(F_rho theta_accum)/dt = F_tot

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Umm, it basically calculates d (F_ρθ_accum) += \int F_{tot} dt. I will clarify the doc, but yes, accumulating via the source function is not ideal. In the CM that's the only easy option we have, but we will be giving feedback to the CMCore team on this to make it less confusing in the new version :).

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

hey! not being critical, just more clarifying for my own sake. I still dont understand the confusion re: using the language source - it seems as though this is exactly the same usage (a RHS term for a differential equation) as we always use it, except in this case the ODE is d(F_rho theta_accum)/dt = F_tot.

I can chat with Valeria more though since she was the one who brought it up!

)
init_state_auxiliary!(m, m.physics.orientation, state_auxiliary, geom)

state_auxiliary.T_sfc = 0
Copy link
Contributor

@kmdeck kmdeck Jun 22, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, at the coupled boundary, aux.T_sfc (of Atmos) is set to be equal to state.T_sfc (of Land) at the beginning of each Atmos coupling cycle. This line just initializes the elements of the MPIStateArray that are not at the boundary (so they should never be used). Technically, this line should not be needed but I was occasionally getting Nans if the whole field was not initialized. 🤷


"""
Boundary conditions
"""
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are these BC methods needed yet? Right now the PDE for land is just an ODE...

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Eeep! Well spotted, I've just removed them! :)

@kmdeck
Copy link
Contributor

kmdeck commented Jun 22, 2021

it seems like a lot of the files in this PR are not specific to the slab land experiment, and looks like some are already in the repo (e.g. utilities/sphere_utils, grid stuff, numerics (filters...), even physics). maybe can be reorganized to have a single copy of these files?

xmax = 10000, # horizontal extent of the SingleStack in the x direction [m]
ymax = 500, # horizontal extent of the SingleStack in the y direction [m]
zmax = 10000, # vertical extent of the SingleStack [m]
g = 9.81, # gravitation acceleration [m /s^2]
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

worth using Clima parameters at this stage?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think for this illustrative example it's ok to leave the params list as is, but it's definitely worth doing for the more formal test cases that we will be publishing as part of CliMA. :)

@LenkaNovak
Copy link
Collaborator Author

it seems like a lot of the files in this PR are not specific to the slab land experiment, and looks like some are already in the repo (e.g. utilities/sphere_utils, grid stuff, numerics (filters...), even physics). maybe can be reorganized to have a single copy of these files?

Good point. These are copied from the RQA interface. We are keeping these folders under each example separate, because the RQA interface was evolving as we were writing these different test cases, so we wanted to have a stable backend for each case. Once the Aleph.jl repo/module (where the RQA code will live) is ready though, we will import these files from there.

Copy link
Contributor

@jb-mackay jb-mackay left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good @LenkaNovak. Once you have addressed @kmdeck's comments feel free to merge. As mentioned, we'll have to re-tidy this and some other experiments when the RQA interface if formalized.

Copy link
Contributor

@kmdeck kmdeck left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

looks good to me @LenkaNovak @jb-mackay. nice job!!

@LenkaNovak LenkaNovak merged commit e743c9f into main Jun 25, 2021
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

3 participants