Regarding FE calculation on Subdomain (example: a linearly elastic domain surrounded by PML)

[ac1512]

Dear MFEM users/ Developers,

I recently started using MFEM to develop a forward solver for linear elastodynamics with PML. Due to the introduction of PML, I have two different element attributes (one for the medium and the other for the PML). The PML formulation is done in such a way that it leads to the introduction of auxiliary variables in the PML subdomain. The formulation is a displacement formulation, and I am using H1 parallel finite element spaces. Displacement is an active variable throughout the domain. Thus, it leads to cross terms between the displacement and auxiliary PML variables. PML variables themselves have mass matrices that are restricted to the PML. As I went through previous threads (#1729), developers suggested using the concept of sub-mesh and defining Fespaces for the variables in the PML and then generating bilinear forms using MixedBilinearForm with different spaces as input for the cross terms. My question is:

The examples show the utility of the MixedBilinearForm for variables defined on the single mesh, but, in my case, the Fespaces are defined on different parts of the mesh, with one mesh being the parent mesh. So, I am a bit confused about whether this approach will still work. Especially during the assembly process, as there will be elements(outside the PML, inside the elastic medium) where few fespaces do not exist.

Another idea that I came up with is to use the same space for every variable (as I have the same regularity assumptions on all the variables in the system) but use the element attribute array as an input to AddDomainIntegrator() so that assembly happens on required subdomains and then use EliminateRowsCols() for the padding up the zero diagonals for deactivated nodes with 1.
In this case, I can generate a dofs list that I need to remove on each processor using the element attribute, but this leads to shared dofs at the partition interface, and EliminateRowsCols() tends not to work and goes into a segfault (Obviously my code can have bugs)?

Hence, I would appreciate any help figuring out the best way to overcome this problem.

Thanks,
Ankit Chakraborty

[tzanio]

ping: @psocratis, @cjvogl

[cjvogl]

Unfortunately, my experience with PML in MFEM doesn't go much further than knowing there's a example of PML for the indefinite Maxwell system (see example 25 in Examples)... a quick read of that example gives me the impression that the auxiliary variables are defined in a finite element space on the whole mesh, with the restriction to the boundary layer coming from the coefficients rather than the finite element basis

[psocratis]

sorry, I should have followed up on this. I met with @ac1512 and we discussed that one way to implement this is to construct transfer operators between the computational domain and the PML domain (without actually using submesh). We discuss how to do it and I think he is testing it. He needs this because he has variables living on the whole mesh and other variables living only in the PML region. So it's not as simple as ex25. @ac1512 let us know if what we discussed works.

Best,

Socratis

[ac1512]

@psocratis, Thanks for the discussion and explaining how to generate a restriction operator to restrict a Matrix or a vector to dofs corresponding to a submesh or subdomain with elements having a specific attribute. I am also attaching the a demo code given by @psocratis which performs this restriction on a vector from the mesh to subdomain and then also pulls back the restricted vector to the full mesh (as we need that for visualization using a fespace defined over the full mesh) for futre references.

#include "mfem.hpp"
#include <fstream>
#include <iostream>

using namespace std;
using namespace mfem;

#include <mpi.h>
#include <iostream>
#include <vector>
#include <map>
#include "mfem.hpp"

using namespace mfem;
using namespace std;

class ElementTdofs {
private:
   const ParFiniteElementSpace* pfes;
   std::vector<int> tdof_offsets;
   MPI_Comm comm;
   int num_procs;

   // Compute offsets for True DOFs across all processors
   void ComputeTdofOffsets() 
   {
      tdof_offsets.resize(num_procs);
      int mytoffset = pfes->GetMyTDofOffset();
      MPI_Allgather(&mytoffset, 1, MPI_INT, tdof_offsets.data(), 1, MPI_INT, comm);
   }

    // Determine which rank owns a given True DOF (tdof)
   int GetRank(int tdof) 
   {
      auto up = std::upper_bound(tdof_offsets.begin(), tdof_offsets.end(), tdof);
      return std::distance(tdof_offsets.begin(), up) - 1;
   }

   // Distribute indices using MPI_Alltoallv
   void DistributeIndices(Array<int>& indices, Array<int>& processors, Array<int>& recv_idx) 
   {
      // Step 1: Prepare data to send
      std::vector<int> send_counts(num_procs, 0);
      std::map<int, std::vector<int>> send_buffers;

      for (int i = 0; i < indices.Size(); ++i) 
      {
         send_buffers[processors[i]].push_back(indices[i]);
      }

      std::vector<int> send_displs(num_procs, 0);
      std::vector<int> send_data;

      for (int i = 0; i < num_procs; ++i) 
      {
         send_counts[i] = send_buffers[i].size();
         send_data.insert(send_data.end(), send_buffers[i].begin(), send_buffers[i].end());
      }

      for (int i = 1; i < num_procs; ++i) 
      {
         send_displs[i] = send_displs[i - 1] + send_counts[i - 1];
      }

      // Step 2: Gather receive counts
      std::vector<int> recv_counts(num_procs, 0);
      MPI_Alltoall(send_counts.data(), 1, MPI_INT, recv_counts.data(), 1, MPI_INT, comm);

      // Step 3: Compute receive displacements
      std::vector<int> recv_displs(num_procs, 0);
      int total_recv_size = 0;

      for (int i = 0; i < num_procs; ++i) 
      {
         recv_displs[i] = total_recv_size;
         total_recv_size += recv_counts[i];
      }

      // Step 4: Allocate buffer and perform MPI_Alltoallv
      recv_idx.SetSize(total_recv_size);
      MPI_Alltoallv(send_data.data(), send_counts.data(), send_displs.data(), MPI_INT,
                    recv_idx.GetData(), recv_counts.data(), recv_displs.data(), MPI_INT, comm);
   }

public:
   ElementTdofs(const ParFiniteElementSpace* pfes_) 
        : pfes(pfes_), comm(pfes_->GetComm()) 
   {
      MPI_Comm_size(comm, &num_procs);
      ComputeTdofOffsets();
   }

    // Extract and distribute True DOFs based on the element attribute
   Array<int> GetTrueDOFs(int element_attribute) 
   {
      Array<int> mydof_list, other_dof_list;

      for (int i = 0; i < pfes->GetNE(); i++) 
      {
         if (pfes->GetParMesh()->GetAttribute(i) == element_attribute) 
         {
            Array<int> dofs;
            pfes->GetElementVDofs(i, dofs);

            for (int j = 0; j < dofs.Size(); j++) 
            {
               int tdof = pfes->GetGlobalTDofNumber(dofs[j]);
               if (tdof >= pfes->GetMyTDofOffset() && 
                   tdof < pfes->GetMyTDofOffset() + pfes->GetTrueVSize()) 
               {
                  mydof_list.Append(tdof);
               } 
               else 
               {
                  other_dof_list.Append(tdof);
               }
            }
         }
      }

      mydof_list.Sort();
      mydof_list.Unique();
      other_dof_list.Sort();
      other_dof_list.Unique();

      // Assign processors for other DOFs
      Array<int> other_dof_proc_list(other_dof_list.Size());
      for (int i = 0; i < other_dof_list.Size(); i++) 
      {
         other_dof_proc_list[i] = GetRank(other_dof_list[i]);
      }

      // Communicate DOFs to appropriate processors
      Array<int> recv_idx;
      DistributeIndices(other_dof_list, other_dof_proc_list, recv_idx);

      // Merge received DOFs into the local list
      mydof_list.Append(recv_idx);
      mydof_list.Sort();
      mydof_list.Unique();

      return mydof_list;
    }
};



real_t test_prolongation(const Vector &x)
{
   return sin(M_PI*(x(0)+x(1))/100);
}

int main(int argc, char *argv[])
{
   // 1. Initialize MPI and HYPRE.
   Mpi::Init(argc, argv);
   int num_procs = Mpi::WorldSize();
   int myid = Mpi::WorldRank();
   Hypre::Init();

   // 2. Parse command line options.
   string mesh_file = "PML_case_b.msh";
   int order = 1;
   int par_ref_levels = 0;

   OptionsParser args(argc, argv);
   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
   args.AddOption(&order, "-o", "--order", "Finite element polynomial degree");
   args.AddOption(&par_ref_levels, "-pr", "--par_ref_levels", "number of parallel refiniement levels");
   args.Parse();
   if (!args.Good())
   {
      if (myid == 0)
      {
         args.PrintUsage(cout);
      }
      return 1;
   }
   if (myid == 0)
   {
      args.PrintOptions(cout);
   }

   Mesh mesh(mesh_file);
   
   ParMesh pmesh(MPI_COMM_WORLD, mesh);
   mesh.Clear(); 
   for (int l = 0; l < par_ref_levels; l++)
   {
      pmesh.UniformRefinement();
   }

   H1_FECollection fec(order, pmesh.Dimension());
   ParFiniteElementSpace fespace(&pmesh, &fec);

   ElementTdofs etdofs1(&fespace);

   Array<int> dof_list1 = etdofs1.GetTrueDOFs(1);
   Array<int> dof_list2 = etdofs1.GetTrueDOFs(2);

   int h1 = dof_list1.Size();
   SparseMatrix St1(h1,fespace.GlobalTrueVSize());

   for (int i = 0; i<h1; i++)
   {
      int col = dof_list1[i];
      St1.Set(i,col,1.0);
   }
   St1.Finalize();
   int rows1[2];
   int cols1[2];
   int nrows1 = St1.Height();

   int row_offset1;
   MPI_Scan(&nrows1,&row_offset1,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);

   row_offset1-=nrows1;
   rows1[0] = row_offset1;
   rows1[1] = row_offset1+nrows1;
   for (int i = 0; i < 2; i++)
   {
      cols1[i] = fespace.GetTrueDofOffsets()[i];
   }
   int glob_nrows1;
   int glob_ncols1 = fespace.GlobalTrueVSize();
   MPI_Allreduce(&nrows1, &glob_nrows1,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);

   mfem::out << "global nrows1 = " << glob_nrows1 << endl;

   HypreParMatrix * P1t = new HypreParMatrix(MPI_COMM_WORLD, nrows1, glob_nrows1,
                                             glob_ncols1, St1.GetI(), St1.GetJ(),
                                             St1.GetData(), rows1,cols1); 
   HypreParMatrix * P1 = P1t->Transpose();


   int h2 = dof_list2.Size();
   SparseMatrix St2(h2,fespace.GlobalTrueVSize());

   for (int i = 0; i<h2; i++)
   {
      int col = dof_list2[i];
      St2.Set(i,col,1.0);
   }
   St2.Finalize();
   int rows2[2];
   int cols2[2];
   int nrows2 = St2.Height();

   int row_offset2;
   MPI_Scan(&nrows2,&row_offset2,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);

   row_offset2-=nrows2;
   rows2[0] = row_offset2;
   rows2[1] = row_offset2+nrows2;
   for (int i = 0; i < 2; i++)
   {
      cols2[i] = fespace.GetTrueDofOffsets()[i];
   }
   int glob_nrows2;
   int glob_ncols2 = fespace.GlobalTrueVSize();
   MPI_Allreduce(&nrows2, &glob_nrows2,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);

   mfem::out << "global nrows2 = " << glob_nrows2 << endl;


   HypreParMatrix * P2t = new HypreParMatrix(MPI_COMM_WORLD, nrows2, glob_nrows2,
                                             glob_ncols2, St2.GetI(), St2.GetJ(),
                                             St2.GetData(), rows2,cols2); 
   HypreParMatrix * P2 = P2t->Transpose();


   ParGridFunction gf(&fespace);
   FunctionCoefficient cf(test_prolongation);
   gf.ProjectCoefficient(cf);

   Vector tv(fespace.GetTrueVSize());
   gf.GetTrueDofs(tv);

   Vector restricted1_tv(P1->Width());
   P1->MultTranspose(tv, restricted1_tv);
   tv = 0.0;
   P1->Mult(restricted1_tv, tv);
   ParGridFunction restricted1_gf(&fespace);
   restricted1_gf.SetFromTrueDofs(tv);
   
   // Visualize both gf and restricted gf in GLvis
   if (1)
   {
      char vishost[] = "localhost";
      int  visport   = 19916;
      socketstream sol_sock(vishost, visport);
      sol_sock << "parallel " << num_procs << " " << myid << "\n";
      sol_sock.precision(8);
      sol_sock << "solution\n" << pmesh << gf << flush;
   }

   if (1)
   {
      char vishost[] = "localhost";
      int  visport   = 19916;
      socketstream sol_sock_r(vishost, visport);
      sol_sock_r << "parallel " << num_procs << " " << myid << "\n";
      sol_sock_r.precision(8);
      sol_sock_r << "solution\n" << pmesh << restricted1_gf 
                 << "window_title 'Domain 1'" << flush;

   }

   gf.GetTrueDofs(tv);
   Vector restricted2_tv(P2->Width());
   P2->MultTranspose(tv, restricted2_tv);
   tv = 0.0;
   P2->Mult(restricted2_tv, tv);
   ParGridFunction restricted2_gf(&fespace);
   restricted2_gf.SetFromTrueDofs(tv);

   if (1)
   {
      char vishost[] = "localhost";
      int  visport   = 19916;
      socketstream sol_sock_r(vishost, visport);
      sol_sock_r << "parallel " << num_procs << " " << myid << "\n";
      sol_sock_r.precision(8);
      sol_sock_r << "solution\n" << pmesh << restricted2_gf 
                 << "window_title 'Domain 2'" << flush;
   }

   return 0;
}