use masking of velocity diffusion based on terrain cells

amr-wind/core/FieldRepo.H

@@ -310,6 +310,12 @@ public:
         int nghost = 0,
         FieldLoc floc = FieldLoc::CELL) const;
 
+    std::unique_ptr<IntScratchField> create_int_scratch_field(
+        const std::string& name,
+        int ncomp = 1,
+        int nghost = 0,
+        FieldLoc floc = FieldLoc::CELL) const;
+
     /** Create a scratch field
      *
      *  ScratchField is a temporary field used to compute and store intermediate
@@ -332,6 +338,9 @@ public:
     std::unique_ptr<ScratchField> create_scratch_field_on_host(
         int ncomp = 1, int nghost = 0, FieldLoc floc = FieldLoc::CELL) const;
 
+    std::unique_ptr<IntScratchField> create_int_scratch_field(
+        int ncomp = 1, int nghost = 0, FieldLoc floc = FieldLoc::CELL) const;
+
     std::unique_ptr<IntScratchField> create_int_scratch_field_on_host(
         const std::string& name,
         int ncomp = 1,

amr-wind/core/FieldRepo.cpp

@@ -368,6 +368,36 @@ std::unique_ptr<ScratchField> FieldRepo::create_scratch_field_on_host(
     return create_scratch_field_on_host(
         "scratch_field_host", ncomp, nghost, floc);
 }
+std::unique_ptr<IntScratchField> FieldRepo::create_int_scratch_field(
+    const std::string& name,
+    const int ncomp,
+    const int nghost,
+    const FieldLoc floc) const
+{
+    BL_PROFILE("amr-wind::FieldRepo::create_int_scratch_field");
+    if (!m_is_initialized) {
+        amrex::Abort(
+            "Integer scratch field creation is not permitted before mesh is "
+            "initialized");
+    }
+    std::unique_ptr<IntScratchField> field(
+        new IntScratchField(*this, name, ncomp, nghost, floc));
+
+    for (int lev = 0; lev <= m_mesh.finestLevel(); ++lev) {
+        const auto ba =
+            amrex::convert(m_mesh.boxArray(lev), field_impl::index_type(floc));
+
+        field->m_data.emplace_back(
+            ba, m_mesh.DistributionMap(lev), ncomp, nghost, amrex::MFInfo(),
+            *(m_leveldata[lev]->m_int_fact));
+    }
+    return field;
+}
+std::unique_ptr<IntScratchField> FieldRepo::create_int_scratch_field(
+    const int ncomp, const int nghost, const FieldLoc floc) const
+{
+    return create_int_scratch_field("scratch_field", ncomp, nghost, floc);
+}
 std::unique_ptr<IntScratchField> FieldRepo::create_int_scratch_field_on_host(
     const std::string& name,
     const int ncomp,

amr-wind/equation_systems/DiffusionOps.cpp

@@ -25,7 +25,7 @@ DiffSolverIface<LinOp>::DiffSolverIface(
     iapply.setMaxCoarseningLevel(0);
 
     const auto& mesh = m_pdefields.repo.mesh();
-    if (!has_overset) {
+    if (!has_overset && !fields.repo.int_field_exists("terrain_blank")) {
         m_solver.reset(new LinOp(
             mesh.Geom(0, mesh.finestLevel()),
             mesh.boxArray(0, mesh.finestLevel()),
@@ -35,7 +35,23 @@ DiffSolverIface<LinOp>::DiffSolverIface(
             mesh.boxArray(0, mesh.finestLevel()),
             mesh.DistributionMap(0, mesh.finestLevel()), iapply));
     } else {
-        auto imask = fields.repo.get_int_field("mask_cell").vec_const_ptrs();
+        auto combo_mask =
+            fields.repo.create_int_scratch_field("combo_mask", 1, 1);
+        combo_mask->setVal(1);
+        if (fields.repo.int_field_exists("terrain_blank")) {
+            auto& terrain_blank = fields.repo.get_int_field("terrain_blank");
+            for (int i = 0; i < mesh.finestLevel() + 1; ++i) {
+                amrex::iMultiFab::Saxpy(
+                    (*combo_mask)(i), -1, terrain_blank(i), 0, 0, 1,
+                    amrex::IntVect(1));
+            }
+        }
+        const auto& overset_mask = fields.repo.get_int_field("mask_cell");
+        for (int i = 0; i < mesh.finestLevel() + 1; ++i) {
+            amrex::iMultiFab::Multiply(
+                (*combo_mask)(i), overset_mask(i), 0, 0, 1, 1);
+        }
+        auto imask = (*combo_mask).vec_const_ptrs();
         m_solver.reset(new LinOp(
             mesh.Geom(0, mesh.finestLevel()),
             mesh.boxArray(0, mesh.finestLevel()),

amr-wind/equation_systems/icns/icns_diffusion.H

@@ -9,6 +9,7 @@
 #include "amr-wind/equation_systems/DiffusionOps.H"
 #include "amr-wind/equation_systems/icns/icns.H"
 #include "amr-wind/utilities/console_io.H"
+#include "amr-wind/core/field_ops.H"
 #include "AMReX_REAL.H"
 
 using namespace amrex::literals;
@@ -72,7 +73,7 @@
         const auto& bchi = diffusion::get_diffuse_tensor_bc(
             m_pdefields.field, amrex::Orientation::high);
 
-        if (!has_overset) {
+        if (!has_overset && !fields.repo.int_field_exists("terrain_blank")) {
 
             m_solver_scalar = std::make_unique<amrex::MLABecLaplacian>(
                 mesh.Geom(0, mesh.finestLevel()),
@@ -89,8 +90,24 @@
                     amrex::FabFactory<amrex::MLABecLaplacian::FAB> const*>(),
                 AMREX_SPACEDIM);
         } else {
-            auto imask =
-                fields.repo.get_int_field("mask_cell").vec_const_ptrs();
+            auto combo_mask =
+                fields.repo.create_int_scratch_field("combo_mask", 1, 1);
+            combo_mask->setVal(1);
+            if (fields.repo.int_field_exists("terrain_blank")) {
+                auto& terrain_blank =
+                    fields.repo.get_int_field("terrain_blank");
+                for (int i = 0; i < mesh.finestLevel() + 1; ++i) {
+                    amrex::iMultiFab::Saxpy(
+                        (*combo_mask)(i), -1, terrain_blank(i), 0, 0, 1,
+                        amrex::IntVect(1));
+                }
+            }
+            const auto& overset_mask = fields.repo.get_int_field("mask_cell");
+            for (int i = 0; i < mesh.finestLevel() + 1; ++i) {
+                amrex::iMultiFab::Multiply(
+                    (*combo_mask)(i), overset_mask(i), 0, 0, 1, 1);
+            }
+            auto imask = (*combo_mask).vec_const_ptrs();
             m_solver_scalar = std::make_unique<amrex::MLABecLaplacian>(
                 mesh.Geom(0, mesh.finestLevel()),
                 mesh.boxArray(0, mesh.finestLevel()),
@@ -290,7 +307,8 @@
             m_pdefields.field, amrex::Orientation::high);
 
         for (int i = 0; i < AMREX_SPACEDIM; ++i) {
-            if (!has_overset) {
+            if (!has_overset &&
+                !fields.repo.int_field_exists("terrain_blank")) {
                 m_solver_scalar[i] = std::make_unique<amrex::MLABecLaplacian>(
                     mesh.Geom(0, mesh.finestLevel()),
                     mesh.boxArray(0, mesh.finestLevel()),
@@ -300,8 +318,25 @@
                     mesh.boxArray(0, mesh.finestLevel()),
                     mesh.DistributionMap(0, mesh.finestLevel()), iapply);
             } else {
-                auto imask =
-                    fields.repo.get_int_field("mask_cell").vec_const_ptrs();
+                auto combo_mask =
+                    fields.repo.create_int_scratch_field("combo_mask", 1, 1);
+                combo_mask->setVal(1);
+                if (fields.repo.int_field_exists("terrain_blank")) {
+                    auto& terrain_blank =
+                        fields.repo.get_int_field("terrain_blank");
+                    for (int i = 0; i < mesh.finestLevel() + 1; ++i) {
+                        amrex::iMultiFab::Saxpy(
+                            (*combo_mask)(i), -1, terrain_blank(i), 0, 0, 1,
+                            amrex::IntVect(1));
+                    }
+                }
+                const auto& overset_mask =
+                    fields.repo.get_int_field("mask_cell");
+                for (int i = 0; i < mesh.finestLevel() + 1; ++i) {
+                    amrex::iMultiFab::Multiply(
+                        (*combo_mask)(i), overset_mask(i), 0, 0, 1, 1);
+                }
+                auto imask = (*combo_mask).vec_const_ptrs();
                 m_solver_scalar[i] = std::make_unique<amrex::MLABecLaplacian>(
                     mesh.Geom(0, mesh.finestLevel()),
                     mesh.boxArray(0, mesh.finestLevel()),