Climatemachine.jl: Boundary condition ideas

[charleskawczynski]

Description

I was sketching out some ideas for changing BCs recently. The basic idea I thought of is:

• Boundary state will stay the same (on the back-end) based on #1727:

for face in faces
    if face == bctag
        boundary_state!(bl, ...)
    end
end

• we add a method, prognostic_state, for users to overload to return each prognostic state:

prognostic_state(::Mass, ::AtmosModel, state) = state.ρ
prognostic_state(::TotalMoisture, ::AtmosModel, state) = state.moisture.ρq_tot

• Then, boundary_state!, can call a new user-facing method set_bc.
• To do this, we iterate over prognostic_vars, grab the prognostic states, and properly assign them based on whether they're Dirichlet/Neumann/Robin:

abstract type BoundaryConditionType end
struct Dirichlet <: BoundaryConditionType end
struct Neumann <: BoundaryConditionType end
struct Robin <: BoundaryConditionType end

function boundary_state!(bl, state⁺, state⁻, ...)
    for prog_var in prognostic_vars(bl)
        prog_state⁻ = prognostic_state(prog_var, bl, state⁻)
        prog_state⁺ = prognostic_state(prog_var, bl, state⁺)
        bctype = bc_type(bl, prog_var, ...) # should be inferrable based on bl and prog_var
        if bctype isa Neumann
            bcval = set_bc(prog_var, bctype, bl, state⁻)
            prog_state⁺ = som_fun(bcval, prog_state⁻)     # I can't ever remember the relation 🤷🏻‍♂️
        elseif bctype isa Dirichlet
            bcval = set_bc(prog_var, bctype, bl, state⁻)
            prog_state⁺ = som_fun(bcval, prog_state⁻)     # I can't ever remember the relation 🤷🏻‍♂️
        end
    end
end

set_bc(::Mass, ::NoFlowThrough, bl, state, n) = 0
set_bc(::Momentum, ::Impenetrable{NoSlip}, bl, state, n) = -state.ρu
set_bc(::Momentum, ::Impenetrable{FreeSlip}, bl, state, n) = - 2 * dot(state.ρu, n) .* SVector(n)

From the user POV, set_bc is now in the spatially continuous form, and no longer sees ⁺/⁻ states.

cc @simonbyrne @blallen @kpamnany @glwagner @thomasgibson

[jkozdon]

for face in faces
    if face == bctag
        boundary_state!(bl, ...)
    end
end

Do we want to assume we are always on a (semi-)structured mesh? This assumes that.

[charleskawczynski]

Ahh, I didn't realize that. I thought that this is already what we're doing. Maybe my pseudocode there isn't accurate.

[jkozdon]

Also, for some bcs to be specified properly it matters which flux you are using. Also for systems, you should really be setting the whole state vector in some way, not just the single components (otherwise you will be setting too many bcs). Maybe this gets captured through you some_fun though?

[charleskawczynski]

I think that's what I had in mind: split the BC definition from managing valid sets of BCs, and how they need to be managed differently between Dirichlet/Neumann in the context of DG.

[glwagner]

Is this issue about the user interface or the way that balance laws are implemented? I can't tell.

The boundary conditions depend on the mesh and domain, as @jkozdon suggests. The idea I had was to use a specification object called BoundaryConditions whose constructors dispatch on the domain. For a RectangularDomain a user might write

u_boundary_conditions = BoundaryConditions(domain, north=DirichletBoundaryCondition(0), bottom=DirichletBoundaryCondition(0))

This calls some function BoundaryConditions(domain::RectangularDomain, ...). Other constructors are needed for other geometries. I don't know how to make the user interface nice for unstructured meshes --- it probably depends on how one generates / loads / specifies an unstructured mesh. That's a much harder problem.

You can also provide NoSlip() as an alias for DirichletBoundaryCondition(0).

These boundary condition specs are then passed to the model. A named tuple is a simple way to organize the boundary condition specification:

model = HydrostaticBoussinesqModel(boundary_conditions=(u=u_bcs, T=T_bcs))

which works because each Model knows the names of its fields.

As @jkozdon mentions, the underlying implementation of a boundary condition depends on the numerics. Thus within a model constructor (at which point things like the numerical flux algorithm has been specified) further processing must occur.

Is it possible to directly specify the flux across the face of a boundary element? This is the most common type of boundary condition for geophysical applications and is seemingly the most natural "boundary condition" in the weak formulation.

[jkozdon]

Is this issue about the user interface or the way that balance laws are implemented? I can't tell.

I think this is a user interface issue. Not a DG implementation issue.

I think that part of my issue is being a numerical analyst an not a physicist. Saying things like Dirichlet for a system is equations has very little meaning to me (since you can easily enforce too many bcs). Similarly, Neumann and Robin don't mean much to me for first order systems of equations.

There are some mechanisms for enforcing the fluxes directly, though once again you have to be careful not to enforce too many boundary conditions.

[charleskawczynski]

I think this is a user interface issue. Not a DG implementation issue.

I agree. I think we should provide users a way to specify BCs as Dirichlet and Neumann, and then we translate these to valid/consistent/well posed BCs in a separate layer. I think, put another way, it's a separation of concerns issue. If I, as a user, want to specify Dirichlet conditions, I should not be concerned with how I do that in the context of DG (as is the case with boundary_state!(...state⁺, state⁻)), that logic (and validation/well-posedness) should get handled in a separate layer.

[glwagner]

At least for the ocean I am pretty sure we will always consider a system of equations in which one boundary condition is specified on each prognostic variable on every boundary.

though once again you have to be careful not to enforce too many boundary conditions.

An example of specifying too many boundary conditions might help.

In a weak formulation, as far as I can tell, boundary conditions boil down to somehow modeling a flux across boundary-adjacent element faces:

• "Dirichlet" means estimate fluxes based on prescribed field _values_ on the other side of the boundary face.
• "Neumann" means estimate the flux based on a prescribed _gradients_ of a field across the boundary face.

It makes sense to me that specifying these depend on the numerical flux scheme. Directly specifying the flux (rather than using a model appropriate for values or gradients) is presumably simpler (but I'm not sure).

From the modeled fluxes we calculate the tendency of the prognostic field to which the flux applies (thus I can't immediately think of a case where too many boundary conditions are applied).

In a finite volume formulation, its rather trivial to implement a prescribed flux since you just add the prescription to the volume average of the tendency in boundary-adjacent cell. Implementing prescribed values or gradients is a bit more complicated, since you have to introduce a model for fluxes, but can also be done.

PS: I personally prefer fewer names of dead people in the code; "value" and "gradient" seem like nice alternatives to "Dirichlet" and "Neumann" (though I'm sure they were great scientists).

[jkozdon]

An example of specifying too many boundary conditions might help.

An advection equation where you try to enforce the value of the solution on all boundaries:

[glwagner]

That's a good example!

Ocean models have second-order diffusion terms essentially 100% of the time (or higher, if its supported...). In Oceananigans (for example), if you specify value/gradient boundary conditions but have no second-order terms, the boundary conditions are not actually respected. This doesn't prevent us from implementing an interface for specifying boundary conditions though... In principle we should throw an error in this invalid (and rare) usage, which is easy to do. (I've included a comment on this case in an open issue we have on this subject.)

Even with no second-order terms, we can specify fluxes, correct?

[jkozdon]

Even with no second-order terms, we can specify fluxes, correct?

Sure. Not sure about the interface for this in the code though

[blallen]

Even with no second-order terms, we can specify fluxes, correct? Sure. Not sure about the interface for this in the code though

I think there might be some confusing terminology here. Greg means can we specify a momentum flux or energy flux in the physical sense without a second order term. We can always do a numerical flux, but can we do something like the energy flux is XXX J/m/s with just first order terms?

[glwagner]

Isn't the numerical flux an approximation of a physical flux? In terms of specifying fluxes across a boundary, presumably the numerical and physical fluxes are the same. This would have to be the case for the volume integrated budgets to be specified correctly. Sorry if I'm confused... :-/

The "prescribed flux" I'm referring to seems to correspond to what Bassi and Rebay refer to as the "weakly prescribed boundary condition"

In the context of Bassi and Rebay's discussion, the only place where there is no ambiguity regarding the "Navier-Stokes flux function" is, in fact, the boundary, when we know what the flux should be. In that case we don't need a model that joins two elements.

In large-scale geophysical models, it is usually the fluxes of quantities themselves that are modeled, rather than their values or gradients at boundaries. A "drag law" is a model for momentum flux.

[christophernhill]

That's a good example!

Ocean models have second-order diffusion terms essentially 100% of the time (or higher, if its supported...). In Oceananigans (for example), if you specify value/gradient boundary conditions but have no second-order terms, the boundary conditions are not actually respected. This doesn't prevent us from implementing an interface for specifying boundary conditions though... In principle we should throw an error in this invalid (and rare) usage, which is easy to do. (I've included a comment on this case in an open issue we have on this subject.)

Even with no second-order terms, we can specify fluxes, correct?

Isn't it legit to impose a boundary and then explore what solutions can work with that. In @jkozdon advection example aren't there some flows for which the bc is fine e.g. a non-divergent vortex in middle of domain? I think we often use models to explore what solutions are consistent with certain bc's. We are OK that the bc's maybe be incompatible with arbitrary solutions, but are interested in whether solutions we see in the real-world maybe all fall in some subset.

Personally I would miss references to Neumann https://en.wikipedia.org/wiki/Carl_Neumann#/media/File:Carl_Gottfried_Neumann.png , but gradient is easier to know what it means.

[jkozdon]

In @jkozdon advection example aren't there some flows for which the bc is fine e.g. a non-divergent vortex in middle of domain?

Yes there certainly are cases, depends on the sign of

[glwagner]

Speaking with @blallen and @christophernhill --- it seems we can write a user interface for specifying value, gradient, and flux boundary conditions already using the HydrostaticBoussinesqModel, without requiring changes to the underlying balance law framework (and thus not requiring any of @jkozdon 's time). So at least the ocean subcomponent does not require a more generic interface.

Note that the ocean's user interface for fluxes would depend on the ocean's specific balance law implementation, which could imply the problem is not simply generalizable to any balance law (in accordance with @jkozdon's arguments). The flux boundary condition user interface is the most important one for us.