diff --git a/amr-wind/equation_systems/AdvOp_Godunov.H b/amr-wind/equation_systems/AdvOp_Godunov.H index e0871b33fb..43231eb05b 100644 --- a/amr-wind/equation_systems/AdvOp_Godunov.H +++ b/amr-wind/equation_systems/AdvOp_Godunov.H @@ -168,6 +168,11 @@ struct AdvectionOp< (godunov_scheme == godunov::scheme::WENO_JS)) { m_allow_inflow_on_outflow = true; } + // if state is NPH, then n and n+1 are known, and only + // spatial extrapolation is performed + const amrex::Real dt_extrap = + (fstate == FieldState::NPH) ? 0.0 : dt; + HydroUtils::ComputeFluxesOnBoxFromState( bx, PDE::ndim, mfi, (PDE::multiply_rho ? rhotrac : tra_arr), @@ -177,7 +182,7 @@ struct AdvectionOp< known_edge_state, u_mac(lev).const_array(mfi), v_mac(lev).const_array(mfi), w_mac(lev).const_array(mfi), divu, - src_term(lev).const_array(mfi), geom[lev], dt, + src_term(lev).const_array(mfi), geom[lev], dt_extrap, dof_field.bcrec(), dof_field.bcrec_device().data(), iconserv.data(), godunov_use_ppm, godunov_use_forces_in_trans, is_velocity, diff --git a/amr-wind/equation_systems/icns/icns.H b/amr-wind/equation_systems/icns/icns.H index 9337fa78c4..f6e11b1fd8 100644 --- a/amr-wind/equation_systems/icns/icns.H +++ b/amr-wind/equation_systems/icns/icns.H @@ -42,8 +42,8 @@ struct ICNS : VectorTransport static constexpr bool multiply_rho = true; static constexpr bool has_diffusion = true; - // No n+1/2 state for velocity for now - static constexpr bool need_nph_state = false; + // n+1/2 velocity for advection iterations + static constexpr bool need_nph_state = true; }; } // namespace amr_wind::pde diff --git a/amr-wind/equation_systems/icns/icns_advection.H b/amr-wind/equation_systems/icns/icns_advection.H index 835e46e293..77d61fabab 100644 --- a/amr-wind/equation_systems/icns/icns_advection.H +++ b/amr-wind/equation_systems/icns/icns_advection.H @@ -190,11 +190,15 @@ struct AdvectionOp m_allow_inflow_on_outflow = true; } + // if state is NPH, then n and n+1 are known, and only + // spatial extrapolation is performed + const amrex::Real dt_extrap = + (fstate == FieldState::NPH) ? 0.0 : dt; for (int lev = 0; lev < repo.num_active_levels(); ++lev) { HydroUtils::ExtrapVelToFaces( dof_field(lev), src_term(lev), u_mac(lev), v_mac(lev), w_mac(lev), dof_field.bcrec(), bcrec_device.data(), - repo.mesh().Geom(lev), dt, godunov_use_ppm, + repo.mesh().Geom(lev), dt_extrap, godunov_use_ppm, godunov_use_forces_in_trans, premac_advection_type, limiter_type, m_allow_inflow_on_outflow); } @@ -329,6 +333,10 @@ struct AdvectionOp m_allow_inflow_on_outflow = true; } + // if state is NPH, then n and n+1 are known, and only + // spatial extrapolation is performed + const amrex::Real dt_extrap = + (fstate == FieldState::NPH) ? 0.0 : dt; #ifdef AMREX_USE_OMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif @@ -346,7 +354,7 @@ struct AdvectionOp known_edge_state, u_mac(lev).const_array(mfi), v_mac(lev).const_array(mfi), w_mac(lev).const_array(mfi), divu, fq.const_array(mfi), - geom[lev], dt, dof_field.bcrec(), + geom[lev], dt_extrap, dof_field.bcrec(), dof_field.bcrec_device().data(), iconserv.data(), godunov_use_ppm, godunov_use_forces_in_trans, is_velocity, fluxes_are_area_weighted, diff --git a/amr-wind/incflo.H b/amr-wind/incflo.H index 6ef1fc65d2..dc8174c4b6 100644 --- a/amr-wind/incflo.H +++ b/amr-wind/incflo.H @@ -96,8 +96,8 @@ public: bool regrid_and_update(); void pre_advance_stage1(); void pre_advance_stage2(); - void do_advance(); - void advance(); + void do_advance(const int fixed_point_iteration); + void advance(const int fixed_point_iteration); void prescribe_advance(); void post_advance_work(); @@ -133,7 +133,9 @@ public: void ComputeDt(bool explicit_diffusion); void ComputePrescribeDt(); - void ApplyPredictor(bool incremental_projection = false); + void ApplyPredictor( + const bool incremental_projection = false, + const int fixed_point_iteration = 0); void ApplyCorrector(); void ApplyPrescribeStep(); @@ -172,6 +174,9 @@ private: // Prescribe advection velocity bool m_prescribe_vel = false; + // Fixed point iterations every timestep + int m_fixed_point_iterations{1}; + //! number of cells on all levels including covered cells amrex::Long m_cell_count{-1}; diff --git a/amr-wind/incflo.cpp b/amr-wind/incflo.cpp index f537155c7a..fbcfa4b656 100644 --- a/amr-wind/incflo.cpp +++ b/amr-wind/incflo.cpp @@ -276,7 +276,10 @@ void incflo::Evolve() amrex::Real time1 = amrex::ParallelDescriptor::second(); // Advance to time t + dt - do_advance(); + for (int fixed_point_iteration = 0; + fixed_point_iteration < m_fixed_point_iterations; + ++fixed_point_iteration) + do_advance(fixed_point_iteration); amrex::Print() << std::endl; amrex::Real time2 = amrex::ParallelDescriptor::second(); @@ -309,15 +312,17 @@ void incflo::Evolve() } } -void incflo::do_advance() +void incflo::do_advance(const int fixed_point_iteration) { if (m_sim.has_overset()) { m_ovst_ops.pre_advance_work(); } - if (m_prescribe_vel) { + if (m_prescribe_vel && fixed_point_iteration == 0) { prescribe_advance(); } else { - advance(); + amrex::Print() << "Fixed point iteration " << fixed_point_iteration + << std::endl; + advance(fixed_point_iteration); } if (m_sim.has_overset()) { m_ovst_ops.post_advance_work(); @@ -372,6 +377,11 @@ void incflo::init_physics_and_pde() if (activate_overset) { m_sim.activate_overset(); } + + pp.query("fixed_point_iterations", m_fixed_point_iterations); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + m_fixed_point_iterations > 0, + "The number of fixed point iterations cannot be less than 1"); } auto& pde_mgr = m_sim.pde_manager(); diff --git a/amr-wind/incflo_advance.cpp b/amr-wind/incflo_advance.cpp index 608851f5fe..a78480f1d2 100644 --- a/amr-wind/incflo_advance.cpp +++ b/amr-wind/incflo_advance.cpp @@ -47,15 +47,20 @@ void incflo::pre_advance_stage2() * * \callgraph */ -void incflo::advance() +void incflo::advance(const int fixed_point_iteration) { BL_PROFILE("amr-wind::incflo::advance"); + if (fixed_point_iteration == 0) { + m_sim.pde_manager().advance_states(); + } - m_sim.pde_manager().advance_states(); - - ApplyPredictor(); + ApplyPredictor(false, fixed_point_iteration); if (!m_use_godunov) { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + fixed_point_iteration == 0, + "Fixed point iterations are not supported for MOL"); + ApplyCorrector(); } } @@ -165,10 +170,10 @@ void incflo::advance() *
  • \ref incflo::ApplyProjection "Apply projection" * */ -void incflo::ApplyPredictor(bool incremental_projection) +void incflo::ApplyPredictor( + const bool incremental_projection, const int fixed_point_iteration) { BL_PROFILE("amr-wind::incflo::ApplyPredictor"); - // We use the new time value for things computed on the "*" state Real new_time = m_time.new_time(); @@ -190,6 +195,16 @@ void incflo::ApplyPredictor(bool incremental_projection) const auto& density_old = density_new.state(amr_wind::FieldState::Old); auto& density_nph = density_new.state(amr_wind::FieldState::NPH); + if (fixed_point_iteration > 0) { + icns().fields().field.fillpatch(m_time.current_time()); + // Get n + 1/2 velocity + amr_wind::field_ops::lincomb( + icns().fields().field.state(amr_wind::FieldState::NPH), 0.5, + icns().fields().field.state(amr_wind::FieldState::Old), 0, 0.5, + icns().fields().field, 0, 0, icns().fields().field.num_comp(), + icns().fields().field.num_grow()); + } + // ************************************************************************************* // Compute viscosity / diffusive coefficients // ************************************************************************************* @@ -268,7 +283,10 @@ void incflo::ApplyPredictor(bool incremental_projection) } // Extrapolate and apply MAC projection for advection velocities - icns().pre_advection_actions(amr_wind::FieldState::Old); + const auto fstate_preadv = (fixed_point_iteration == 0) + ? amr_wind::FieldState::Old + : amr_wind::FieldState::NPH; + icns().pre_advection_actions(fstate_preadv); // For scalars only first // ************************************************************************************* @@ -336,11 +354,14 @@ void incflo::ApplyPredictor(bool incremental_projection) } // With scalars computed, compute advection of momentum - icns().compute_advection_term(amr_wind::FieldState::Old); + const auto fstate = (fixed_point_iteration == 0) + ? amr_wind::FieldState::Old + : amr_wind::FieldState::NPH; + icns().compute_advection_term(fstate); // ************************************************************************************* - // Define (or if use_godunov, re-define) the forcing terms and viscous terms - // independently for the right hand side, without 1/rho term + // Define (or if use_godunov, re-define) the forcing terms and viscous + // terms independently for the right hand side, without 1/rho term // ************************************************************************************* icns().compute_source_term(amr_wind::FieldState::NPH); icns().compute_diffusion_term(amr_wind::FieldState::Old); @@ -397,8 +418,25 @@ void incflo::ApplyPredictor(bool incremental_projection) if (m_verbose > 2) { PrintMaxVelLocations("after nodal projection"); } + // ScratchField to store the old np1 + if (fixed_point_iteration > 0 && m_verbose > 0) { + auto vel_np1_old = m_repo.create_scratch_field( + "vel_np1_old", AMREX_SPACEDIM, 1, amr_wind::FieldLoc::CELL); + + auto vel_diff = m_repo.create_scratch_field( + "vel_diff", AMREX_SPACEDIM, 1, amr_wind::FieldLoc::CELL); + // lincomb to get old np1 + amr_wind::field_ops::lincomb( + (*vel_np1_old), -1.0, + icns().fields().field.state(amr_wind::FieldState::Old), 0, 2, + icns().fields().field.state(amr_wind::FieldState::NPH), 0, 0, + icns().fields().field.num_comp(), 1); + + amr_wind::io::print_nonlinear_residual(m_sim, *vel_diff, *vel_np1_old); + vel_np1_old.reset(); + vel_diff.reset(); + } } - // // Apply corrector: // diff --git a/amr-wind/setup/init.cpp b/amr-wind/setup/init.cpp index b669a8b2af..5ad166e110 100644 --- a/amr-wind/setup/init.cpp +++ b/amr-wind/setup/init.cpp @@ -103,7 +103,7 @@ void incflo::InitialIterations() amrex::Print() << "In initial_iterations: iter = " << iter << "\n"; } - ApplyPredictor(true); + ApplyPredictor(true, 0); { auto& vel = icns().fields().field; diff --git a/amr-wind/utilities/console_io.H b/amr-wind/utilities/console_io.H index 974c92b3dc..6b3cdabdc6 100644 --- a/amr-wind/utilities/console_io.H +++ b/amr-wind/utilities/console_io.H @@ -3,6 +3,7 @@ #include #include "AMReX_MLMG.H" +#include "amr-wind/CFDSim.H" namespace amr_wind::io { @@ -20,6 +21,9 @@ void print_mlmg_info(const std::string& solve_name, const amrex::MLMG& mlmg); void print_tpls(std::ostream& /*out*/); +void print_nonlinear_residual( + const CFDSim& sim, ScratchField& vel_diff, const ScratchField& vel_star); + } // namespace amr_wind::io #endif /* CONSOLE_IO_H */ diff --git a/amr-wind/utilities/console_io.cpp b/amr-wind/utilities/console_io.cpp index bff2a9a267..2c16191c20 100644 --- a/amr-wind/utilities/console_io.cpp +++ b/amr-wind/utilities/console_io.cpp @@ -4,6 +4,7 @@ #include "amr-wind/AMRWindVersion.H" #include "AMReX.H" #include "AMReX_OpenMP.H" +#include "amr-wind/CFDSim.H" #ifdef AMR_WIND_USE_NETCDF #include "netcdf.h" @@ -214,4 +215,60 @@ void print_tpls(std::ostream& out) } } +void print_nonlinear_residual( + const CFDSim& sim, ScratchField& vel_diff, const ScratchField& vel_star) +{ + const int nlevels = sim.repo().num_active_levels(); + const auto& mesh = sim.mesh(); + + const auto& velocity_new = sim.pde_manager().icns().fields().field; + const auto& oset_mask = sim.repo().get_int_field("mask_cell"); + + for (int lev = 0; lev < nlevels; ++lev) { + + amrex::iMultiFab level_mask; + if (lev < nlevels - 1) { + level_mask = makeFineMask( + mesh.boxArray(lev), mesh.DistributionMap(lev), + mesh.boxArray(lev + 1), mesh.refRatio(lev), 1, 0); + } else { + level_mask.define( + mesh.boxArray(lev), mesh.DistributionMap(lev), 1, 0, + amrex::MFInfo()); + level_mask.setVal(1); + } + + const auto& velnew_arr = velocity_new(lev).const_arrays(); + const auto& velstar_arr = vel_star(lev).const_arrays(); + const auto& veldiff_arr = vel_diff(lev).arrays(); + const auto& levelmask_arr = level_mask.const_arrays(); + const auto& osetmask_arr = oset_mask(lev).const_arrays(); + + amrex::ParallelFor( + velocity_new(lev), amrex::IntVect(0), AMREX_SPACEDIM, + [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k, int n) noexcept { + if (osetmask_arr[nbx](i, j, k) == 0) { + veldiff_arr[nbx](i, j, k, n) = 0.; + } else { + veldiff_arr[nbx](i, j, k, n) = + (velnew_arr[nbx](i, j, k, n) - + velstar_arr[nbx](i, j, k, n)) * + levelmask_arr[nbx](i, j, k); + } + }); + } + + amrex::Array rms_vel = {0.0}; + + for (int lev = 0; lev < nlevels; ++lev) { + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + rms_vel[idim] += vel_diff(lev).norm2(idim); + } + } + + amrex::Print() << "Norm of change over fixed point iterations in u: " + << rms_vel[0] << ", v: " << rms_vel[1] + << ", w: " << rms_vel[2] << std::endl; +} + } // namespace amr_wind::io diff --git a/test/test_files/abl_godunov_wenoz_fixedpt/abl_godunov_wenoz_fixedpt.inp b/test/test_files/abl_godunov_wenoz_fixedpt/abl_godunov_wenoz_fixedpt.inp new file mode 100644 index 0000000000..6e8ab161cf --- /dev/null +++ b/test/test_files/abl_godunov_wenoz_fixedpt/abl_godunov_wenoz_fixedpt.inp @@ -0,0 +1,74 @@ +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# SIMULATION STOP # +#.......................................# +#time.stop_time = 22000.0 # Max (simulated) time to evolve +time.stop_time = 50.0 # Max (simulated) time to evolve +time.max_step = 100 # Max number of time steps + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# TIME STEP COMPUTATION # +#.......................................# +time.fixed_dt = 0.5 # Use this constant dt if > 0 +time.cfl = 0.95 # CFL factor + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# INPUT AND OUTPUT # +#.......................................# +time.plot_interval = 10 # Steps between plot files +time.checkpoint_interval = -5 # Steps between checkpoint files + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# PHYSICS # +#.......................................# +incflo.gravity = 0. 0. -9.81 # Gravitational force (3D) +incflo.density = 1.0 # Reference density +incflo.use_godunov = 1 +incflo.godunov_type = "weno_z" +incflo.diffusion_type = 2 +incflo.fixed_point_iterations = 10 +transport.viscosity = 1.0e-5 +transport.laminar_prandtl = 0.7 +transport.turbulent_prandtl = 0.3333 +turbulence.model = Smagorinsky +Smagorinsky_coeffs.Cs = 0.135 + + +incflo.physics = ABL +ICNS.source_terms = BoussinesqBuoyancy CoriolisForcing ABLForcing +BoussinesqBuoyancy.reference_temperature = 300.0 +ABL.reference_temperature = 300.0 +CoriolisForcing.latitude = 41.3 +ABLForcing.abl_forcing_height = 90 + +incflo.velocity = 6.128355544951824 5.142300877492314 0.0 + +ABL.temperature_heights = 650.0 750.0 1000.0 +ABL.temperature_values = 300.0 308.0 308.75 + +ABL.kappa = .41 +ABL.surface_roughness_z0 = 0.15 + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# ADAPTIVE MESH REFINEMENT # +#.......................................# +amr.n_cell = 48 48 48 # Grid cells at coarsest AMRlevel +amr.max_level = 0 # Max AMR level in hierarchy + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# GEOMETRY # +#.......................................# +geometry.prob_lo = 0. 0. 0. # Lo corner coordinates +geometry.prob_hi = 1000. 1000. 1000. # Hi corner coordinates +geometry.is_periodic = 1 1 0 # Periodicity x y z (0/1) + +# Boundary conditions +zlo.type = "wall_model" + +zhi.type = "slip_wall" +zhi.temperature_type = "fixed_gradient" +zhi.temperature = 0.003 # tracer is used to specify potential temperature gradient + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# VERBOSITY # +#.......................................# +incflo.verbose = 0 # incflo_level diff --git a/unit_tests/equation_systems/test_pde.cpp b/unit_tests/equation_systems/test_pde.cpp index af23217c34..907e0ac0b9 100644 --- a/unit_tests/equation_systems/test_pde.cpp +++ b/unit_tests/equation_systems/test_pde.cpp @@ -19,7 +19,7 @@ TEST_F(PDETest, test_pde_create_godunov) EXPECT_EQ(pde_mgr.scalar_eqns().size(), 1); - EXPECT_EQ(mesh().field_repo().num_fields(), 21); + EXPECT_EQ(mesh().field_repo().num_fields(), 22); } TEST_F(PDETest, test_pde_create_mol) @@ -35,7 +35,7 @@ TEST_F(PDETest, test_pde_create_mol) EXPECT_EQ(pde_mgr.scalar_eqns().size(), 1); - EXPECT_EQ(mesh().field_repo().num_fields(), 25); + EXPECT_EQ(mesh().field_repo().num_fields(), 26); } } // namespace amr_wind_tests