diff --git a/dynamics/spam/Dycore.h b/dynamics/spam/Dycore.h index 973175d2..628b4c95 100644 --- a/dynamics/spam/Dycore.h +++ b/dynamics/spam/Dycore.h @@ -243,26 +243,24 @@ class Dycore { // adjust crm_per_phys and dtcrm if needed - if(params.adjust_crm_per_phys_using_vert_cfl) - { - tendencies.adjust_crm_per_phys_using_vert_cfl(params, prognostic_vars); - } - //tendencies.a - // Time stepping loop + if (params.adjust_crm_per_phys_using_vert_cfl) { + tendencies.adjust_crm_per_phys_using_vert_cfl(params, prognostic_vars); + } + // tendencies.a + // Time stepping loop debug_print("start time stepping loop", par.masterproc); for (uint nstep = 0; nstep < params.crm_per_phys; nstep++) { yakl::fence(); time_integrator->step_forward(params.dtcrm); - if(params.clip_vertical_velocities) - { -tendencies.clip_vertical_velocities(params, prognostic_vars); -} + if (params.clip_vertical_velocities) { + tendencies.clip_vertical_velocities(params, prognostic_vars); + } - if (params.clip_negative_densities) { - tendencies.clip_negative_densities(prognostic_vars); - } + if (params.clip_negative_densities) { + tendencies.clip_negative_densities(prognostic_vars); + } #if defined PAMC_AN || defined PAMC_MAN if (params.check_anelastic_constraint) { @@ -299,7 +297,6 @@ tendencies.clip_vertical_velocities(params, prognostic_vars); } prevstep += params.crm_per_phys; - // convert dynamics state to Coupler state tendencies.convert_dynamics_to_coupler_state( coupler, prognostic_vars, constant_vars, params.couple_wind, diff --git a/dynamics/spam/src/models/extrudedmodel.h b/dynamics/spam/src/models/extrudedmodel.h index 0226665e..77f1953d 100644 --- a/dynamics/spam/src/models/extrudedmodel.h +++ b/dynamics/spam/src/models/extrudedmodel.h @@ -396,8 +396,8 @@ class ModelTendencies : public ExtrudedTendencies { #endif } -void clip_vertical_velocities(ModelParameters ¶ms, FieldSet &x) -{ + void clip_vertical_velocities(ModelParameters ¶ms, + FieldSet &x) { real max_w = params.max_w; const auto &primal_topology = primal_geometry.topology; @@ -412,16 +412,19 @@ void clip_vertical_velocities(ModelParameters ¶ms, FieldSet &x) SimpleBounds<4>(primal_topology.nl, primal_topology.n_cells_y, primal_topology.n_cells_x, primal_topology.nens), YAKL_LAMBDA(int k, int j, int i, int n) { - real wedgelength = primal_geometry.get_area_01entity(k + pks, j + pjs, i + pis, n); - wvar(0, k + pks, j + pjs, i + pis, n) = - std::max(-max_w * wedgelength, std::min(max_w * wedgelength, wvar(0, k + pks, j + pjs, i + pis, n))); + real wedgelength = + primal_geometry.get_area_01entity(k + pks, j + pjs, i + pis, n); + wvar(0, k + pks, j + pjs, i + pis, n) = + std::max(-max_w * wedgelength, + std::min(max_w * wedgelength, + wvar(0, k + pks, j + pjs, i + pis, n))); }); } -void adjust_crm_per_phys_using_vert_cfl(ModelParameters ¶ms, FieldSet &x) -{ -//cfl = c * dt / dx -//dt = cfl * dx / c + void adjust_crm_per_phys_using_vert_cfl(ModelParameters ¶ms, + FieldSet &x) { + // cfl = c * dt / dx + // dt = cfl * dx / c const auto &primal_topology = primal_geometry.topology; const auto wvar = x.fields_arr[WVAR].data; @@ -429,7 +432,8 @@ void adjust_crm_per_phys_using_vert_cfl(ModelParameters ¶ms, FieldSet(primal_topology.nl, primal_topology.n_cells_y, primal_topology.n_cells_x, primal_topology.nens), YAKL_LAMBDA(int k, int j, int i, int n) { - real wedgelength = primal_geometry.get_area_01entity(k + pks, j + pjs, i + pis, n); - real wval = std::abs(wvar(0, k + pks, j + pjs, i + pis, n)) / wedgelength; - min_dt_4d(k, j, i, n) = target_cfl * wedgelength / (wval + eps); + real wedgelength = + primal_geometry.get_area_01entity(k + pks, j + pjs, i + pis, n); + real wval = + std::abs(wvar(0, k + pks, j + pjs, i + pis, n)) / wedgelength; + min_dt_4d(k, j, i, n) = target_cfl * wedgelength / (wval + eps); }); - real min_dt = yakl::intrinsics::minval( min_dt_4d ); - real adj_min_dt = std::min(params.dtcrm, min_dt); + real min_dt = yakl::intrinsics::minval(min_dt_4d); + real adj_min_dt = std::min(params.dtcrm, min_dt); - params.crm_per_phys = (int) std::ceil(params.dt_crm_phys / adj_min_dt); - params.dtcrm = params.dt_crm_phys / params.crm_per_phys; + params.crm_per_phys = (int)std::ceil(params.dt_crm_phys / adj_min_dt); + params.dtcrm = params.dt_crm_phys / params.crm_per_phys; - // std::cout << "min_dt: " << min_dt << "\n"; - // std::cout << "adj_min_dt: " << adj_min_dt << "\n"; - // std::cout << "dtcrm: " << params.dtcrm << "\n"; - // std::cout << "dt_crm_phys: " << params.dt_crm_phys << "\n"; - // std::cout << "crm per phys: " << params.crm_per_phys << "\n"; - -} + // std::cout << "min_dt: " << min_dt << "\n"; + // std::cout << "adj_min_dt: " << adj_min_dt << "\n"; + // std::cout << "dtcrm: " << params.dtcrm << "\n"; + // std::cout << "dt_crm_phys: " << params.dt_crm_phys << "\n"; + // std::cout << "crm per phys: " << params.crm_per_phys << "\n"; + } void convert_dynamics_to_coupler_state(PamCoupler &coupler, const FieldSet &prog_vars, diff --git a/dynamics/spam/src/models/layermodel.h b/dynamics/spam/src/models/layermodel.h index 011f5e97..27588660 100644 --- a/dynamics/spam/src/models/layermodel.h +++ b/dynamics/spam/src/models/layermodel.h @@ -129,11 +129,11 @@ class ModelTendencies : public Tendencies { } } -void clip_vertical_velocities(ModelParameters ¶ms, FieldSet &x) -{} + void clip_vertical_velocities(ModelParameters ¶ms, + FieldSet &x) {} -void adjust_crm_per_phys_using_vert_cfl(ModelParameters ¶ms, FieldSet &x) -{} + void adjust_crm_per_phys_using_vert_cfl(ModelParameters ¶ms, + FieldSet &x) {} void compute_constants(FieldSet &const_vars, FieldSet &x) override {}