Bulk formula implementation

yairchn commented 3 years ago

Description

In 'momentum_bc.jl' the bulk formula is operating on 'state_int⁻' which is the first DG point above z=0. The computed flux from the formula is applied however at the lower level 'state⁻'. This flux is too large for this point and this creating problems that are mostly evident if turbulence is shut off and +1m/s horizontal wind at the surface can turn into -200m/s winds within 30min simulation.
A consistent implementation as suggested by this branch

original bulk formula code

u1⁻ = state_int⁻.ρu / state_int⁻.ρ
u_int⁻_tan = u1⁻ - dot(u1⁻, n) .* SVector(n)
normu_int⁻_tan = norm(u_int⁻_tan)
-- NOTE: difference from design docs since normal points outwards
C = bc_momentum.drag.fn(state⁻, aux⁻, t, normu_int⁻_tan)
τn = C * normu_int⁻_tan * u_int⁻_tan

Consistent implementation (in consistent with the assumptions of bulk formula though)

u1⁻ = state⁻.ρu / state⁻.ρ
u⁻_tan = u1⁻ - dot(u1⁻, n) .* SVector(n)
normu⁻_tan = norm(u⁻_tan)
# NOTE: difference from design docs since normal points outwards
C = bc_momentum.drag.fn(state⁻, aux⁻, t, normu⁻_tan)
τn = C * normu⁻_tan * u⁻_tan
fluxᵀn.ρu += state⁻.ρ * τn
fluxᵀn.energy.ρe += state⁻.ρu' * τn

does ensure the vanishing of the horizontal velocities at the surface but is somewhat inconsistent with the approach of a bulk formula that should operate at some finite height above the surface.

yairchn commented 3 years ago

The issue can be reproduced by calling include("test/Atmos/EDMF/stable_bl_anelastic1d.jl") in a Julia RPEL and after its is done calling in the same RPEL "" include("test/Atmos/EDMF/compute_mse.jl") dons_arr = get_dons_arr(diag_arr) using Plots z = dons_arr[end]["z"] ρu_1 = dons_arr[end]["prog_ρu_1"] plot(ρu_1,z,xlabel = "ρu_1",ylabel = "z (m)") "" The setup is now a non rotating, single column of dry air with some stable temperature profile without any heat fluxes at the surface a sponge at the top and with initial u=1m/s and v=0. In this setup I have no turbulence (i.e. turbulence = ConstantKinematicViscosity(FT(0)),). I would expect in this setup the profiles of u, v, e to stay unchanged in time.

Note - Since it is a 1D Anelastic SingleStack there are no soundwaves and time steps can be large (1sec). Thus I can easily run 0.5h simulation which I am showing here. After 0.5h I see a value of u=-200m/s.

tapios commented 3 years ago

@akshaysridhar suggested raising the lower boundary to a nominal height (e.g., \Delta z = 10 m). Then we can use the corresponding \Delta z for the calculation of fluxes in the similarity theory and apply them as numerical fluxes on the face of the mesh element facing the boundary.

yairchn commented 3 years ago

If I understand correctly this is a change in the manner in which the transfer coefficient 'C' is computed: C = bc_momentum.drag.fn(state⁻, aux⁻, t, normu_int⁻_tan)

but does it entail a change in of how 'τn' is computed as well: τn = C normu_int⁻_tan u_int⁻_tan

namely in this expression, where are 'normu_int⁻_tan' and 'u_int⁻_tan' computed? at the z=0 interior point or at the z>0 interior point?

tapios commented 3 years ago

If I understand correctly this is a change in the manner in which the transfer coefficient 'C' is computed: C = bc_momentum.drag.fn(state⁻, aux⁻, t, normu_int⁻_tan)

but does it entail a change in of how 'τn' is computed as well: τn = C normu_int⁻_tan u_int⁻_tan

namely in this expression, where are 'normu_int⁻_tan' and 'u_int⁻_tan' computed? at the z=0 interior point or at the z>0 interior point?

Both would be computed at the interface point (z=0, but we'd nominally interpret as, e.g., z=10 m).

yairchn commented 3 years ago

Both would be computed at the interface point (z=0, but we'd nominally interpret as, e.g., z=10 m).

sounds like a good plan! I am happy to help if needed (@akshaysridhar).

yairchn commented 3 years ago

@akshaysridhar if you are changing the bulk formula notice that the same issues arises in the energy_bc.jl https://github.com/CliMA/ClimateMachine.jl/blob/6c11bec9c4f8e54a9928973ce821c438efe598c6/src/Atmos/Model/bc_energy.jl#L204 and moisture_bc.jl https://github.com/CliMA/ClimateMachine.jl/blob/6c11bec9c4f8e54a9928973ce821c438efe598c6/src/Atmos/Model/bc_moisture.jl#L115

CliMA / ClimateMachine.jl