From 304eb2a1467330014b435ca9c415226a4236985d Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Mon, 18 May 2026 09:38:00 -0700 Subject: [PATCH 1/2] more precision changes --- src/embedded_boundaries/eb_annulus.cpp | 2 +- src/embedded_boundaries/eb_box.cpp | 18 +++++++++++------- src/embedded_boundaries/eb_cylinder.cpp | 2 +- src/embedded_boundaries/eb_sphere.cpp | 2 +- src/setup/init.cpp | 2 +- src/utilities/io.cpp | 4 ++-- 6 files changed, 17 insertions(+), 13 deletions(-) diff --git a/src/embedded_boundaries/eb_annulus.cpp b/src/embedded_boundaries/eb_annulus.cpp index f8c96a83a..b7eed7152 100644 --- a/src/embedded_boundaries/eb_annulus.cpp +++ b/src/embedded_boundaries/eb_annulus.cpp @@ -39,7 +39,7 @@ void incflo::make_eb_annulus() // Compute distance between cylinder centres Real offset = 0.0; for(int i = 0; i < AMREX_SPACEDIM; i++) - offset += std::pow(outer_center[i] - inner_center[i], 2); + offset += amrex::Math::powi<2>(outer_center[i] - inner_center[i]); offset = std::sqrt(offset); // Check that the inner cylinder is fully contained in the outer one diff --git a/src/embedded_boundaries/eb_box.cpp b/src/embedded_boundaries/eb_box.cpp index 3c361e1d8..9a372ae09 100644 --- a/src/embedded_boundaries/eb_box.cpp +++ b/src/embedded_boundaries/eb_box.cpp @@ -32,7 +32,11 @@ void incflo::make_eb_box() ************************************************************************/ Vector boxLo(AMREX_SPACEDIM), boxHi(AMREX_SPACEDIM); - Real offset = 1.0e-15; +#ifdef AMREX_USE_FLOAT + Real offset = Real(1e-8); +#else + Real offset = Real(1e-15); +#endif bool inside = true; for(int i = 0; i < AMREX_SPACEDIM; i++) @@ -54,8 +58,8 @@ void incflo::make_eb_box() // putting them one domain width away if(geom[0].isPeriodic(0)) { - xlo = 2.0 * geom[0].ProbLo(0) - geom[0].ProbHi(0); - xhi = 2.0 * geom[0].ProbHi(0) - geom[0].ProbLo(0); + xlo = Real(2) * geom[0].ProbLo(0) - geom[0].ProbHi(0); + xhi = Real(2) * geom[0].ProbHi(0) - geom[0].ProbLo(0); } Real ylo = boxLo[1] + offset; @@ -65,8 +69,8 @@ void incflo::make_eb_box() // putting them one domain width away if(geom[0].isPeriodic(1)) { - ylo = 2.0 * geom[0].ProbLo(1) - geom[0].ProbHi(1); - yhi = 2.0 * geom[0].ProbHi(1) - geom[0].ProbLo(1); + ylo = Real(2) * geom[0].ProbLo(1) - geom[0].ProbHi(1); + yhi = Real(2) * geom[0].ProbHi(1) - geom[0].ProbLo(1); } #if (AMREX_SPACEDIM > 2) @@ -77,8 +81,8 @@ void incflo::make_eb_box() // putting them one domain width away if(geom[0].isPeriodic(2)) { - zlo = 2.0 * geom[0].ProbLo(2) - geom[0].ProbHi(2); - zhi = 2.0 * geom[0].ProbHi(2) - geom[0].ProbLo(2); + zlo = Real(2) * geom[0].ProbLo(2) - geom[0].ProbHi(2); + zhi = Real(2) * geom[0].ProbHi(2) - geom[0].ProbLo(2); } #endif diff --git a/src/embedded_boundaries/eb_cylinder.cpp b/src/embedded_boundaries/eb_cylinder.cpp index 02afa9815..cc5ba9c03 100644 --- a/src/embedded_boundaries/eb_cylinder.cpp +++ b/src/embedded_boundaries/eb_cylinder.cpp @@ -33,7 +33,7 @@ void incflo::make_eb_cylinder() pp.getarr("center", centervec, 0, 3); Array center = {AMREX_D_DECL(centervec[0], centervec[1], centervec[2])}; - rotation = (rotation/180.)*M_PI; + rotation = (rotation/180.)*Real(M_PI); // Print info about cylinder amrex::Print() << " " << "\n"; diff --git a/src/embedded_boundaries/eb_sphere.cpp b/src/embedded_boundaries/eb_sphere.cpp index cdf735e06..09ef3563b 100644 --- a/src/embedded_boundaries/eb_sphere.cpp +++ b/src/embedded_boundaries/eb_sphere.cpp @@ -16,7 +16,7 @@ void incflo::make_eb_sphere() { // Initialise sphere parameters bool inside = true; - Real radius = 0.0002; + Real radius = Real(0.0002); Vector centervec(3); // Get sphere information from inputs file. * diff --git a/src/setup/init.cpp b/src/setup/init.cpp index 2a3db0709..db7259667 100644 --- a/src/setup/init.cpp +++ b/src/setup/init.cpp @@ -222,7 +222,7 @@ void incflo::ReadParameters () amrex::Real tol_deg(0.); pp_eb_flow.query("normal_tol", tol_deg); - m_eb_flow.normal_tol = tol_deg*M_PI/amrex::Real(180.); + m_eb_flow.normal_tol = tol_deg*amrex::Real(M_PI)/amrex::Real(180); } if (m_advect_tracer && m_eb_flow.enabled && m_eb_flow.tracer.empty()) { diff --git a/src/utilities/io.cpp b/src/utilities/io.cpp index 801011e10..4f3096bb9 100644 --- a/src/utilities/io.cpp +++ b/src/utilities/io.cpp @@ -175,7 +175,7 @@ void incflo::ReadCheckpointFile() int i = 0; while(lis >> word) { - prob_lo[i++] = std::stod(word); // NOLINT(clang-analyzer-security.ArrayBound) + prob_lo[i++] = amrex::Real(std::stod(word)); // NOLINT(clang-analyzer-security.ArrayBound) } } @@ -186,7 +186,7 @@ void incflo::ReadCheckpointFile() int i = 0; while(lis >> word) { - prob_hi[i++] = std::stod(word); // NOLINT(clang-analyzer-security.ArrayBound) + prob_hi[i++] = amrex::Real(std::stod(word)); // NOLINT(clang-analyzer-security.ArrayBound) } } From 25d4ec75f6359051d8b76024bc0600cd8a8316bb Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Mon, 21 Sep 2026 11:49:56 -0700 Subject: [PATCH 2/2] Fix issues #233-#239 #233 compute_MAC_projected_velocities: the direction_dependent inflow tests read a zeroed scratch MultiFab, whose ghost cells hold the interior copy (0) on the outflow part of an in/out face. The non-strict >= / <= tests therefore treated that as zero inflow and overwrote the extrapolated outflow face velocity with 0, after which Hydro aborts with "no outflow from the direction dependent boundaries". Make the tests strict. #234 IncfloVelFill: the Burggraf (probtype 16) lid profile 16 x^2 (1-x)^2 is the tangential velocity u on the high-y wall, not the normal velocity v. The PR #122 refactor assigned it to norm_vel, so benchmark.burggraf ran a uniform-speed lid with a spurious normal velocity through it. Restore the pre-#122 behaviour. #235 init_plane_poiseuille: probtype 43 shared probtype 42's zero initial velocity, so the outflow part of its direction_dependent x faces started at 0 and InitialProjection's in/out solvability check saw no outflux. Give probtype 43 the zero-net-flux profile u = 6y(1-y)-1 that IncfloVelFill imposes on those faces. #236 ComputeDt: the explicit-diffusion limit used the constant m_mu, ignoring the strain-rate dependent viscosity and the tracer/temperature diffusivities, which are advanced explicitly by the same m_diff_type switch. Build the effective diffusivity per cell from compute_viscosity_at_level, max_n mu_s[n] and mu_T/cp, divide by rho, and take its max. #237 particleData.Redistribute() inside RemakeLevel / MakeNewLevelFromCoarse runs before AmrCore::regrid installs the new BoxArrays, DistributionMappings and finest_level, so it redistributes against the old layout and the first step after a regrid indexes u_mac with stale grid indices. Redistribute once after regrid() returns instead. #238 ApplyPredictor advected the tracer particles during InitialIterations. The fields are restored after each initial iteration but the particle positions are not, so the particles started the run m_initial_iterations*dt ahead of the fluid. Skip the advection when incremental_projection is set. #239 IncfloDenFill / IncfloTracFill / IncfloTempFill never wrote the outflow part of a direction_dependent face (filcc has no case for that BC type and the functors were inflow-only), and decided inflow/outflow from the raw input velocity rather than the probtype profile that IncfloVelFill uses. Factor the prescribed normal velocity into incflo_bc_normal_velocity(), use it in all four functors, and give the scalar functors the same "outflow: copy the first interior cell" branch IncfloVelFill has. Co-Authored-By: Claude Opus 5 (1M context) --- ...ncflo_compute_MAC_projected_velocities.cpp | 18 +- src/incflo.cpp | 7 + src/incflo_apply_predictor.cpp | 5 +- src/incflo_compute_dt.cpp | 76 +++-- src/incflo_regrid.cpp | 8 - src/prob/prob_bc.H | 314 ++++++++++++------ src/prob/prob_init_fluid.cpp | 21 +- 7 files changed, 305 insertions(+), 144 deletions(-) diff --git a/src/convection/incflo_compute_MAC_projected_velocities.cpp b/src/convection/incflo_compute_MAC_projected_velocities.cpp index a40bb292b..d05f68dab 100644 --- a/src/convection/incflo_compute_MAC_projected_velocities.cpp +++ b/src/convection/incflo_compute_MAC_projected_velocities.cpp @@ -208,6 +208,12 @@ incflo::compute_MAC_projected_velocities ( velBC_MF = make_BC_MF(lev, m_bcrec_velocity_d, "velocity"); } + // On a direction_dependent face IncfloVelFill writes the prescribed velocity + // into the ghost cell only where it points into the domain; on the outflow + // part it copies the interior cell, which is 0 here because + // time_dep_inflow_vel was zeroed above. The inflow tests below must + // therefore be strict, otherwise the outflow part of the face is treated as + // zero inflow and the extrapolated outflow velocity is lost. for (MFIter mfi(time_dep_inflow_vel,false); mfi.isValid(); ++mfi) { Box const& bx = mfi.validbox(); @@ -222,11 +228,11 @@ incflo::compute_MAC_projected_velocities ( int n = 0; const auto bc = HydroBC::getBC(i, j, k, n, domain, bc_vel_d, velbc_arr); if (i == dlo.x && ( bc.lo(0) == BCType::ext_dir || - (bc.lo(0) == BCType::direction_dependent && cc_arr(i-1,j,k,0) >= Real(0.0)) ) ) { + (bc.lo(0) == BCType::direction_dependent && cc_arr(i-1,j,k,0) > Real(0.0)) ) ) { umac_arr(i,j,k) = cc_arr(i-1,j,k,0); } if (i == dhi.x && ( bc.hi(0) == BCType::ext_dir || - (bc.hi(0) == BCType::direction_dependent && cc_arr(i+1,j,k,0) <= Real(0.0)) ) ) { + (bc.hi(0) == BCType::direction_dependent && cc_arr(i+1,j,k,0) < Real(0.0)) ) ) { umac_arr(i+1,j,k) = cc_arr(i+1,j,k,0); } }); @@ -236,11 +242,11 @@ incflo::compute_MAC_projected_velocities ( int n = 1; const auto bc = HydroBC::getBC(i, j, k, n, domain, bc_vel_d, velbc_arr); if (j == dlo.y && ( bc.lo(1) == BCType::ext_dir|| - (bc.lo(1) == BCType::direction_dependent && cc_arr(i,j-1,k,1) >= Real(0.0)) ) ) { + (bc.lo(1) == BCType::direction_dependent && cc_arr(i,j-1,k,1) > Real(0.0)) ) ) { vmac_arr(i,j,k) = cc_arr(i,j-1,k,1); } if (j == dhi.y && ( bc.hi(1) == BCType::ext_dir || - (bc.hi(1) == BCType::direction_dependent && cc_arr(i,j+1,k,1) <= Real(0.0)) ) ) { + (bc.hi(1) == BCType::direction_dependent && cc_arr(i,j+1,k,1) < Real(0.0)) ) ) { vmac_arr(i,j+1,k) = cc_arr(i,j+1,k,1); } }); @@ -251,11 +257,11 @@ incflo::compute_MAC_projected_velocities ( int n = 2; const auto bc = HydroBC::getBC(i, j, k, n, domain, bc_vel_d, velbc_arr); if (k == dlo.z && ( bc.lo(2) == BCType::ext_dir || - (bc.lo(2) == BCType::direction_dependent && cc_arr(i,j,k-1,2) >= Real(0.0)) ) ) { + (bc.lo(2) == BCType::direction_dependent && cc_arr(i,j,k-1,2) > Real(0.0)) ) ) { wmac_arr(i,j,k) = cc_arr(i,j,k-1,2); } if (k == dhi.z && ( bc.hi(2) == BCType::ext_dir || - (bc.hi(2) == BCType::direction_dependent && cc_arr(i,j,k+1,2) <= Real(0.0)) ) ) { + (bc.hi(2) == BCType::direction_dependent && cc_arr(i,j,k+1,2) < Real(0.0)) ) ) { wmac_arr(i,j,k+1) = cc_arr(i,j,k+1,2); } }); diff --git a/src/incflo.cpp b/src/incflo.cpp index 7a29256df..ed47f4e5b 100644 --- a/src/incflo.cpp +++ b/src/incflo.cpp @@ -170,6 +170,13 @@ void incflo::Evolve() { if (m_verbose > 0) amrex::Print() << "Regridding...\n"; regrid(0, m_cur_time); +#ifdef INCFLO_USE_PARTICLES + // This must be done after regrid() returns: inside RemakeLevel / + // MakeNewLevelFromCoarse the ParGDB still sees the old BoxArray, + // DistributionMapping and finest_level, so a Redistribute there would + // map the particles onto the pre-regrid layout. + particleData.Redistribute(); +#endif if (m_verbose > 0 && ParallelDescriptor::IOProcessor()) { printGridSummary(amrex::OutStream(), 0, finest_level); } diff --git a/src/incflo_apply_predictor.cpp b/src/incflo_apply_predictor.cpp index 114e3a7e3..5f1d0847c 100644 --- a/src/incflo_apply_predictor.cpp +++ b/src/incflo_apply_predictor.cpp @@ -229,7 +229,10 @@ void incflo::ApplyPredictor (bool incremental_projection) // ************************************************************************************** // Update the particle positions // ************************************************************************************** - if (m_advection_type != "MOL") { + // Not during the initial pressure iterations: the fields are reset after each + // of those, but the particle positions would not be, so the particles would + // start the run m_initial_iterations*dt ahead of the fluid. + if (m_advection_type != "MOL" && !incremental_projection) { evolveTracerParticles(AMREX_D_DECL(GetVecOfConstPtrs(u_mac), GetVecOfConstPtrs(v_mac), GetVecOfConstPtrs(w_mac))); } diff --git a/src/incflo_compute_dt.cpp b/src/incflo_compute_dt.cpp index b4b5f1ca5..f96269ed2 100644 --- a/src/incflo_compute_dt.cpp +++ b/src/incflo_compute_dt.cpp @@ -52,6 +52,51 @@ void incflo::ComputeDt (int initialization, bool explicit_diffusion, double cur_ compute_vel_forces_on_level (lev, vel_forces, vel, rho, tra_o, tra); + // Explicit-diffusion bound: the largest diffusivity that is applied + // explicitly, divided by rho. We must use the strain-rate dependent + // viscosity here (not just the constant m_mu) and we must cover the tracer + // and temperature diffusivities as well, since those are updated explicitly + // with the same m_diff_type switch. Covered cells are set to zero. + MultiFab nu; + if (explicit_diffusion) { + nu.define(grids[lev], dmap[lev], 1, 0, MFInfo(), Factory(lev)); + compute_viscosity_at_level(lev, &nu, &m_leveldata[lev]->density, + &m_leveldata[lev]->velocity, geom[lev], + m_cur_time, 0); + Real mu_s_max = Real(0.0); + if (m_advect_tracer) { + for (int n = 0; n < m_ntrac; ++n) { + mu_s_max = amrex::max(mu_s_max, m_mu_s[n]); + } + } + Real mu_T_eff = m_use_temperature ? m_mu_T / m_cp : Real(0.0); +#ifdef AMREX_USE_EB + auto const& dt_flags = EBFactory(lev).getMultiEBCellFlagFab(); +#endif +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(nu,TilingIfNotGPU()); mfi.isValid(); ++mfi) { + Box const& bx = mfi.tilebox(); + Array4 const& nu_a = nu.array(mfi); + Array4 const& r = rho.const_array(mfi); +#ifdef AMREX_USE_EB + Array4 const& f = dt_flags.const_array(mfi); +#endif + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { +#ifdef AMREX_USE_EB + if (f(i,j,k).isCovered()) { nu_a(i,j,k) = Real(0.0); return; } +#endif + Real rinv = Real(1.0)/r(i,j,k); + // Velocity and temperature diffuse with eta/rho and mu_T/(rho cp); + // a non-conservative tracer diffuses with mu_s itself. + nu_a(i,j,k) = amrex::max(amrex::max(nu_a(i,j,k), mu_T_eff)*rinv, + mu_s_max*amrex::max(Real(1.0), rinv)); + }); + } + } + #ifdef AMREX_USE_EB if (!vel.isAllRegular()) { auto const& flag = EBFactory(lev).getMultiEBCellFlagFab(); @@ -72,22 +117,7 @@ void incflo::ComputeDt (int initialization, bool explicit_diffusion, double cur_ return mx; }); if (explicit_diffusion) { - diff_lev = amrex::ReduceMax(rho, flag, 0, - [=] AMREX_GPU_HOST_DEVICE (Box const& b, - Array4 const& r, - Array4 const& f) -> Real - { - Real mx = Real(-1.0); - amrex::Loop(b, [=,&mx] (int i, int j, int k) noexcept - { - if (!f(i,j,k).isCovered()) { - Real rho_inv = Real(1.0)/r(i,j,k); - mx = amrex::max(rho_inv, mx); - } - }); - return mx; - }); - diff_lev *= m_mu; + diff_lev = nu.max(0, 0, true); } // Forcing term -- old way of computing @@ -131,19 +161,7 @@ void incflo::ComputeDt (int initialization, bool explicit_diffusion, double cur_ }); if (explicit_diffusion) { - diff_lev = amrex::ReduceMax(rho, 0, - [=] AMREX_GPU_HOST_DEVICE (Box const& b, - Array4 const& r) -> Real - { - Real mx = Real(-1.0); - amrex::Loop(b, [=,&mx] (int i, int j, int k) noexcept - { - Real rho_inv = Real(1.0)/r(i,j,k); - mx = amrex::max(rho_inv, mx); - }); - return mx; - }); - diff_lev *= m_mu; + diff_lev = nu.max(0, 0, true); } // Forcing term -- old way of computing diff --git a/src/incflo_regrid.cpp b/src/incflo_regrid.cpp index 11ef44c7b..37adfc7de 100644 --- a/src/incflo_regrid.cpp +++ b/src/incflo_regrid.cpp @@ -59,10 +59,6 @@ void incflo::MakeNewLevelFromCoarse (int lev, #else macproj = std::make_unique(Geom(0,lev)); #endif - -#ifdef INCFLO_USE_PARTICLES - particleData.Redistribute(); -#endif } // Remake an existing level using provided BoxArray and DistributionMapping and @@ -120,10 +116,6 @@ void incflo::RemakeLevel (int lev, Real time, const BoxArray& ba, #else macproj = std::make_unique(Geom(0,finest_level)); #endif - -#ifdef INCFLO_USE_PARTICLES - particleData.Redistribute(); -#endif } // Rebuild macproj on demand. Note that this must not be done inside ClearLevel: diff --git a/src/prob/prob_bc.H b/src/prob/prob_bc.H index 29382e683..d85af5810 100644 --- a/src/prob/prob_bc.H +++ b/src/prob/prob_bc.H @@ -5,6 +5,95 @@ #include #include +// Normal velocity prescribed on a domain face, including the probtype-specific +// inflow profiles. Every fill functor uses this so that the direction_dependent +// inflow/outflow decision is the same for velocity, density, tracer and +// temperature. (Probtype 16 is a moving wall, not an inflow, and so is handled +// in IncfloVelFill directly.) +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +amrex::Real incflo_bc_normal_velocity (int probtype, amrex::Orientation face, + int i, int j, int k, + amrex::Box const& domain_box, amrex::Real time, + amrex::GpuArray, AMREX_SPACEDIM*2> const& bcv_vel) +{ + amrex::ignore_unused(i,j,k,time); + const int dir = face.coordDir(); + amrex::Real norm_vel = bcv_vel[face][dir]; + + if (dir == 0 && face.isLow()) + { + if (42 == probtype) + { + norm_vel = time; + } + else if (31 == probtype) + { + amrex::Real y = amrex::Real(j+0.5)*(amrex::Real(1)/domain_box.length(1)); + norm_vel = amrex::Real(6.) * y * (amrex::Real(1)-y); + } + else if (43 == probtype) + { + amrex::Real y = amrex::Real(j+0.5)*(amrex::Real(1.0)/domain_box.length(1)); + norm_vel = amrex::Real(6) * y * (amrex::Real(1)-y) - amrex::Real(1); + } +#if (AMREX_SPACEDIM == 3) + else if (311 == probtype) + { + amrex::Real z = amrex::Real(k+0.5)*(amrex::Real(1)/domain_box.length(2)); + norm_vel = amrex::Real(6.) * z * (amrex::Real(1)-z); + } + else if (41 == probtype) + { + amrex::Real z = amrex::Real(k+0.5)*(amrex::Real(1)/domain_box.length(2)); + norm_vel = amrex::Real(0.5) * z; + } +#endif + } + else if (dir == 0 && face.isHigh()) + { + if (42 == probtype) + { + norm_vel = time; + } + else if (43 == probtype) + { + amrex::Real y = amrex::Real(j+0.5)*(amrex::Real(1.0)/domain_box.length(1)); + norm_vel = amrex::Real(6) * y * (amrex::Real(1.0)-y) - amrex::Real(1); + } + } + else if (dir == 1 && face.isLow()) + { +#if (AMREX_SPACEDIM == 3) + if (32 == probtype) + { + amrex::Real z = amrex::Real(k+0.5)*(amrex::Real(1.0)/domain_box.length(2)); + norm_vel *= amrex::Real(6.) * z * (amrex::Real(1.0)-z); + } +#endif + if (322 == probtype) + { + amrex::Real x = amrex::Real(i+0.5)*(amrex::Real(1.0)/domain_box.length(0)); + norm_vel *= amrex::Real(6.) * x * (amrex::Real(1.0)-x); + } + } +#if (AMREX_SPACEDIM == 3) + else if (dir == 2 && face.isLow()) + { + if (33 == probtype) + { + amrex::Real x = amrex::Real(i+0.5)*(amrex::Real(1.0)/domain_box.length(0)); + norm_vel *= amrex::Real(6.0) * x * (amrex::Real(1.0)-x); + } + else if (333 == probtype) + { + amrex::Real y = amrex::Real(j+0.5)*(amrex::Real(1.0)/domain_box.length(1)); + norm_vel *= amrex::Real(6.0) * y * (amrex::Real(1.0)-y); + } + } +#endif + return norm_vel; +} + struct IncfloVelFill { int probtype; @@ -44,35 +133,9 @@ struct IncfloVelFill if (i < domain_box.smallEnd(0)) { int dir = 0; - amrex::Real norm_vel = bcv_vel[amrex::Orientation(amrex::Direction::x,amrex::Orientation::low)][dir]; - - // This may modify the normal velocity for specific problems - if (42 == probtype) - { - norm_vel = time; - } - else if (31 == probtype) - { - amrex::Real y = amrex::Real(j+0.5)*(amrex::Real(1)/domain_box.length(1)); - norm_vel = amrex::Real(6.) * y * (amrex::Real(1)-y); - } - else if (43 == probtype) - { - amrex::Real y = amrex::Real(j+0.5)*(amrex::Real(1.0)/domain_box.length(1)); - norm_vel = amrex::Real(6) * y * (amrex::Real(1)-y) - amrex::Real(1); - } -#if (AMREX_SPACEDIM == 3) - else if (311 == probtype) - { - amrex::Real z = amrex::Real(k+0.5)*(amrex::Real(1)/domain_box.length(2)); - norm_vel = amrex::Real(6.) * z * (amrex::Real(1)-z); - } - else if (41 == probtype) - { - amrex::Real z = amrex::Real(k+0.5)*(amrex::Real(1)/domain_box.length(2)); - norm_vel = amrex::Real(0.5) * z; - } -#endif + amrex::Real norm_vel = incflo_bc_normal_velocity(probtype, + amrex::Orientation(amrex::Direction::x,amrex::Orientation::low), + i, j, k, domain_box, time, bcv_vel); // This is a special case -- all the logic is contained here if (1101 == probtype) @@ -103,18 +166,9 @@ struct IncfloVelFill if (i > domain_box.bigEnd(0)) { int dir = 0; - amrex::Real norm_vel = bcv_vel[amrex::Orientation(amrex::Direction::x,amrex::Orientation::high)][dir]; - - // This may modify the normal velocity for specific problems - if (42 == probtype) - { - norm_vel = time; - } - else if (43 == probtype) - { - amrex::Real y = amrex::Real(j+0.5)*(amrex::Real(1.0)/domain_box.length(1)); - norm_vel = amrex::Real(6) * y * (amrex::Real(1.0)-y) - amrex::Real(1); - } + amrex::Real norm_vel = incflo_bc_normal_velocity(probtype, + amrex::Orientation(amrex::Direction::x,amrex::Orientation::high), + i, j, k, domain_box, time, bcv_vel); // This is a special case -- all the logic is contained here if (1101 == probtype) @@ -145,21 +199,9 @@ struct IncfloVelFill if (j < domain_box.smallEnd(1)) { int dir = 1; - amrex::Real norm_vel = bcv_vel[amrex::Orientation(amrex::Direction::y,amrex::Orientation::low)][dir]; - - // This may modify the normal velocity for specific problems -#if (AMREX_SPACEDIM == 3) - if (32 == probtype) - { - amrex::Real z = amrex::Real(k+0.5)*(amrex::Real(1.0)/domain_box.length(2)); - norm_vel *= amrex::Real(6.) * z * (amrex::Real(1.0)-z); - } -#endif - if (322 == probtype) - { - amrex::Real x = amrex::Real(i+0.5)*(amrex::Real(1.0)/domain_box.length(0)); - norm_vel *= amrex::Real(6.) * x * (amrex::Real(1.0)-x); - } + amrex::Real norm_vel = incflo_bc_normal_velocity(probtype, + amrex::Orientation(amrex::Direction::y,amrex::Orientation::low), + i, j, k, domain_box, time, bcv_vel); if ( (bc.lo(dir) == amrex::BCType::ext_dir) || (bc.lo(dir) == amrex::BCType::direction_dependent && norm_vel >= 0.) ) @@ -182,14 +224,9 @@ struct IncfloVelFill if (j > domain_box.bigEnd(1)) { int dir = 1; - amrex::Real norm_vel = bcv_vel[amrex::Orientation(amrex::Direction::y,amrex::Orientation::high)][dir]; - - // This may modify the normal velocity for specific problems - if (16 == probtype) - { - amrex::Real x = amrex::Real(i+0.5)*(amrex::Real(1.0)/domain_box.length(0)); - norm_vel = amrex::Real(16) * (x*x*x*x - amrex::Real(2) * x*x*x + x*x); - } + amrex::Real norm_vel = incflo_bc_normal_velocity(probtype, + amrex::Orientation(amrex::Direction::y,amrex::Orientation::high), + i, j, k, domain_box, time, bcv_vel); #if (AMREX_SPACEDIM == 3) if (1102 == probtype) @@ -207,6 +244,12 @@ struct IncfloVelFill { if ( orig_comp+nc == dir ) { vel(i,j,k,dcomp+nc) = norm_vel; + } else if (16 == probtype && orig_comp+nc == 0) { + // Burggraf lid: the prescribed profile is the *tangential* + // velocity u = 16 x^2 (1-x)^2 on the high-y wall; the normal + // velocity keeps the wall value. + amrex::Real x = amrex::Real(i+0.5)*(amrex::Real(1.0)/domain_box.length(0)); + vel(i,j,k,dcomp+nc) = amrex::Real(16) * (x*x*x*x - amrex::Real(2) * x*x*x + x*x); } else { vel(i,j,k,dcomp+nc) = bcv_vel[amrex::Orientation(amrex::Direction::y,amrex::Orientation::high)][orig_comp+nc]; } @@ -225,19 +268,9 @@ struct IncfloVelFill if (k < domain_box.smallEnd(2)) { int dir = 2; - amrex::Real norm_vel = bcv_vel[amrex::Orientation(amrex::Direction::z,amrex::Orientation::low)][dir]; - - // This may modify the normal velocity for specific problems - if (33 == probtype) - { - amrex::Real x = amrex::Real(i+0.5)*(amrex::Real(1.0)/domain_box.length(0)); - norm_vel *= amrex::Real(6.0) * x * (amrex::Real(1.0)-x); - } - else if (333 == probtype) - { - amrex::Real y = amrex::Real(j+0.5)*(amrex::Real(1.0)/domain_box.length(1)); - norm_vel *= amrex::Real(6.0) * y * (amrex::Real(1.0)-y); - } + amrex::Real norm_vel = incflo_bc_normal_velocity(probtype, + amrex::Orientation(amrex::Direction::z,amrex::Orientation::low), + i, j, k, domain_box, time, bcv_vel); // This is a special case -- all the logic is contained here if (1100 == probtype) @@ -270,7 +303,9 @@ struct IncfloVelFill if (k > domain_box.bigEnd(2)) { int dir = 2; - amrex::Real norm_vel = bcv_vel[amrex::Orientation(amrex::Direction::z,amrex::Orientation::high)][dir]; + amrex::Real norm_vel = incflo_bc_normal_velocity(probtype, + amrex::Orientation(amrex::Direction::z,amrex::Orientation::high), + i, j, k, domain_box, time, bcv_vel); // This is a special case -- all the logic is contained here if (1100 == probtype) @@ -314,7 +349,7 @@ struct IncfloDenFill AMREX_GPU_DEVICE void operator() (const amrex::IntVect& iv, amrex::Array4 const& rho, const int /*dcomp*/, const int /*numcomp*/, - amrex::GeometryData const& geom, const amrex::Real /*time*/, + amrex::GeometryData const& geom, const amrex::Real time, const amrex::BCRec* bcr, const int bcomp, const int /*orig_comp*/) const { @@ -381,48 +416,75 @@ struct IncfloDenFill else if ( (i < domain_box.smallEnd(0)) && ( (bc.lo(0) == amrex::BCType::ext_dir) || (bc.lo(0) == amrex::BCType::direction_dependent && - bcv_vel[amrex::Orientation(amrex::Direction::x,amrex::Orientation::low)][0] >= 0.) ) ) + incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::x,amrex::Orientation::low), i, j, k, domain_box, time, bcv_vel) >= 0.) ) ) { rho(i,j,k) = bcv_den[amrex::Orientation(amrex::Direction::x,amrex::Orientation::low)]; } else if ( (i > domain_box.bigEnd(0)) && ( (bc.hi(0) == amrex::BCType::ext_dir) || (bc.hi(0) == amrex::BCType::direction_dependent && - bcv_vel[amrex::Orientation(amrex::Direction::x,amrex::Orientation::high)][0] <= 0.) ) ) + incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::x,amrex::Orientation::high), i, j, k, domain_box, time, bcv_vel) <= 0.) ) ) { rho(i,j,k) = bcv_den[amrex::Orientation(amrex::Direction::x,amrex::Orientation::high)]; } + // Outflow part of a direction_dependent face: extrapolate the interior cell. + // Nobody else fills these ghost cells -- filcc has no case for + // BCType::direction_dependent -- so without this they keep stale data. + else if ( (i < domain_box.smallEnd(0)) && (bc.lo(0) == amrex::BCType::direction_dependent) ) + { + rho(i,j,k) = rho(domain_box.smallEnd(0),j,k); + } + else if ( (i > domain_box.bigEnd(0)) && (bc.hi(0) == amrex::BCType::direction_dependent) ) + { + rho(i,j,k) = rho(domain_box.bigEnd(0),j,k); + } if ( (j < domain_box.smallEnd(1)) && ( (bc.lo(1) == amrex::BCType::ext_dir) || (bc.lo(1) == amrex::BCType::direction_dependent && - bcv_vel[amrex::Orientation(amrex::Direction::y,amrex::Orientation::low)][1] >= 0.) ) ) + incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::y,amrex::Orientation::low), i, j, k, domain_box, time, bcv_vel) >= 0.) ) ) { rho(i,j,k) = bcv_den[amrex::Orientation(amrex::Direction::y,amrex::Orientation::low)]; } else if ( (j > domain_box.bigEnd(1)) && ( (bc.hi(1) == amrex::BCType::ext_dir) || (bc.hi(1) == amrex::BCType::direction_dependent && - bcv_vel[amrex::Orientation(amrex::Direction::y,amrex::Orientation::high)][1] <= 0.) ) ) + incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::y,amrex::Orientation::high), i, j, k, domain_box, time, bcv_vel) <= 0.) ) ) { rho(i,j,k) = bcv_den[amrex::Orientation(amrex::Direction::y,amrex::Orientation::high)]; } + else if ( (j < domain_box.smallEnd(1)) && (bc.lo(1) == amrex::BCType::direction_dependent) ) + { + rho(i,j,k) = rho(i,domain_box.smallEnd(1),k); + } + else if ( (j > domain_box.bigEnd(1)) && (bc.hi(1) == amrex::BCType::direction_dependent) ) + { + rho(i,j,k) = rho(i,domain_box.bigEnd(1),k); + } #if (AMREX_SPACEDIM == 3) if ( (k < domain_box.smallEnd(2)) && ( (bc.lo(2) == amrex::BCType::ext_dir) || (bc.lo(2) == amrex::BCType::direction_dependent && - bcv_vel[amrex::Orientation(amrex::Direction::z,amrex::Orientation::low)][2] >= 0.) ) ) + incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::z,amrex::Orientation::low), i, j, k, domain_box, time, bcv_vel) >= 0.) ) ) { rho(i,j,k) = bcv_den[amrex::Orientation(amrex::Direction::z,amrex::Orientation::low)]; } + else if ( (k < domain_box.smallEnd(2)) && (bc.lo(2) == amrex::BCType::direction_dependent) ) + { + rho(i,j,k) = rho(i,j,domain_box.smallEnd(2)); + } if ( (k > domain_box.bigEnd(2)) && ( (bc.hi(2) == amrex::BCType::ext_dir) || (bc.hi(2) == amrex::BCType::direction_dependent && - bcv_vel[amrex::Orientation(amrex::Direction::z,amrex::Orientation::high)][2] <= 0.) ) ) + incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::z,amrex::Orientation::high), i, j, k, domain_box, time, bcv_vel) <= 0.) ) ) { rho(i,j,k) = bcv_den[amrex::Orientation(amrex::Direction::z,amrex::Orientation::high)]; } + else if ( (k > domain_box.bigEnd(2)) && (bc.hi(2) == amrex::BCType::direction_dependent) ) + { + rho(i,j,k) = rho(i,j,domain_box.bigEnd(2)); + } #endif } }; @@ -442,7 +504,7 @@ struct IncfloTracFill AMREX_GPU_DEVICE void operator() (const amrex::IntVect& iv, amrex::Array4 const& tracer, const int /*dcomp*/, const int /*numcomp*/, - amrex::GeometryData const& geom, const amrex::Real /*time*/, + amrex::GeometryData const& geom, const amrex::Real time, const amrex::BCRec* bcr, const int bcomp, const int /*orig_comp*/) const { @@ -511,47 +573,74 @@ struct IncfloTracFill else if ( (i < domain_box.smallEnd(0)) && ( (bc.lo(0) == amrex::BCType::ext_dir) || (bc.lo(0) == amrex::BCType::direction_dependent && - bcv_vel[amrex::Orientation(amrex::Direction::x,amrex::Orientation::low)][0] >= 0.) ) ) + incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::x,amrex::Orientation::low), i, j, k, domain_box, time, bcv_vel) >= 0.) ) ) { tracer(i,j,k,n) = bcv_tra[amrex::Orientation(amrex::Direction::x,amrex::Orientation::low)][n]; } else if ( (i > domain_box.bigEnd(0)) && ( (bc.hi(0) == amrex::BCType::ext_dir) || (bc.hi(0) == amrex::BCType::direction_dependent && - bcv_vel[amrex::Orientation(amrex::Direction::x,amrex::Orientation::high)][0] <= 0.) ) ) + incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::x,amrex::Orientation::high), i, j, k, domain_box, time, bcv_vel) <= 0.) ) ) { tracer(i,j,k,n) = bcv_tra[amrex::Orientation(amrex::Direction::x,amrex::Orientation::high)][n]; } + // Outflow part of a direction_dependent face: extrapolate the interior cell. + // Nobody else fills these ghost cells -- filcc has no case for + // BCType::direction_dependent -- so without this they keep stale data. + else if ( (i < domain_box.smallEnd(0)) && (bc.lo(0) == amrex::BCType::direction_dependent) ) + { + tracer(i,j,k,n) = tracer(domain_box.smallEnd(0),j,k,n); + } + else if ( (i > domain_box.bigEnd(0)) && (bc.hi(0) == amrex::BCType::direction_dependent) ) + { + tracer(i,j,k,n) = tracer(domain_box.bigEnd(0),j,k,n); + } if ( (j < domain_box.smallEnd(1)) && ( (bc.lo(1) == amrex::BCType::ext_dir) || (bc.lo(1) == amrex::BCType::direction_dependent && - bcv_vel[amrex::Orientation(amrex::Direction::y,amrex::Orientation::low)][1] >= 0.) ) ) + incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::y,amrex::Orientation::low), i, j, k, domain_box, time, bcv_vel) >= 0.) ) ) { tracer(i,j,k,n) = bcv_tra[amrex::Orientation(amrex::Direction::y,amrex::Orientation::low)][n]; } else if ( (j > domain_box.bigEnd(1)) && ( (bc.hi(1) == amrex::BCType::ext_dir) || (bc.hi(1) == amrex::BCType::direction_dependent && - bcv_vel[amrex::Orientation(amrex::Direction::y,amrex::Orientation::high)][1] <= 0.) ) ) + incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::y,amrex::Orientation::high), i, j, k, domain_box, time, bcv_vel) <= 0.) ) ) { tracer(i,j,k,n) = bcv_tra[amrex::Orientation(amrex::Direction::y,amrex::Orientation::high)][n]; } + else if ( (j < domain_box.smallEnd(1)) && (bc.lo(1) == amrex::BCType::direction_dependent) ) + { + tracer(i,j,k,n) = tracer(i,domain_box.smallEnd(1),k,n); + } + else if ( (j > domain_box.bigEnd(1)) && (bc.hi(1) == amrex::BCType::direction_dependent) ) + { + tracer(i,j,k,n) = tracer(i,domain_box.bigEnd(1),k,n); + } #if (AMREX_SPACEDIM == 3) if ( (k < domain_box.smallEnd(2)) && ( (bc.lo(2) == amrex::BCType::ext_dir) || (bc.lo(2) == amrex::BCType::direction_dependent && - bcv_vel[amrex::Orientation(amrex::Direction::z,amrex::Orientation::low)][2] >= 0.) ) ) + incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::z,amrex::Orientation::low), i, j, k, domain_box, time, bcv_vel) >= 0.) ) ) { tracer(i,j,k,n) = bcv_tra[amrex::Orientation(amrex::Direction::z,amrex::Orientation::low)][n]; } else if ( (k > domain_box.bigEnd(2)) && ( (bc.hi(2) == amrex::BCType::ext_dir) || (bc.hi(2) == amrex::BCType::direction_dependent && - bcv_vel[amrex::Orientation(amrex::Direction::z,amrex::Orientation::high)][2] <= 0.) ) ) + incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::z,amrex::Orientation::high), i, j, k, domain_box, time, bcv_vel) <= 0.) ) ) { tracer(i,j,k,n) = bcv_tra[amrex::Orientation(amrex::Direction::z,amrex::Orientation::high)][n]; } + else if ( (k < domain_box.smallEnd(2)) && (bc.lo(2) == amrex::BCType::direction_dependent) ) + { + tracer(i,j,k,n) = tracer(i,j,domain_box.smallEnd(2),n); + } + else if ( (k > domain_box.bigEnd(2)) && (bc.hi(2) == amrex::BCType::direction_dependent) ) + { + tracer(i,j,k,n) = tracer(i,j,domain_box.bigEnd(2),n); + } #endif } } @@ -572,7 +661,7 @@ struct IncfloTempFill AMREX_GPU_DEVICE void operator() (const amrex::IntVect& iv, amrex::Array4 const& temperature, const int /*dcomp*/, const int /*numcomp*/, - amrex::GeometryData const& geom, const amrex::Real /*time*/, + amrex::GeometryData const& geom, const amrex::Real time, const amrex::BCRec* bcr, const int bcomp, const int /*orig_comp*/) const { @@ -639,47 +728,74 @@ struct IncfloTempFill else if ( (i < domain_box.smallEnd(0)) && ( (bc.lo(0) == amrex::BCType::ext_dir) || (bc.lo(0) == amrex::BCType::direction_dependent && - bcv_vel[amrex::Orientation(amrex::Direction::x,amrex::Orientation::low)][0] >= 0.) ) ) + incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::x,amrex::Orientation::low), i, j, k, domain_box, time, bcv_vel) >= 0.) ) ) { temperature(i,j,k) = bcv_tem[amrex::Orientation(amrex::Direction::x,amrex::Orientation::low)]; } else if ( (i > domain_box.bigEnd(0)) && ( (bc.hi(0) == amrex::BCType::ext_dir) || (bc.hi(0) == amrex::BCType::direction_dependent && - bcv_vel[amrex::Orientation(amrex::Direction::x,amrex::Orientation::high)][0] <= 0.) ) ) + incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::x,amrex::Orientation::high), i, j, k, domain_box, time, bcv_vel) <= 0.) ) ) { temperature(i,j,k) = bcv_tem[amrex::Orientation(amrex::Direction::x,amrex::Orientation::high)]; } + // Outflow part of a direction_dependent face: extrapolate the interior cell. + // Nobody else fills these ghost cells -- filcc has no case for + // BCType::direction_dependent -- so without this they keep stale data. + else if ( (i < domain_box.smallEnd(0)) && (bc.lo(0) == amrex::BCType::direction_dependent) ) + { + temperature(i,j,k) = temperature(domain_box.smallEnd(0),j,k); + } + else if ( (i > domain_box.bigEnd(0)) && (bc.hi(0) == amrex::BCType::direction_dependent) ) + { + temperature(i,j,k) = temperature(domain_box.bigEnd(0),j,k); + } if ( (j < domain_box.smallEnd(1)) && ( (bc.lo(1) == amrex::BCType::ext_dir) || (bc.lo(1) == amrex::BCType::direction_dependent && - bcv_vel[amrex::Orientation(amrex::Direction::y,amrex::Orientation::low)][1] >= 0.) ) ) + incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::y,amrex::Orientation::low), i, j, k, domain_box, time, bcv_vel) >= 0.) ) ) { temperature(i,j,k) = bcv_tem[amrex::Orientation(amrex::Direction::y,amrex::Orientation::low)]; } else if ( (j > domain_box.bigEnd(1)) && ( (bc.hi(1) == amrex::BCType::ext_dir) || (bc.hi(1) == amrex::BCType::direction_dependent && - bcv_vel[amrex::Orientation(amrex::Direction::y,amrex::Orientation::high)][1] <= 0.) ) ) + incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::y,amrex::Orientation::high), i, j, k, domain_box, time, bcv_vel) <= 0.) ) ) { temperature(i,j,k) = bcv_tem[amrex::Orientation(amrex::Direction::y,amrex::Orientation::high)]; } + else if ( (j < domain_box.smallEnd(1)) && (bc.lo(1) == amrex::BCType::direction_dependent) ) + { + temperature(i,j,k) = temperature(i,domain_box.smallEnd(1),k); + } + else if ( (j > domain_box.bigEnd(1)) && (bc.hi(1) == amrex::BCType::direction_dependent) ) + { + temperature(i,j,k) = temperature(i,domain_box.bigEnd(1),k); + } #if (AMREX_SPACEDIM == 3) if ( (k < domain_box.smallEnd(2)) && ( (bc.lo(2) == amrex::BCType::ext_dir) || (bc.lo(2) == amrex::BCType::direction_dependent && - bcv_vel[amrex::Orientation(amrex::Direction::z,amrex::Orientation::low)][2] >= 0.) ) ) + incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::z,amrex::Orientation::low), i, j, k, domain_box, time, bcv_vel) >= 0.) ) ) { temperature(i,j,k) = bcv_tem[amrex::Orientation(amrex::Direction::z,amrex::Orientation::low)]; } else if ( (k > domain_box.bigEnd(2)) && ( (bc.hi(2) == amrex::BCType::ext_dir) || (bc.hi(2) == amrex::BCType::direction_dependent && - bcv_vel[amrex::Orientation(amrex::Direction::z,amrex::Orientation::high)][2] <= 0.) ) ) + incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::z,amrex::Orientation::high), i, j, k, domain_box, time, bcv_vel) <= 0.) ) ) { temperature(i,j,k) = bcv_tem[amrex::Orientation(amrex::Direction::z,amrex::Orientation::high)]; } + else if ( (k < domain_box.smallEnd(2)) && (bc.lo(2) == amrex::BCType::direction_dependent) ) + { + temperature(i,j,k) = temperature(i,j,domain_box.smallEnd(2)); + } + else if ( (k > domain_box.bigEnd(2)) && (bc.hi(2) == amrex::BCType::direction_dependent) ) + { + temperature(i,j,k) = temperature(i,j,domain_box.bigEnd(2)); + } #endif } }; diff --git a/src/prob/prob_init_fluid.cpp b/src/prob/prob_init_fluid.cpp index 5762b87cd..b2aef967b 100644 --- a/src/prob/prob_init_fluid.cpp +++ b/src/prob/prob_init_fluid.cpp @@ -1035,7 +1035,7 @@ void incflo::init_plane_poiseuille (Box const& vbx, Box const& /*gbx*/, if (nt > 2 && i <= dhi.x*3/4) tracer(i,j,k,2) = 3.0; }); } - else if (42 == m_probtype || 43 == m_probtype) + else if (42 == m_probtype) { ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { @@ -1046,6 +1046,25 @@ void incflo::init_plane_poiseuille (Box const& vbx, Box const& /*gbx*/, } }); } + else if (43 == m_probtype) + { + // Start from the zero-net-flux in/out profile that IncfloVelFill imposes on + // the x faces (prob_bc.H), so that the outflow part of the direction_dependent + // boundary has a non-zero velocity when the initial projection enforces + // in/out solvability. + ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real y = Real(j+0.5)*dyinv; + AMREX_D_TERM(vel(i,j,k,0) = Real(6.0) * y * (Real(1.0)-y) - Real(1.0);, + vel(i,j,k,1) = Real(0.0);, + vel(i,j,k,2) = Real(0.0);); + + const int nt = tracer.nComp(); + for (int n = 0; n < nt; ++n) { + tracer(i,j,k,n) = 0.0; + } + }); + } else if (32 == m_probtype) { Real v = m_ic_v;