From 287c82149076e0593cca1c2a333b6c990642adb9 Mon Sep 17 00:00:00 2001 From: Brandon Allen Date: Thu, 18 Jun 2020 17:03:58 -0400 Subject: [PATCH] fill out communication functions and add necessary aux vriables Co-Authored-By: Jean-Michel Campin --- .../HydrostaticBoussinesqModel.jl | 2 + src/Ocean/ShallowWater/ShallowWaterModel.jl | 8 +- src/Ocean/SplitExplicit/SplitExplicitModel.jl | 152 +++++++++++++++--- 3 files changed, 138 insertions(+), 24 deletions(-) diff --git a/src/Ocean/HydrostaticBoussinesq/HydrostaticBoussinesqModel.jl b/src/Ocean/HydrostaticBoussinesq/HydrostaticBoussinesqModel.jl index f721db6ab61..ca884961772 100644 --- a/src/Ocean/HydrostaticBoussinesq/HydrostaticBoussinesqModel.jl +++ b/src/Ocean/HydrostaticBoussinesq/HydrostaticBoussinesqModel.jl @@ -165,6 +165,8 @@ function vars_state_auxiliary(m::HBModel, T) w::T # ∫(-∇⋅u) pkin::T # ∫(-αᵀθ) wz0::T # w at z=0 + uᵈ::SVector{2, T} # velocity deviation from vertical mean + ΔGᵘ::SVector{2, T} # vertically averaged tendency end end diff --git a/src/Ocean/ShallowWater/ShallowWaterModel.jl b/src/Ocean/ShallowWater/ShallowWaterModel.jl index 7616c5a34aa..462e81bb790 100644 --- a/src/Ocean/ShallowWater/ShallowWaterModel.jl +++ b/src/Ocean/ShallowWater/ShallowWaterModel.jl @@ -56,8 +56,8 @@ SWModel = ShallowWaterModel function vars_state_conservative(m::SWModel, T) @vars begin - U::SVector{3, T} η::T + U::SVector{3, T} end end @@ -65,6 +65,12 @@ function vars_state_auxiliary(m::SWModel, T) @vars begin f::SVector{3, T} τ::SVector{3, T} # value includes τₒ, g, and ρ + Gᵁ::SVector{3, T} # integral of baroclinic tendency + ηᶜ::T # cummulate η value over fast steps + Uᶜ::SVector{3, T} # cummulate U value over fast steps + ηˢ::T # starting η field value + Uˢ::SVector{3, T} # starting U field value + Δu::SVector{3, T} # reconciliation Δu = 1/H * (Ū - ∫u) end end diff --git a/src/Ocean/SplitExplicit/SplitExplicitModel.jl b/src/Ocean/SplitExplicit/SplitExplicitModel.jl index cbf05ee265e..116d85c423a 100644 --- a/src/Ocean/SplitExplicit/SplitExplicitModel.jl +++ b/src/Ocean/SplitExplicit/SplitExplicitModel.jl @@ -9,28 +9,134 @@ import ..DGMethods: cummulate_fast_solution!, reconcile_from_fast_to_slow! -@inline initialize_states!( - slow::HydrostaticBoussinesqModel, - fast::ShallowWaterModel, - _..., -) = nothing - -@inline tendency_from_slow_to_fast!( - slow::HydrostaticBoussinesqModel, - fast::ShallowWaterModel, - _..., -) = nothing - -@inline cummulate_fast_solution!( - slow::HydrostaticBoussinesqModel, - fast::ShallowWaterModel, - _..., -) = nothing - -@inline reconcile_from_fast_to_slow!( - slow::HydrostaticBoussinesqModel, - fast::ShallowWaterModel, - _..., -) = nothing +@inline function initialize_states!( + baroclinic::HydrostaticBoussinesqModel, + barotropic::ShallowWaterModel, + model_bc, + model_bt, + state_bc, + state_bt, +) + model_bt.state_auxiliary.ηᶜ .= -0 + model_bt.state_auxiliary.Uᶜ .= (@SVector [-0, -0, -0])' + + state_bt.η .= model_bt.state_auxiliary.ηˢ + state_bt.U .= model_bt.state_auxiliary.Uˢ + + model_bc.state_auxiliary.ΔGᵘ .= -0 + + return nothing +end + +@inline function tendency_from_slow_to_fast!( + baroclinic::HydrostaticBoussinesqModel, + barotropic::ShallowWaterModel, + model_bc, + model_bt, + state_bc, + state_bt, + forcing_tendency, +) + # integrate the tendency + tendency_model = model_bc.modeldata.tendency_model + update_auxiliary_state!( + tendency_model, + tendency_model.balance_law, + forcing_tendency, + 0, + ) + + Nq, Nqk, _, _, nelemv, _, nelemh, _ = basic_grid_info(model_bc) + + ## get top value (=integral over full depth) of ∫du + ∫du = tendency_model.state_auxiliary.∫du + boxy_∫du = reshape(∫du, Nq^2, Nqk, 2, nelemv, nelemh) + flat_∫du = @view boxy_∫du[:, end, :, end, :] + + ## copy into Gᵁ of dgFast + model_bt.state_auxiliary.Gᵁ .= reshape(flat_∫du, Nq^2, 2, nelemh) + + ## scale by -1/H and copy back to ΔGu + # note: since tendency_dg.auxstate.∫du is not used after this, could be + # re-used to store a 3-D copy of "-Gu" + ΔGᵘ = model_bc.state_auxiliary.ΔGᵘ + boxy_ΔGᵘ = reshape(ΔGᵘ, Nq^2, Nqk, 2, nelemv, nelemh) + boxy_ΔGᵘ .= -reshape(flat_∫du, Nq^2, 1, 2, 1, nelemh) / baroclinic.problem.H + + return nothing +end + +@inline function cummulate_fast_solution!( + baroclinic::HydrostaticBoussinesqModel, + barotropic::ShallowWaterModel, + model_bt, + state_bt, + fast_time, + fast_dt, + substep, +) + ### for now no averaging, just save the final time step + model_bt.state_auxiliary.ηᶜ .= state_bt.η + model_bt.state_auxiliary.Uᶜ .= state_bt.U + + ### technically unnecessary right now, but for future changes + model_bt.state_auxiliary.ηˢ .= state_bt.η + model_bt.state_auxiliary.Uˢ .= state_bt.U + + return nothing +end + +@inline function reconcile_from_fast_to_slow!( + baroclinic::HydrostaticBoussinesqModel, + barotropic::ShallowWaterModel, + model_bc, + model_bt, + state_bc, + state_bt, +) + Nq, Nqk, _, _, nelemv, _, nelemh, _ = basic_grid_info(model_bc) + + # need to calculate int_u using integral kernels + # u_slow := u_slow + (1/H) * (u_fast - \int_{-H}^{0} u_slow) + + # Compute: \int_{-H}^{0} u_slow) + ### need to make sure this is stored into aux.∫u + + # integrate vertically horizontal velocity + flow_model = model_bc.modeldata.flow_model + update_auxiliary_state!(flow_model, flow_model.balance_law, state_bc, 0) + + ## get top value (=integral over full depth) + ∫u = flow_model.state_auxiliary.∫u + boxy_∫u = reshape(∫u, Nq^2, Nqk, 2, nelemv, nelemh) + flat_∫u = @view boxy_∫u[:, end, :, end, :] + + ## get time weighted averaged out of cumulative arrays + ## no weights right now + state_bt.state_auxiliary.Uᶜ .*= 1 # / fast_time_rec[2] + state_bt.state_auxiliary.ηᶜ .*= 1 # / fast_time_rec[2] + + ### substract ∫u from U and divide by H + + ### Δu is a place holder for 1/H * (Ū - ∫u) + Δu = model_bt.state_auxiliary.Δu + Δu .= 1 / baroclinic.problem.H * (model_bt.state_auxiliary.Uᶜ - flat_∫u) + + ### copy the 2D contribution down the 3D solution + ### need to reshape these things for the broadcast + boxy_u = reshape(state_bc.u, Nq^2, Nqk, 2, nelemv, nelemh) + boxy_Δu = reshape(Δu, Nq^2, 1, 3, 1, nelemh) + ### this works, we tested it + # currently SWModel has a 3-vector for U + # need to shrink to 2-vector at some point + boxy_u .+= boxy_Δu[:, :, 1:2, :, :] + + ### copy 2D eta over to 3D model + boxy_η_3D = reshape(state_bc.η, Nq^2, Nqk, nelemv, nelemh) + boxy_η_2D = reshape(model_bt.state_auxiliary.ηᶜ, Nq^2, 1, 1, nelemh) + boxy_η_3D .= boxy_η_2D + + return nothing +end end