Skip to content

Commit

Permalink
Merge #111
Browse files Browse the repository at this point in the history 
111: unify timesteps in the model r=yairchn a=yairchn

This PR unifies the time steps in the model across the updrafts and the grid mean the environment and the microphysics. The tilmestep is set to 3 sec which works. Next step would be to add a CFL condition that combines updraft velocity and rainfall velocity to get the optimal tilmestep. 

** behavioral changes !! ** 

take a close look

Co-authored-by: yairchn <[email protected]>
  • Loading branch information
@bors
bors[bot] and yairchn authored Aug 16, 2021
2 parents 5b345ea + ad7b0f5 commit f034f07
Show file tree
Hide file tree
Showing 4 changed files with 61 additions and 83 deletions.
16 changes: 8 additions & 8 deletions integration_tests/GABLS.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,14 +12,14 @@ using .NameList

best_mse = OrderedDict()

best_mse["updraft_thetal"] = 5.0248694008424222e+00
best_mse["v_mean"] = 4.4630163253493214e+00
best_mse["u_mean"] = 9.6412500983148846e+00
best_mse["tke_mean"] = 2.4683562232142533e+00
best_mse["temperature_mean"] = 8.8678308948818863e-06
best_mse["thetal_mean"] = 8.7900622849842488e-06
best_mse["Hvar_mean"] = 1.2891667586703637e+01
best_mse["QTvar_mean"] = 4.4424755799649029e-01
best_mse["updraft_thetal"] = 5.0248696023347037e+00
best_mse["v_mean"] = 4.4593534457868529e+00
best_mse["u_mean"] = 9.6414943665200035e+00
best_mse["tke_mean"] = 2.4674095133951375e+00
best_mse["temperature_mean"] = 8.8584843672667532e-06
best_mse["thetal_mean"] = 8.7856734759460943e-06
best_mse["Hvar_mean"] = 1.2892749042279126e+01
best_mse["QTvar_mean"] = 4.4456710317999498e-01

@testset "GABLS" begin
case_name = "GABLS"
Expand Down
26 changes: 13 additions & 13 deletions integration_tests/TRMM_LBA.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,19 +14,19 @@ using .NameList
CLIMAParameters.Planet.T_freeze(::EarthParameterSet) = 100.0

best_mse = OrderedDict()
best_mse["qt_mean"] = 2.1123818327762196e+00
best_mse["updraft_area"] = 2.2936100671363765e+04
best_mse["updraft_w"] = 9.7871834647382809e+02
best_mse["updraft_qt"] = 3.0602076184487512e+01
best_mse["updraft_thetal"] = 1.1001487658208509e+02
best_mse["v_mean"] = 2.9250821023986430e+02
best_mse["u_mean"] = 1.6873198920581733e+03
best_mse["tke_mean"] = 9.3648065006898219e+02
best_mse["temperature_mean"] = 8.1826613372337135e-04
best_mse["ql_mean"] = 7.2432912110535301e+02
best_mse["thetal_mean"] = 8.2748526945473598e-03
best_mse["Hvar_mean"] = 3.5701761092936244e+03
best_mse["QTvar_mean"] = 1.7865858561907439e+03
best_mse["qt_mean"] = 2.1180570149373197e+00
best_mse["updraft_area"] = 2.2911123097125790e+04
best_mse["updraft_w"] = 9.9122631422628353e+02
best_mse["updraft_qt"] = 3.0750107437154107e+01
best_mse["updraft_thetal"] = 1.1001770046016014e+02
best_mse["v_mean"] = 2.9250578870751696e+02
best_mse["u_mean"] = 1.6873177041648653e+03
best_mse["tke_mean"] = 9.3810901861977175e+02
best_mse["temperature_mean"] = 8.1897112533893871e-04
best_mse["ql_mean"] = 7.3150520783875129e+02
best_mse["thetal_mean"] = 8.2746696205323478e-03
best_mse["Hvar_mean"] = 3.5185010182273427e+03
best_mse["QTvar_mean"] = 1.7745546315637387e+03

@testset "TRMM_LBA" begin
case_name = "TRMM_LBA"
Expand Down
58 changes: 21 additions & 37 deletions src/EDMF_Rain.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -120,51 +120,35 @@ function solve_rain_fall(
term_vel[k] = terminal_velocity(Rain.C_drag, Rain.MP_n_0, QR.values[k], self.Ref.rho0_half[k])
end

# calculate the allowed timestep (CFL_limit >= v dt / dz)
# TODO: report bug: dt_rain has no value in else-case
if !(maximum(term_vel) ≈ 0)
dt_rain = min(dt_model, CFL_limit * grid.dz / maximum(term_vel))
else
dt_rain = dt_model
end

# rain falling through the domain
while t_elapsed < dt_model
# TODO: verify translation
@inbounds for k in reverse(real_center_indicies(grid))
CFL_out = dt_rain / dz * term_vel[k]

if is_toa_center(grid, k)
CFL_in = 0.0
else
CFL_in = dt_rain / dz * term_vel[k + 1]
end

rho_frac = self.Ref.rho0_half[k + 1] / self.Ref.rho0_half[k]
area_frac = 1.0 # RainArea.values[k] / RainArea.new[k]

QR.new[k] = (QR.values[k] * (1 - CFL_out) + QR.values[k + 1] * CFL_in * rho_frac) * area_frac
if QR.new[k] != 0.0
RainArea.new[k] = 1.0
end

term_vel_new[k] = terminal_velocity(Rain.C_drag, Rain.MP_n_0, QR.new[k], self.Ref.rho0_half[k])
end
@inbounds for k in reverse(real_center_indicies(grid))
CFL_out = dt_model / dz * term_vel[k]

t_elapsed += dt_rain
if is_toa_center(grid, k)
CFL_in = 0.0
else
CFL_in = dt_model / dz * term_vel[k + 1]
end

QR.values .= QR.new
RainArea.values .= RainArea.new
if max(CFL_in, CFL_out) > CFL_limit
error("Time step is too large for rain fall velocity!")
end

term_vel .= term_vel_new
rho_frac = self.Ref.rho0_half[k + 1] / self.Ref.rho0_half[k]
area_frac = 1.0 # RainArea.values[k] / RainArea.new[k]

if maximum(abs.(term_vel)) > eps(Float64)
dt_rain = min(dt_model - t_elapsed, CFL_limit * grid.dz / maximum(term_vel))
else
dt_rain = dt_model - t_elapsed
QR.new[k] = (QR.values[k] * (1 - CFL_out) + QR.values[k + 1] * CFL_in * rho_frac) * area_frac
if QR.new[k] != 0.0
RainArea.new[k] = 1.0
end

term_vel_new[k] = terminal_velocity(Rain.C_drag, Rain.MP_n_0, QR.new[k], self.Ref.rho0_half[k])
end

QR.values .= QR.new
RainArea.values .= RainArea.new

term_vel .= term_vel_new
return
end

Expand Down
44 changes: 19 additions & 25 deletions src/Turbulence_PrognosticTKE.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -316,29 +316,24 @@ function compute_prognostic_updrafts(
set_old_with_values(self.UpdVar)

set_updraft_surface_bc(self, GMV, Case)
self.dt_upd = min(TS.dt, 0.5 * get_grid(self).dz / max(maximum(self.UpdVar.W.values), 1e-10))

clear_precip_sources(self.UpdThermo)

while time_elapsed < TS.dt
compute_updraft_closures(self, GMV, Case)
solve_updraft_velocity_area(self)
solve_updraft_scalars(self, GMV)
microphysics(self.UpdThermo, self.UpdVar, self.Rain, TS.dt) # causes division error in dry bubble first time step

set_values_with_new(self.UpdVar)
zero_area_fraction_cleanup(self, GMV)
time_elapsed += self.dt_upd
self.dt_upd = min(TS.dt - time_elapsed, 0.5 * get_grid(self).dz / max(maximum(self.UpdVar.W.values), 1e-10))
# (####)
# TODO - see comment (###)
# It would be better to have a simple linear rule for updating environment here
# instead of calling EnvThermo saturation adjustment scheme for every updraft.
decompose_environment(self, GMV, "values")
saturation_adjustment(self.EnvThermo, self.EnvVar)
buoyancy(self.UpdThermo, self.UpdVar, self.EnvVar, GMV, self.extrapolate_buoyancy)
set_subdomain_bcs(self)
end
compute_updraft_closures(self, GMV, Case)
solve_updraft_velocity_area(self, TS)
solve_updraft_scalars(self, GMV, TS)
microphysics(self.UpdThermo, self.UpdVar, self.Rain, TS.dt) # causes division error in dry bubble first time step

set_values_with_new(self.UpdVar)
zero_area_fraction_cleanup(self, GMV)
# (####)
# TODO - see comment (###)
# It would be better to have a simple linear rule for updating environment here
# instead of calling EnvThermo saturation adjustment scheme for every updraft.
decompose_environment(self, GMV, "values")
saturation_adjustment(self.EnvThermo, self.EnvVar)
buoyancy(self.UpdThermo, self.UpdVar, self.EnvVar, GMV, self.extrapolate_buoyancy)
set_subdomain_bcs(self)

update_total_precip_sources(self.UpdThermo)
return
Expand Down Expand Up @@ -1058,13 +1053,13 @@ function set_subdomain_bcs(self::EDMF_PrognosticTKE)
return
end

function solve_updraft_velocity_area(self::EDMF_PrognosticTKE)
function solve_updraft_velocity_area(self::EDMF_PrognosticTKE, TS::TimeStepping)
grid = get_grid(self)
ref_state = reference_state(self)
kc_surf = kc_surface(grid)
kf_surf = kf_surface(grid)
dzi = grid.dzi
dti_ = 1.0 / self.dt_upd
dti_ = 1.0 / TS.dt
dt_ = 1.0 / dti_

@inbounds for i in xrange(self.n_updrafts)
Expand Down Expand Up @@ -1161,11 +1156,11 @@ function solve_updraft_velocity_area(self::EDMF_PrognosticTKE)
return
end

function solve_updraft_scalars(self::EDMF_PrognosticTKE, GMV::GridMeanVariables)
function solve_updraft_scalars(self::EDMF_PrognosticTKE, GMV::GridMeanVariables, TS::TimeStepping)
grid = get_grid(self)
ref_state = reference_state(self)
dzi = grid.dzi
dti_ = 1.0 / self.dt_upd
dti_ = 1.0 / TS.dt
sa = eos_struct()

@inbounds for i in xrange(self.n_updrafts)
Expand Down Expand Up @@ -2023,7 +2018,6 @@ function compute_covariance_rain(self::EDMF_PrognosticTKE, TS::TimeStepping, GMV
return
end


function compute_covariance_dissipation(self::EDMF_PrognosticTKE, Covar::EnvironmentVariable_2m)
grid = get_grid(self)
ae = 1 .- self.UpdVar.Area.bulkvalues
Expand Down

0 comments on commit f034f07

Please sign in to comment.