Docs update

.gitignore

@@ -8,11 +8,22 @@
 *~
 TAGS
 
+# Figs
+*.png
+*.jpg
+*.jpeg
+*.svg
+
+# Movies
+*.gif
+*.mp4
+
 # Docs
 docs/build/
 docs/src/generated/
 !docs/src/assets/*.png
 !docs/src/assets/*.svg
+!docs/src/images/*.png
 
 # Data
 *.vtk
@@ -22,15 +33,5 @@ docs/src/generated/
 *.nc
 *.jld2
 
-# Figs
-*.png
-*.jpg
-*.jpeg
-*.svg
-
-# Movies
-*.gif
-*.mp4
-
 # Julia System Images
 *.so

docs/Manifest.toml

@@ -145,7 +145,7 @@ uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
 version = "1.1.0"
 
 [[CouplerMachine]]
-deps = ["CLIMAParameters", "ClimateMachine", "Dates", "JLD2", "LinearAlgebra", "MPI", "Oceananigans", "PrettyTables", "StaticArrays", "Statistics", "Unitful"]
+deps = ["CLIMAParameters", "ClimateMachine", "Dates", "GaussQuadrature", "JLD2", "KernelAbstractions", "LinearAlgebra", "MPI", "Oceananigans", "PrettyTables", "StaticArrays", "Statistics", "Unitful"]
 path = ".."
 uuid = "4ade58fe-a8da-486c-bd89-46df092ec0c7"
 version = "0.1.0"

docs/make.jl

@@ -7,6 +7,7 @@ const OUTPUT_DIR      = joinpath(@__DIR__, "..", "docs/src/generated")
 
 experiments = [
                "DesignTests/simple_2testcomp.jl",
+               "AdvectionDiffusion/run_script_v2.jl"
               ]
 
 for experiment in experiments
@@ -15,14 +16,16 @@ for experiment in experiments
 end
 
 experiment_pages = [
-                    "Simple Two Component Test" => "generated/simple_2testcomp.md",
+                    "Vertical Column Heat Diffusion" => "generated/simple_2testcomp.md",
+                    "Advection-diffusion on a Sphere" => "generated/run_script_v2.md",
                    ]
 
+interface_pages = ["couplerstate.md", "timestepping.md", "coupledmodel.md"]
+
 pages = Any[
     "Home" => "index.md",
     "Examples" => experiment_pages,
-    "Coupler Object" => "couplerstate.md",
-    "Coupled Timestepping" => "timestepping.md",
+    "Coupler Interface" => interface_pages,
 ]
 
 

docs/src/coupledmodel.md

@@ -0,0 +1,9 @@
+# Coupled Model Components
+
+CouplerMachine provides a wrapper for model components so that they may
+connect to the coupler. The [`CplModel`](@ref) struct packages a component
+model with the information needed for the [`CplSolver`](@ref) to evolve it.
+
+```@docs
+    CouplerMachine.CplModel
+```

docs/src/couplerstate.md

@@ -1,9 +1,14 @@
-# Coupler Object
+# Coupler State
 
-The CouplerMachine defines a type ```CplState``` for a _container_ variable that holds information about the field 
-values that are being used to couple between components. Components can use a ```coupler_put!``` operation to 
-export a set of field values to a ```CplState``` variable. A ```coupler_get``` operation is used to retrieve
-a set field values from a ```CplState``` variable.
+The coupler provides a space to store information being passed between coupled model components at their boundaries. During this exchange, the coupler manages ancillary operations such as regridding, unit conversions, filtering, etc.
+
+The CouplerMachine defines a type [`CplState`](@ref) for a _container_ variable
+that holds information about the field boundary values that are being used to
+couple components. Components can use a [`coupler_put!`](@ref) operation to 
+export a set of field values to a `CplState` variable. A [`coupler_get`](@ref)
+operation is used to retrieve a set of field values from a `CplState` variable.
+During this exchange, the coupler manages ancillary operations such as 
+regridding, unit conversions, or filtering.
 
 ## Coupler Object API
 

docs/src/index.md

@@ -1,13 +1,11 @@
 # CouplerMachine.jl
 
-![](favicon-32x32.png) *Coupling CliMA Models*
-
+*Coupling CliMA Models*
 
 ```@meta
 CurrentModule = CouplerMachine
 ```
 
-
 CouplerMachine.jl provides a means to couple climate model components from and within
 [ClimateMachine.jl](https://github.com/CliMA/ClimateMachine.jl) and [Oceananigans.jl](https://github.com/CliMA/Oceananigans.jl). 
 It is designed to provide a flexible way to map boundary fluxes of quantities, like moisture and heat, that leave one component 
@@ -21,14 +19,10 @@ Functionality includes:
 The CouplerMachine supports coupling components that are all within the same process or coupling components (using MPI) that
 are running on different processes.
 
+| ![Coupler Scheme](images/cplsetup.png) |
+|:--:|
+| *CouplerMachine allows for independent development of interchangeable component models.* |
+
 ```@docs
     CouplerMachine
 ```
-
-```@contents
-Pages = [
-    "examples.md",
-    "timestepping.md",
-    "couplerstate.md",
-]
-```

docs/src/timestepping.md

@@ -1,5 +1,14 @@
 # Coupled Timestepping
 
+`CouplerMachine` currently uses a sequential timestepping framework in which one
+component steps forward before passing its updated state to another. An outer
+"coupled" timestep determines when component models synchronize and coordinate
+with the coupler. Within this coupled timestep, components take an integer number
+of substeps, and evolve independently from each other.
+
+`CplSolver` extends the ODE solver API of
+[ClimateMachine.ODESolvers](https://clima.github.io/ClimateMachine.jl/latest/APIs/Numerics/ODESolvers/ODESolvers/).
+
 ```@docs
 CouplerMachine.CplSolver
 ```

experiments/AdvectionDiffusion/run_script_v2.jl

@@ -1,3 +1,54 @@
+# # Advection-diffusion on a Sphere
+# ## Setup
+# - Passive advection of an anomaly using non-divergent horizontal flow on the sphere (currently 300km/s but switching to Williamson stream function)
+# - Homogeneous tracer concentration of the outer shell couples to an inhomogeneous concentration of the inner shell
+
+# ## Equations
+# The RHS evaluation follows:
+# ```math
+# \frac{\partial \theta}{\partial t} = - \nabla \cdot (\vec{u} \theta + \kappa \nabla \theta) \qquad
+# ```
+# where
+#  -  $\theta$ is the tracer (e.g. potential temperature)
+#  -  $\vec{u}(u,0,0)$ is a constant zonal velocity
+#  -  $\kappa$ is the diffusivity tensor.
+# **Boundary Conditions**
+#  - External boundary conditions (Neumann):
+# ```math
+# \vec{n} \cdot \kappa \nabla \theta = 0 \qquad \text{at } z = 4\text{km}, -4\text{km}
+# ```
+# where $\vec{n}$ is the upward pointing unit normal vector at the boundary.
+#  - Primary coupled boundary condition (forced by secondary domain):
+# ```math
+# \vec{n} \cdot \kappa_a \nabla \theta =f_a= \lambda (\theta - \theta_{secondary}) \qquad \text{at atmos } z = 0
+# ```
+# where $\lambda$ is a constant but will be replaced by a roughness coefficient factor in more complex test cases.
+# $f_a$ and $\kappa_a$ are, respectively, the flux and diffusivity of the atmosphere component.
+#  - Secondary coupled boundary condition (accumulated flux by primary domain or prescribed flux):
+# If accumulated,
+# ```math
+# \vec{n} \cdot \kappa_o\nabla \theta = f_o \qquad \text{at ocean } z = 0 \text{ (or pseudo subsurface)}
+# and
+# ```
+# ```math
+# f_o = f_{accum} = \int_{\tau_0}^{\tau} f_a dt
+# ```
+# where $τ-τ_0$ gives the coupling cycle length.
+# $f_o$ and $\kappa_o$ are, respectively, the flux and diffusivity of the atmosphere component, with
+# $f_{accum}$ being the accumulated flux.
+
+# ```
+# ATMOS
+# ------------- θ @z=0 [lowest atmos level]
+#   A
+#  / \          f = n⋅κ_a ∇θ = λ(θ−θ_secondary​)    
+#   |
+# ------------- θ_secondary​ @z=0 [highest ocean level]
+# OCEAN
+# ```
+
+
+# # Import packages
 using ClimateMachine, MPI
 using ClimateMachine.DGMethods.NumericalFluxes
 using ClimateMachine.DGMethods
@@ -39,18 +90,19 @@ const param_set = EarthParameterSet()
 ClimateMachine.init()
 FT = Float64
 
-# Shared functions
+## Shared functions
 include("domains.jl")
 include("abstractions.jl")
 include("callbacks.jl")
 
-# Main balance law and its components
+## Main balance law and its components
 include("CplMainBL.jl")
 
+# # Set simulation parameters
 nstepsA = 10
 nstepsO = 5
 
-#  Background atmos and ocean diffusivities
+##  Background atmos and ocean diffusivities
 const κᵃʰ = FT(1e4) * 0.0
 const κᵃᶻ = FT(1e-1)
 const κᵒʰ = FT(1e3) * 0.0
@@ -61,38 +113,42 @@ const λ_airsea = FT(L_airsea / τ_airsea)
 coupling_lambda() = (λ_airsea)
 const u_max = FT(1e-5) # max. advective velocity in radians
 
+# # Set up coupled model
+# Define component models and initialize the coupler
+
 function main(::Type{FT}) where {FT}
 
-    # Domain
-    # confusing name - better might be to use something like DeepSphericalShellDomain directly?
+    ## Domain - Spherical Shells
+    ## confusing name - better might be to use something like DeepSphericalShellDomain directly?
     ΩO = DeepSphericalShellDomain(radius = FT(planet_radius(param_set)) - FT(4e3), height = FT(4e3))
     ΩA = DeepSphericalShellDomain(radius = FT(planet_radius(param_set)) , height = FT(4e3))
 
-    # Grid
+    ## Grid
     nelem = (;horizontal = 8, vertical = 4)
     polynomialorder = (;horizontal = 5, vertical = 5)
     overintegrationorder = (;horizontal = 1, vertical = 1)
 
     gridA = DiscontinuousSpectralElementGrid(ΩA, nelem, polynomialorder)
     gridO = DiscontinuousSpectralElementGrid(ΩO, nelem, polynomialorder)
 
-    # Numerics-specific options
-    numerics = (NFfirstorder = CentralNumericalFluxFirstOrder(), NFsecondorder = PenaltyNumFluxDiffusive(), overint_params = (overintegrationorder, polynomialorder) ) #, NFsecondorder = CentralNumericalFluxSecondOrder() )#PenaltyNumFluxDiffusive() )#, overint_params = (overintegrationorder, polynomialorder) ) 
+    ## Numerics-specific options
+    numerics = (NFfirstorder = CentralNumericalFluxFirstOrder(), NFsecondorder = PenaltyNumFluxDiffusive(),
+                overint_params = (overintegrationorder, polynomialorder))
 
-    # Timestepping
+    ## Timestepping
     Δt_ = calculate_dt(gridA, wavespeed = u_max*(ΩA.radius), diffusivity = maximum([κᵃʰ, κᵃᶻ]) ) 
 
     t_time, end_time = ( 0  , 20Δt_ )
 
-    # Collect spatial info, timestepping, balance law and DGmodel for the two components
+    ## Collect spatial info, timestepping, balance law and DGmodel for the two components
     boundary_mask( xc, yc, zc ) = @. ( xc^2 + yc^2 + zc^2 )^0.5 ≈ planet_radius(param_set)
-    
-    # 1. Atmos component
+#+   
+# **1.** Atmos component
 
     ## Prop atmos functions to override defaults
     atmos_structure(λ, ϕ, r) = FT(30)#30.0 + 10.0 * cos(ϕ) * sin(5λ)
     atmos_θⁱⁿⁱᵗ(npt, el, x, y, z) = atmos_structure( lon(x,y,z), lat(x,y,z), rad(x,y,z) )                # Set atmosphere initial state function
-    #atmos_θ_shadowflux(θᵃ, θᵒ, npt, el, xc, yc, zc) = FT(0.0)
+    ##atmos_θ_shadowflux(θᵃ, θᵒ, npt, el, xc, yc, zc) = FT(0.0)
     atmos_θ_shadowflux(θᵃ, θᵒ, npt, el, xc, yc, zc) = is_surface(xc,yc,zc) ? (1.0 / τ_airsea) * (θᵃ - θᵒ) : 0.0 # Set atmosphere shadow boundary flux function
     atmos_calc_kappa_diff(_...) = κᵃʰ, κᵃʰ, κᵃᶻ               # Set atmos diffusion coeffs
     atmos_source_θ(θᵃ, npt, el, xc, yc, zc, θᵒ) = FT(0.0)     # Set atmos source!
@@ -129,8 +185,8 @@ function main(::Type{FT}) where {FT}
         timestepper = LSRK54CarpenterKennedy,
         numerics...,
     )
-
-    # 2. Ocean component
+#+
+# **2.** Ocean component
     ## Prop ocean functions to override defaults
     tropical_heating_1(λ, ϕ, r) = 30.0 + 10.0 * cos(ϕ) * sin(5λ)
     tropical_heating_2(λ, ϕ, r) = 30.0 + 10.0 * cos(ϕ) + 1 * sin(5λ) * cos(ϕ)
@@ -165,18 +221,16 @@ function main(::Type{FT}) where {FT}
         timestepper = LSRK54CarpenterKennedy,
         numerics...,
     )
-
-    # Create a Coupler State object for holding import/export fields.
-    # Try using Dict here - not sure if that will be OK with GPU
+#+
+# Create the coupler object for holding import/export fields and performs mappings
+# and instantiate the coupled timestepper:
     coupler = CplState()
     coupler_register!(coupler, :Ocean_SST, deepcopy(mO.state.θ[mO.boundary]), mO.grid, DateTime(0), u"°C")
     coupler_register!(coupler, :Atmos_MeanAirSeaθFlux, deepcopy(mA.state.F_accum[mA.boundary]), mA.grid, DateTime(0), u"°C")
 
-
     # Instantiate a coupled timestepper that steps forward the components and
     # implements mapings between components export bondary states and
     # other components imports.
-
     compA = (pre_step = preatmos, component_model = mA, post_step = postatmos)
     compO = (pre_step = preocean, component_model = mO, post_step = postocean)
     component_list = (atmosphere = compA, ocean = compO)
@@ -186,8 +240,10 @@ function main(::Type{FT}) where {FT}
         coupling_dt = Δt_,
         t0 = 0.0,
     )
-
-    # For now applying callbacks only to atmos.
+
+#+
+# Creat callbacks and bundle the simulation:
+    ## For now applying callbacks only to atmos.
     callbacks = (
         ExponentialFiltering(),
         VTKOutput((
@@ -213,7 +269,7 @@ end
 
 
 function run(cpl_solver, numberofsteps, cbvector)
-    # Run the model
+    ## Run the model
     solve!(
         nothing,
         cpl_solver;
@@ -222,6 +278,8 @@ function run(cpl_solver, numberofsteps, cbvector)
     )
 end
 
+# # Run simulation
+
 simulation = main(Float64);
 nsteps = Int(simulation.simtime[2] / simulation.coupled_odesolver.dt)
 cbvector = create_callbacks(simulation, simulation.odesolver)