Embedded boundaries with Fortran and C interfaces

[drjosephbrock]

I am trying to extend the Advection_octree_F example to use embedded boundaries. My main solver is written in fortran and I want to implement embedded boundaries, which as of right now is not available through the fortran interface APIs. My initial idea is to set up the domain and options through the fortran code (similar to Advection_octree_F) and then include an interface to initialize the EB to tag cells. Something like;

interface
subroutine EB_setup(geometry, ba, dm, max_coarsening_level) bind(C, name="EB_setup")
    import c_ptr
    type(c_ptr), intent(IN) :: geometry
    type(c_ptr), intent(IN) :: ba
    type(c_ptr), intent(IN) :: dm
    integer, intent(IN) :: max_coarsening_level
end subroutine EB_setup
end interface

and,

`
#ifdef __cplusplus
extern "C" {
#endif

void EB_setup(const amrex::Geometry& geom, const amrex::BoxArray& ba, const amrex::DistributionMapping& dm, const int& max_level);

#ifdef __cplusplus
}
#endif
`

the call from the main routine would look like;

    call EB_setup(amrex_geom(0)%p, ba%p, dm%p, max_level)

Is this approach possible?

[WeiqunZhang]

In principle, yes.

[drjosephbrock]

Do I need to loop through amrex_geom array and initialize for each or should I pass the 0th or max level?

[WeiqunZhang]

Usually one calls EB2:Build with the Geometry at the max level at the beginning. This is done only once. This step does not use any BoxArrays or DistributionMappings. (See https://amrex-codes.github.io/amrex/docs_html/EB.html#initializing-the-geometric-database)

Then one builds EB data for each AMR level. During regridding, this needs to be rebuilt. (See https://amrex-codes.github.io/amrex/docs_html/EB.html#ebfarrayboxfactory)

[drjosephbrock]

hmmm. Seems like I am getting seg faults when trying to build EB data for each AMR level.

! Make a new level from scratch and put the data in phi_new.
! Note tha phi_old contains no valid data after this.
subroutine my_make_new_level_from_scratch (lev, time, pba, pdm) bind(c)
use prob_module, only : init_prob_data
integer, intent(in), value :: lev
real(amrex_real), intent(in), value :: time
type(c_ptr), intent(in), value :: pba, pdm

type(amrex_boxarray) :: ba
type(amrex_distromap) :: dm
type(amrex_mfiter) :: mfi
type(amrex_box) :: bx
real(amrex_real), contiguous, pointer :: phi(:,:,:,:)

ba = pba
dm = pdm

print*, "new_level_from_scratch",lev
call getCutCellDataStructures(amrex_geom(lev)%p, pba, pdm)

where my getCutCellDataStructures code looks like;

 void getCutCellDataStructures(amrex::Geometry& geom, const amrex::BoxArray& ba, const amrex::DistributionMapping& dm){

          amrex::Print() << " -- Build a fab array with embedded boundary data...\n";
           const amrex::EB2::IndexSpace& eb_is = amrex::EB2::IndexSpace::top();
           const amrex::EB2::Level& eb_level = eb_is.getLevel(geom);

           amrex::Print() << " -- set support...\n";
           amrex::EBSupport ebs = amrex::EBSupport::full;

           amrex::Print() << " -- nGhosts...\n";
           amrex::Vector<int> ng_ebs = {0, 0, 0};

           amrex::Print() << " -- Build factory to access EB database...\n";
           amrex::EBFArrayBoxFactory factory(eb_level, geom, ba, dm, ng_ebs, ebs);
           amrex::Print() << "done creating factory\n";

[WeiqunZhang]

Try to build in debug mode and see what you get.