@@ -20,6 +20,8 @@ public:
2020
2121 using cMF = FabArray<BaseFab<GpuComplex<typename MF ::value_type>>>;
2222
23+ Stokes () = default ;
24+
2325 /* *
2426 * \brief Construct a Stokes solver with explicit boundary types.
2527 *
@@ -31,41 +33,61 @@ public:
3133 Stokes (Geometry const & geom,
3234 Array<std::pair<Boundary,Boundary>,AMREX_SPACEDIM > const & bc)
3335 requires (IsFabArray_v<MF >)
34- : m_domain_lo(geom.Domain().smallEnd()),
35- m_geom (detail::shift_geom(geom)),
36- m_bc(bc)
3736 {
37+ define (geom, bc);
38+ }
39+
40+ /* *
41+ * \brief Construct a purely periodic Stokes solver.
42+ *
43+ * \param geom Periodic geometry (all directions must be periodic).
44+ */
45+ explicit Stokes (Geometry const & geom)
46+ requires (IsFabArray_v<MF >)
47+ {
48+ define (geom);
49+ }
50+
51+ /* *
52+ * \brief Define a Stokes solver with explicit boundary types.
53+ *
54+ * Only periodic BCs are supported.
55+ *
56+ * \param geom Geometry describing the domain and metric information.
57+ * \param bc Pair of boundary descriptors (low/high) per coordinate direction.
58+ */
59+ void define (Geometry const & geom,
60+ Array<std::pair<Boundary,Boundary>,AMREX_SPACEDIM > const & bc)
61+ requires (IsFabArray_v<MF >)
62+ {
63+ m_domain_lo = geom.Domain ().smallEnd ();
64+ m_geom = detail::shift_geom (geom);
65+ m_bc = bc;
66+
3867 bool all_periodic = true ;
3968 for (int idim = 0 ; idim < AMREX_SPACEDIM ; ++idim) {
4069 all_periodic = all_periodic
4170 && (bc[idim].first == Boundary::periodic)
4271 && (bc[idim].second == Boundary::periodic);
4372 }
44- if (all_periodic) {
45- m_r2c = std::make_unique<R2C <typename MF ::value_type>>(m_geom.Domain ());
46- } else {
73+ if (!all_periodic || !m_geom.isAllPeriodic ()) {
4774 amrex::Abort (" FFT::Stokes: only supports periodic BC" );
4875 }
76+
77+ m_r2c = std::make_unique<R2C <typename MF ::value_type>>(m_geom.Domain ());
4978 }
5079
5180 /* *
52- * \brief Construct a purely periodic Stokes solver.
81+ * \brief Define a purely periodic Stokes solver.
5382 *
5483 * \param geom Periodic geometry (all directions must be periodic).
5584 */
56- explicit Stokes (Geometry const & geom)
85+ void define (Geometry const & geom)
5786 requires (IsFabArray_v<MF >)
58- : m_domain_lo(geom.Domain().smallEnd()),
59- m_geom(detail::shift_geom(geom)),
60- m_bc{AMREX_D_DECL (std::make_pair (Boundary::periodic,Boundary::periodic),
61- std::make_pair (Boundary::periodic,Boundary::periodic),
62- std::make_pair (Boundary::periodic,Boundary::periodic))}
6387 {
64- if (m_geom.isAllPeriodic ()) {
65- m_r2c = std::make_unique<R2C <typename MF ::value_type>>(m_geom.Domain ());
66- } else {
67- amrex::Abort (" FFT::Stokes: only supports periodic BC" );
68- }
88+ define (geom, {AMREX_D_DECL (std::make_pair (Boundary::periodic,Boundary::periodic),
89+ std::make_pair (Boundary::periodic,Boundary::periodic),
90+ std::make_pair (Boundary::periodic,Boundary::periodic))});
6991 }
7092
7193 /* *
@@ -120,6 +142,10 @@ void Stokes<MF>::solve (MF& U,
120142
121143 using T = typename MF ::value_type;
122144
145+ if (!m_r2c) {
146+ amrex::Abort (" FFT::Stokes::solve called before define()" );
147+ }
148+
123149 AMREX_ASSERT (p.ixType () == IndexType::TheCellType () &&
124150 AMREX_D_TERM (U.ixType () == IndexType (IntVect::TheDimensionVector (0 )) &&
125151 rhsx.ixType () == U.ixType (),
0 commit comments