diff --git a/amr-wind/core/FieldBCOps.H b/amr-wind/core/FieldBCOps.H index 98b605e4a0..99b6609dbd 100644 --- a/amr-wind/core/FieldBCOps.H +++ b/amr-wind/core/FieldBCOps.H @@ -11,6 +11,18 @@ namespace amr_wind { +namespace { +template +struct FaceOp +{ + AMREX_GPU_DEVICE AMREX_FORCE_INLINE int operator()(const int idir) + { + return (idir == face_dir) ? 1 : 0; + } +}; + +} // namespace + class Field; /** Interface to setup up boundary conditions on a field @@ -42,6 +54,8 @@ public: struct FieldBCNoOp { using FunctorType = FieldBCNoOp; + template + using FunctorFType = FieldBCNoOp; AMREX_GPU_HOST constexpr FieldBCNoOp() = default; @@ -50,6 +64,11 @@ struct FieldBCNoOp constexpr explicit FieldBCNoOp(const Field& /*unused*/) {} FunctorType operator()() const { return *this; } + template + T opf() const + { + return *this; + } AMREX_GPU_DEVICE void operator()( @@ -110,7 +129,7 @@ struct ConstDirichlet /** Sets constant values on specified boundaries * \ingroup field_fillpatch */ -template +template struct DirichletOp { using InflowOpType = InflowOp; @@ -168,10 +187,13 @@ struct DirichletOp continue; } - // Check if the point in question is on the boundary + // Check if the point in question is on the boundary, which + // depends on the location of the field variable + const int bigEnd_for_floc = + domain_box.bigEnd(idir) + FaceOp()(idir); const bool is_bndry = (ori.isLow() ? (iv[idir] < domain_box.smallEnd(idir)) - : (iv[idir] > domain_box.bigEnd(idir))); + : (iv[idir] > bigEnd_for_floc)); if (!is_bndry) { continue; } @@ -204,10 +226,17 @@ struct DirichletOp struct FieldBCDirichlet { using FunctorType = DirichletOp; + template + using FunctorFType = DirichletOp; explicit FieldBCDirichlet(const Field& fld) : m_field(fld) {} inline FunctorType operator()() const { return FunctorType(m_field); } + template + inline T opf() const + { + return T(m_field); + } const Field& m_field; }; @@ -218,6 +247,8 @@ struct BCOpCreator using InflowOpType = typename InflowOp::DeviceType; using WallOpType = typename WallOp::DeviceType; using FunctorType = DirichletOp; + template + using FunctorFType = DirichletOp; explicit BCOpCreator(const Field& fld) : m_field(fld), m_inflow_op(fld), m_wall_op(fld) @@ -235,6 +266,14 @@ struct BCOpCreator m_wall_op.device_instance()}; } + template + inline T opf() const + { + return T{ + m_field, m_inflow_op.device_instance(), + m_wall_op.device_instance()}; + } + const Field& m_field; InflowOp m_inflow_op; WallOp m_wall_op; diff --git a/amr-wind/core/FieldFillPatchOps.H b/amr-wind/core/FieldFillPatchOps.H index c2870d07ff..21421cd071 100644 --- a/amr-wind/core/FieldFillPatchOps.H +++ b/amr-wind/core/FieldFillPatchOps.H @@ -167,6 +167,8 @@ class FieldFillPatchOps : public FieldFillPatchOpsBase { public: using Functor = typename BCOpCreator::FunctorType; + template + using FunctorF = typename BCOpCreator::template FunctorFType; /** * @param field Field whose patches are filled by this instance * @param mesh The mesh instance to determine amrex::Geometry at a level @@ -268,15 +270,33 @@ public: const amrex::Vector& bcrec, const FieldState /*fstate = FieldState::New*/) override { - if (lev == 0) { - amrex::PhysBCFunct> physbc( - m_mesh.Geom(lev), bcrec, bc_functor()); - for (int i = 0; i < static_cast(mfabs.size()); i++) { - amrex::FillPatchSingleLevel( - *mfabs[i], nghost, time, {ffabs[i]}, {time}, 0, 0, 1, - m_mesh.Geom(lev), physbc, i); - } + AMREX_D_TERM( + auto physbcx = std::make_pair( + 0, amrex::PhysBCFunct>>( + m_mesh.Geom(lev), bcrec, + bc_functor_face>())); + , + auto physbcy = std::make_pair( + 1, amrex::PhysBCFunct>>( + m_mesh.Geom(lev), bcrec, + bc_functor_face>())); + , + auto physbcz = std::make_pair( + 2, amrex::PhysBCFunct>>( + m_mesh.Geom(lev), bcrec, + bc_functor_face>()));); + auto physbcs = + std::make_tuple(AMREX_D_DECL(physbcx, physbcy, physbcz)); + std::apply( + [=](auto&&... args) { + ((amrex::FillPatchSingleLevel( + *mfabs[args.first], nghost, time, {ffabs[args.first]}, + {time}, 0, 0, 1, m_mesh.Geom(lev), args.second, + args.first)), + ...); + }, + physbcs); } else { amrex::PhysBCFunct> cphysbc( m_mesh.Geom(lev - 1), bcrec, bc_functor()); @@ -386,6 +406,11 @@ public: protected: Functor bc_functor() { return m_op(); } + template + T bc_functor_face() + { + return m_op.template opf(); + } void check_face_mapper() { diff --git a/amr-wind/equation_systems/icns/icns_advection.H b/amr-wind/equation_systems/icns/icns_advection.H index be5218929c..2351336730 100644 --- a/amr-wind/equation_systems/icns/icns_advection.H +++ b/amr-wind/equation_systems/icns/icns_advection.H @@ -301,7 +301,7 @@ struct AdvectionOp // MAC projection m_macproj_op(fstate, dt); - // Fill mac velocities using velocity BCs + std::cout << "// Fill mac velocities using velocity BCs\n"; if (fvm::Godunov::nghost_state > 0) { amrex::Array mac_vel = { AMREX_D_DECL(&u_mac, &v_mac, &w_mac)}; diff --git a/amr-wind/physics/udfs/Rankine.H b/amr-wind/physics/udfs/Rankine.H index 221e092992..f61375e41a 100644 --- a/amr-wind/physics/udfs/Rankine.H +++ b/amr-wind/physics/udfs/Rankine.H @@ -45,6 +45,27 @@ struct Rankine vel_ref[2] + 0.0)}; field(iv[0], iv[1], iv[2], dcomp + comp) = vel[orig_comp + comp]; + + amrex::IntVect ivhh(AMREX_D_DECL(41, 12, 1)); + amrex::IntVect ivhi(AMREX_D_DECL(40, 12, 1)); + amrex::IntVect ivlo(AMREX_D_DECL(0, 12, 1)); + amrex::IntVect ivlm(AMREX_D_DECL(-1, 12, 1)); + if (iv == ivlo) + amrex::Print() + << "Rankine.H op: filling " << vel[orig_comp + comp] + << " at i=0" << std::endl; + if (iv == ivlm) + amrex::Print() + << "Rankine.H op: filling " << vel[orig_comp + comp] + << " at i=-1" << std::endl; + if (iv == ivhi) + amrex::Print() + << "Rankine.H op: filling " << vel[orig_comp + comp] + << " at i=40" << std::endl; + if (iv == ivhh) + amrex::Print() + << "Rankine.H op: filling " << vel[orig_comp + comp] + << " at i=41" << std::endl; } }; using DeviceType = DeviceOp; diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 030c3d2330..c368057145 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -241,6 +241,7 @@ add_test_re(vortex_dipole_wall_collision) add_test_re(burggraf_flow) add_test_re(abl_godunov_rayleigh_damping) add_test_re(rankine) +add_test_re(rankine-sym) if(NOT AMR_WIND_ENABLE_CUDA) add_test_re(ctv_godunov_plm) diff --git a/test/test_files/rankine-sym/rankine-sym.inp b/test/test_files/rankine-sym/rankine-sym.inp new file mode 100644 index 0000000000..60ca17bb27 --- /dev/null +++ b/test/test_files/rankine-sym/rankine-sym.inp @@ -0,0 +1,80 @@ +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# SIMULATION STOP # +#.......................................# +time.stop_time = 450.0 # vort init at -1250; domain length = 2000; to travel 4500 at 10 m/s, need 450 s +time.max_step = 100000 +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# TIME STEP COMPUTATION # +#.......................................# +time.fixed_dt = 1.0 # Use this constant dt if > 0 +time.cfl = 0.95 # CFL factor +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# INPUT AND OUTPUT # +#.......................................# +time.plot_interval = 1 # Steps between plot files +time.checkpoint_interval = -1 # 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.diffusion_type = 2 +incflo.do_initial_proj = false +incflo.initial_iterations = 0 +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 = CoriolisForcing +#CoriolisForcing.east_vector = 1.0 0.0 0.0 +#CoriolisForcing.north_vector = 0.0 1.0 0.0 +#CoriolisForcing.latitude = 90.0 +#CoriolisForcing.rotational_time_period = 125663.706143592 +incflo.velocity = 10.0 0.0 0.0 +ABL.reference_temperature = 300.0 +ABL.temperature_heights = 0.0 2000.0 +ABL.temperature_values = 300.0 300.0 +ABL.perturb_temperature = false +ABL.perturb_velocity = false +#ABL.cutoff_height = 50.0 +#ABL.perturb_ref_height = 50.0 +#ABL.Uperiods = 4.0 +#ABL.Vperiods = 4.0 +#ABL.deltaU = 1.0 +#ABL.deltaV = 1.0 +#ABL.kappa = .41 +#ABL.surface_roughness_z0 = 0.01 + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# ADAPTIVE MESH REFINEMENT # +#.......................................# +amr.n_cell = 40 60 4 # Grid cells at coarsest AMRlevel +amr.max_level = 0 # Max AMR level in hierarchy +amr.blocking_factor = 4 +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# GEOMETRY # +#.......................................# +geometry.prob_lo = 0. -1500. 0. # Lo corner coordinates +geometry.prob_hi = 2000. 1500. 200. # Hi corner coordinates +geometry.is_periodic = 0 0 0 # Periodicity x y z (0/1) + +# Boundary conditions +xlo.type = "mass_inflow" +xlo.density = 1.0 +xlo.velocity.inflow_type = Rankine +xlo.temperature = 300.0 +xhi.type = "mass_inflow" +xhi.density = 1.0 +xhi.velocity.inflow_type = Rankine +xhi.temperature = 300.0 +ylo.type = "slip_wall" +yhi.type = "slip_wall" +zlo.type = "slip_wall" +zhi.type = "slip_wall" +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# VERBOSITY # +#.......................................# +incflo.verbose = 0 # incflo_level