Mfem: Calculating jumps across elements of a GridFunction in parallel

[salazardetroya]

I have an L2 function of order 0. I want to calculate the gradient by looking at the jump of a GridFunction across the element faces. So I iterate over the faces (pmesh->GetNumFaces();) and at each face I calculate the jump (difference in value) of a GridFunction. Then I add that jump to a piecewise constant GridFunction

This works in serial, but it does not in parallel (I obtain nan values). How could I make it work in parallel?

The function in question:

std::shared_ptr<ParGridFunction> geometry_gradient(std::shared_ptr<ParGridFunction> geometry,
        std::shared_ptr<ParFiniteElementSpace> geom_grad_fespace) {

    ParFiniteElementSpace *fes = geometry->ParFESpace();
    ParMesh * pmesh = fes->GetParMesh();
    const unsigned int dim = pmesh->Dimension();
    Array<int> vdofs1, vdofs2;
    Vector rho1, rho2;

    auto geometry_gradient = std::make_shared<ParGridFunction>(geom_grad_fespace.get());

    const int num_faces = pmesh->GetNumFaces();
    for (int i = 0; i < num_faces; i++) {
        FaceElementTransformations *FTr = pmesh->GetInteriorFaceTransformations(i);
        if (FTr == NULL) {
            continue;
        }

        fes->GetElementVDofs(FTr->Elem1No, vdofs1);
        fes->GetElementVDofs(FTr->Elem2No, vdofs2);
        geometry->GetSubVector(vdofs1, rho1);
        geometry->GetSubVector(vdofs2, rho2);


        Vector center1(2), center2(2);;
        pmesh->GetElementCenter(FTr->Elem1No, center1);
        pmesh->GetElementCenter(FTr->Elem2No, center2);
        double distance = center1.Add(-1.0, center2).Norml2();

        double face_error = (rho1(0) - rho2(0)) / distance;
        face_error = std::abs(face_error);

        geometry_gradient->operator()(FTr->Elem1No) += 0.5 * face_error;
        geometry_gradient->operator()(FTr->Elem2No) += 0.5 * face_error;
    }

    return geometry_gradient;
}

Whole code to reproduce

#include "mfem.hpp"
#include <memory>
using namespace mfem;

class CrossSectionCoefficient : public Coefficient {
public:
    CrossSectionCoefficient(double epsilon, double bar_length) : epsilon_(epsilon), bar_length_(bar_length) {}

    virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip) {
        double x[2];
        Vector transip(x, 2);
        T.Transform(ip, transip);

        return heaviside(transip[0]);
    }

    double heaviside(double distance) {
        return (1.0/2.0*(1.0 + tanh(epsilon_ * (distance - bar_length_ / 2.0)) / (tanh(epsilon_ * (bar_length_ / 2.0)))));
    }

private:
    double epsilon_, bar_length_;
};

std::shared_ptr<ParGridFunction> geometry_gradient(std::shared_ptr<ParGridFunction> geometry,
        std::shared_ptr<ParFiniteElementSpace> geom_grad_fespace) {

    ParFiniteElementSpace *fes = geometry->ParFESpace();
    ParMesh * pmesh = fes->GetParMesh();
    const unsigned int dim = pmesh->Dimension();
    Array<int> vdofs1, vdofs2;
    Vector rho1, rho2;

    auto geometry_gradient = std::make_shared<ParGridFunction>(geom_grad_fespace.get());

    const int num_faces = pmesh->GetNumFaces();
    for (int i = 0; i < num_faces; i++) {
        FaceElementTransformations *FTr = pmesh->GetInteriorFaceTransformations(i);
        if (FTr == NULL) {
            continue;
        }

        fes->GetElementVDofs(FTr->Elem1No, vdofs1);
        fes->GetElementVDofs(FTr->Elem2No, vdofs2);
        geometry->GetSubVector(vdofs1, rho1);
        geometry->GetSubVector(vdofs2, rho2);


        Vector center1(2), center2(2);;
        pmesh->GetElementCenter(FTr->Elem1No, center1);
        pmesh->GetElementCenter(FTr->Elem2No, center2);
        double distance = center1.Add(-1.0, center2).Norml2();

        double face_error = (rho1(0) - rho2(0)) / distance;
        face_error = std::abs(face_error);

        geometry_gradient->operator()(FTr->Elem1No) += 0.5 * face_error;
        geometry_gradient->operator()(FTr->Elem2No) += 0.5 * face_error;
    }

    return geometry_gradient;
}

int main(int argc, char *argv[])
{
    int num_procs, myid;
    MPI_Init(&argc, &argv);
    MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
    MPI_Comm_rank(MPI_COMM_WORLD, &myid);

    const char *mesh_file = "beam-quad.mesh";
    OptionsParser args(argc, argv);
    int ref_levels = 5;
    double heaviside_epsilon = 200;
    args.AddOption(&ref_levels, "-rs", "--refine-serial",
                   "Number of times to refine the mesh uniformly in serial,"
                   " -1 for auto.");
    args.AddOption(&heaviside_epsilon, "-e", "--epsilon",
                   "Heaviside approximation parameter."
                   " Negative values are with a mesh dependent heuristic.");
    args.Parse();

    Mesh *mesh = new Mesh(mesh_file, 1, 1);
    int dim = mesh->Dimension();
    for (int l = 0; l < ref_levels; l++) {
        mesh->UniformRefinement();
    }

    ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
    delete mesh;

    auto l2_fec = std::make_shared<L2_FECollection>(0, dim);
    auto l2_fespace =
        std::make_shared<ParFiniteElementSpace>(pmesh, l2_fec.get());
    auto l2_fespace_flux =
        std::make_shared<ParFiniteElementSpace>(pmesh, l2_fec.get(), dim);

    CrossSectionCoefficient cross_section(heaviside_epsilon, 8.0);
    auto geometry = std::make_shared<ParGridFunction>(l2_fespace.get());
    geometry->ProjectCoefficient(cross_section);
    auto geometry_gradient_func = geometry_gradient(geometry, l2_fespace_flux);


    ParaViewDataCollection paraview_dc_bf("test_boundary", pmesh);
    paraview_dc_bf.SetLevelsOfDetail(1);
    paraview_dc_bf.SetCycle(1);
    paraview_dc_bf.SetTime(0.0);
    paraview_dc_bf.RegisterField("Geometry gradient", geometry_gradient_func.get());
    paraview_dc_bf.RegisterField("Geometry ", geometry.get());
    paraview_dc_bf.Save();

    delete pmesh;
}

[salazardetroya]

Closing this for now. I need to give a better thought to the method I'm using here.

[pazner]

Closing this for now. I need to give a better thought to the method I'm using here.

I hope you won't mind if I chime in on this anyway.

It is relatively straightforward to get your code to work in parallel. First, you need to make sure that the "face neighbor data" is communicated in parallel. You can do this by calling geometry->ExchangeFaceNbrData() before you begin your computations.

The loop you currently have over the faces will skip any face that is on a processor boundary. This is because GetInteriorFaceTransformations will return NULL for such a face. You need to add a second loop over the faces that are shared between two MPI processes. Something like:

for (int i = 0; i < pmesh->GetNSharedFaces(); i++)

and then use GetSharedFaceTransformations to get the FaceElementTransformations object. In this case, FTr->Elem2No refers to a "face neighbor element", and you need to get its vdofs and extract the grid function subvector like this:

fes->GetFaceNbrElementVDofs(FTr->Elem2No, vdofs2);
geometry->FaceNbrData().GetSubVector(vdofs2, rho2);

You'll also have to make your own version of GetElementCenter to make it work with face neighbor elements. This is easy to do using the ElementTransformation.
In this second loop, you don't need to modify geometry_gradient corresponding to the neighboring element, since the MPI process who owns that element will take care of that operation.

Hope this helps...

[salazardetroya]

Thanks @pazner. I ended up using a different method to calculate the gradient. This method I showed here is actually not accurate. In any case, it is good to know how to deal with shared faces so I appreciate your comment.