Would like some idea on how to integrate EBGeometry on a EB_CNS -based AMReX code.

[milin333martin]

I have been trying to add the EBGeometry capability to our AMReX code, that is based of EB_CNS. The simulation we essentially want to run is a flow from one corner of the domain (Base) propagates downstream and enters a narrow passage with some holes made of a few shapes. This narrow passage is made using an STL file. I was facing several issues doing it using the native STL capability of AMReX and wanted to try adapting the EBGeometry.

The issue that I come across now is that I tried adding the lines of code based on AMRex_PaintEB into a function called init_eb_geometry (Based on EB_CNS). The code compiles correctly but it doesnt proceed after the input file is read, the code just remains stagnant. Here, are the lines of the code so that you maybe able to give me a nudge in the correct direction.

/*!
  @brief This is an AMReX-capable version of the EBGeometry BVH accelerator. It
  is templated as T, BV, K which indicate the EBGeometry precision, bounding volume, 
  tree degree, and meta-data type for the triangles (defaults to an integer)
*/
```
```
template <class T, class BV, size_t K, class Meta = double>
class AMReXSDF
{
public:
  /*!
    @brief Shortcut for DCEL face type
  */
  using Face = EBGeometry::DCEL::FaceT<T, Meta>;

  /*!
    @brief Full constructor.
    @param[in] a_filename File name. Must be an STL file.
  */
  AMReXSDF(const std::string a_filename)
  {
    // Read in the mesh into a DCEL mesh and partition it into a bounding volume hierarchy
    auto mesh = EBGeometry::Parser::readIntoDCEL<T, Meta>(a_filename);
    m_sdf     = std::make_shared<EBGeometry::FastCompactMeshSDF<T, Meta, BV, K>>(mesh);

    // Set the meta-data for all facets to their "index", i.e. position in the list of facets
    auto& faces = mesh->getFaces();
    for (size_t i = 0; i < faces.size(); i++) {
      faces[i]->getMetaData() = 1.0 * i;
    }
  }

  /*!
    @brief Copy constructor.
    @param[in] a_other Other SDF.
  */
  AMReXSDF(const AMReXSDF& a_other)
  {
    this->m_sdf = a_other.m_sdf;
  }

  /*!
    @brief AMReX's implicit function definition.
  */
  Real
  operator()(AMREX_D_DECL(Real x, Real y, Real z)) const noexcept
  {
    return m_sdf->value(EBGeometry::Vec3T<T>(x, y, z));
  };

  /*!
    @brief Also an AMReX implicit function implementation
  */
  inline Real
  operator()(const RealArray& p) const noexcept
  {
    return this->operator()(AMREX_D_DECL(p[0], p[1], p[2]));
  }

  /*!
    @brief Get the face(s) that are closest to the input point. 
  */
  inline std::vector<std::pair<std::shared_ptr<const Face>, T>>
  getClosestFaces(AMREX_D_DECL(Real x, Real y, Real z)) const noexcept
  {
    return m_sdf->getClosestFaces(EBGeometry::Vec3T<T>(x, y, z), true);
  }

protected:
  /*!
    @brief DCEL mesh represented as a BVH of its facets, exposed as an implicit function. 
  */
  std::shared_ptr<EBGeometry::FastCompactMeshSDF<T, Meta, BV, K>> m_sdf;
};


AMREX_GPU_HOST
AMREX_FORCE_INLINE
void
init_eb_geometry (const amrex::Geometry& geom, const int max_coarsening_level)
{

	std::string filename = "./Shapes1.stl";
	bool use_bvh         = true;
	int num_coarsen_opt  = 0;
	
	int n_cell = 64;
	
	//int nCells = prob.n_cells

	ParmParse pp("eb2");
	
	pp.query("use_bvh", use_bvh);
	pp.query("num_coarsen_opt", num_coarsen_opt);

	// Define RealBox: Adjust this as needed for your case
	RealBox rb({0.11051, 0.0085021, 0.0085021}, 
	       {0.17051, 0.0115021, 0.0115021});

	Array<int, AMREX_SPACEDIM> is_periodic{false, false, false};
	Geometry::Setup(&rb, 0, is_periodic.data());

	// Define the problem domain (Ensure n_cell is properly set elsewhere)
	//Box domain(IntVect(0), IntVect(n_cell - 1));
	//geom.define(domain);

  	// Create our signed distance function. K is the tree degree while T is the
	// EBGeometry precision.
	constexpr int K = 4;

	using T    = float;
	using Vec3 = EBGeometry::Vec3T<T>;
	using BV   = EBGeometry::BoundingVolumes::AABBT<T>;

	AMReXSDF<T, BV, K> sdf(filename);

	// Create the EB geometry
	auto gshop = EB2::makeShop(sdf);
	EB2::Build(gshop, geom, 0, 0, 1, true, true, num_coarsen_opt);
	

}`
```

Thank you so much for your time.

[rmrsk]

Hi, I'll take a look at it, but can you also share the STL file?

[milin333martin]

Thank you so much for the quick response. I am attaching my STL file here.
Shapes1.zip

Once again, thank you

[milin333martin]

One thing that I would like to mention is that I have previously tried using this STL file using the native STL features of AMReX. One issue that I keep coming across is the error of "amrex::Abort::0::Failed to build required coarse EB level 1 !!!"
I sometimes get away with this by making my base mesh finer, but for my current case it becomes severely memory restrictive.

[rmrsk]

Hi, and sorry, but work is nuts right now and I just havent had time to check this out properly....

That being said, have you checked that the STL file is a proper geometry. I.e., that it is watertight with consistently oriented faces and no self-intersections?

[milin333martin]

Thank you so much for making time amidst your busy schedule.
It was a great comment on the 'properness' of the STL geometry. I realized that there where several self-interactions. I am currently trying ways to rectify that and then try the code again. I will update you on how things went.

[milin333martin]

Hi again,

B3.zip

This STL is now a watertight and proper geometry with no self-intersections. I tried using this with the previously posted code but I get the same problem, where the code is not proceeding to the time loop. Secondly, I also noticed that if I use this STL file directly using the native STL tools of AMReX, then too the code remains idle after showing the number of triangles and the Min and Max positions.

AMREX_GPU_HOST
AMREX_FORCE_INLINE
void init_eb_geometry(const amrex::Geometry& geom_in, const int max_coarsening_level)
{
    // Create a modifiable copy of the input geometry.
    amrex::Geometry geom = geom_in;

    std::string filename = "./B3.stl";
    bool use_bvh = true;
    int num_coarsen_opt = 0;
    //int n_cell = 128;

    ParmParse pp("eb2");
    //pp.get("n_cell", n_cell);
    pp.query("use_bvh", use_bvh);
    pp.query("num_coarsen_opt", num_coarsen_opt);

    // Setup the RealBox (adjust these values as needed).
    RealBox rb({0.11051, 0.0085021, 0.0085021},
               {0.17051, 0.0115021, 0.0115021});
    Print() << "Working so far 1...\n";

    Array<int, AMREX_SPACEDIM> is_periodic{false, false, false};
    Print() << "Working so far 2...\n";

    // Configure global geometry settings.
    Geometry::Setup(&rb, 0, is_periodic.data());
    Print() << "Working so far 3...\n";

    // Retrieve the domain from the copied geometry instance.
    const Box& domain = geom.Domain();
    // Now call define() on the non-const copy.
    geom.define(domain);

    // Set up the EB geometry from the STL file.
    constexpr int K = 4;
    using T = float;
    using BV = EBGeometry::BoundingVolumes::AABBT<T>;
    Print() << "Working so far 4...\n";

    AMReXSDF<T, BV, K> sdf(filename);
    Print() << "Working so far 6...\n";

    auto gshop = EB2::makeShop(sdf);
    Print() << "Working so far 7...\n";

    EB2::Build(gshop, geom, 0, 0, 1, true, true, num_coarsen_opt);
    Print() << "Working so far 8...\n";
}

mpirun -np 48 ./CNS3d.gnu.MPI.ex inputs 
Initializing AMReX (24.02-30-g2ecafcff4013-dirty)...
MPI initialized with 48 MPI processes
MPI initialized with thread support level 0
AMReX (24.02-30-g2ecafcff4013-dirty) initialized
Successfully read inputs file ... 
Working so far 1...
Working so far 2...
Working so far 3...
Working so far 4...
Working so far 6...
Working so far 7...

.....Idles

[rmrsk]

Hi again,

First: I'm very sorry for responding so ridiculously late - I mean well and then drop the ball...

I did take a look at this and could generate the geometry rather easily. Maybe you can verify the same steps:

1. I added the following lines to EBGeometry/Exec/Examples/AMReX_DCEL, just below the snippet that would read which_geom=7

      else if (which_geom == 8) { // B3 geometry
        rb       = RealBox({-0.01, -0.001, -0.001}, {0.025, 0.003, 0.003});
        filename = "../Resources/B3.stl";
      }

2. Add the B3.stl file to the Resources folder containing the other STLs.

After that I was able to generate the geometry using the instructions in the README.md file. Specifically, I ran

mpirun -np 8 main3d.gnu.TEST.MPI.ex eb2.cover_multiple_cuts=1 which_geom=8

Here is an image of the geometry

Maybe there is something wrong other than the STL itself (or the method of generation)....? Can you check if the above steps resolved anything before we proceed further?

[milin333martin]

Thank you so much for your suggestions @rmrsk , I am grateful for the time that you spare for me :)

I was successfully able to build the STL using the AMReX_DCEL. Actually even before I was able to build it. Which is why I wanted to use that strategy to build the geometry in the EB_CNS based code that we have.

One thing that i forgot to mention previously is that we tried the conventional way by which EB_CNS cases read the STL, as shown here:

eb2.geom_type = stl
eb2.stl_file = B3.stl # CylinderTube.stl #
eb2.stl_scale = 1 # Adjust as needed
eb2.stl_center = 0.045 0.0001 0.0001 # Adjust to your desired center
eb2.stl_reverse_normal = 0 # Set to 1 if you need to reverse the normals
eb2.cover_multiple_cuts = 1
eb2.num_coarsen_opt = 1

But when I did that way I would either end up getting things like:
AMReX EB: Iter. 32 fixed 26 small cells
4::Assertion `iter < maxiter' failed, file "/scratch/user/milin333martin/Milin/amrex_new//Src/EB/AMReX_EB2_Level.H", line 501, Msg: "EB: failed to fix small cells" !!!

or things like:
amrex::Abort::0::Failed to build required coarse EB level 1 !!!

Based on other github discussions, I tried fiddling around and when the build seems to work, it usually leads to a deadlock during build. Which is why I wanted to try your method.

The way we do it is by calling the cns_init_eb2.cpp as shown here:
CNS_init_eb2.zip

Although I have been able to successfully compile this part after including the EBGeometry's functionalities, I still end up with EB building deadlocks. It would be great if you could direct me to an example in the original EB_CNS test case wherein the EBGeometry is used to build STLs, if possible (Again, I am so sorry for using up so much of your time, I really appreciate the guidance :))

[rmrsk]

The error messages are due to AMReX itself as it fails to coarsen your geometry. There might be ways around this, and you should probably reach out to the AMReX team regarding the deadlocks and error messages. So far, I can't see anything that has to do with an "improper geometry".

The deeper issue is that AMReX doesn't permit multi-valued cells, which appear during coarsening. Unfortunately I'm not very hopeful that resolving the deadlock will get you anywhere, because failure to coarsen geometries will 1) break geometric multigrid, and 2) prevent proper subcycling, which might affect your code/results.