Generalizing boundary condition interface

[tomidakn]

We now have many physics modules besides MHD, and I am planning to merge Pinghui Huang's dust field module. However, currently we do not have a general interface for (physical) boundary conditions. For example, the passive scalar module has no interface to set its boundary conditions. As some physics modules are coupled to each other (e.g., we may want to change the passive scalar according to the flow direction), it is not very practical to have separate interfaces for each physics. However, it is also not practical to keep adding new physics to the boundary condition functions, as it cannot maintain backward compatibility.

So, I would like to propose a new interface to register variables for boundary conditions using std::vector, which should look like:

int bcidx_scalar;

Mesh::InitUserMeshdata(ParameterInput *pin) {
  EnrollUserBoundaryFunction(BoundaryFace::inner_x1, MyBoundary_ix1);
}

MeshBlock::InitUserMeshBlockData(ParameterInput *pin) {
  bcidx_scalar = RegisterBoundaryVariable(&(pscalar->r), &(pscalar->coarse_r_));
}

void MyBoundary_ix1(MeshBlock *pmb, Coordinates *pco,
                    AthenaArray<Real> &prim, FaceField &b,
                    std::vector<AthenaArray<Real>*> bcv,
                    Real time, Real dt,
                    int il, int iu, int jl, int ju, int kl, int ku, int ngh) {
  AthenaArray<Real> &r = *(bcv[bcidx_scalar]);
  for (int k=kl; k<=ku; ++k) {
    for (int j=jl; j<=ju; ++j) {
      for (int i=1; i<=ngh; ++i) {
        prim(IDN,k,j,il-i) = prim(IDN,k,j,il);
        prim(IVX,k,j,il-i) = prim(IVX,k,j,il);
        prim(IVY,k,j,il-i) = prim(IVY,k,j,il);
        prim(IVZ,k,j,il-i) = prim(IVZ,k,j,il);
        if (NON_BAROTROPIC_EOS)
          prim(IPR,k,j,il-i) = prim(IPR,k,j,il);
        for (int s = 0; s < NSCALARS; ++s)
          r(s,k,j,il-i) = r(s,k,j,il);
      }
    }
  }
}

I admit it is a bit awkward to register it in InitUserMeshBlockData, but I could not find any other way.

In fact, something similar is already implemented for mesh refinement and boundary communications, so we may integrate them, although there are some differences (e.g., physical BC is optional, but communications are necessary).

Note that this is for physical processes included in the main time integrator. Physics implemented as a separate operator (e.g. Multigrid) should have its own boundary condition interface as is.

How do you think?