Bugfix

GEOS-DEV · Jan 30, 2025 · acfcb44 · acfcb44
1 parent 252da89
commit acfcb44
Show file tree

Hide file tree

Showing 2 changed files with 25 additions and 22 deletions.
diff --git a/...csSolvers/wavePropagation/dg/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationDG.cpp b/...csSolvers/wavePropagation/dg/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationDG.cpp
@@ -518,6 +518,8 @@ void AcousticWaveEquationDG::computeUnknowns( real64 const & time_n,
           using FE_TYPE = TYPEOFREF( finiteElement );
 
 
+          printf("before call");
+
           AcousticWaveEquationDGKernels::
           PressureComputationKernel::
           pressureComputation< FE_TYPE, EXEC_POLICY, ATOMIC_POLICY >
@@ -559,27 +561,28 @@ void AcousticWaveEquationDG::synchronizeUnknowns( real64 const & time_n,
                                                    real64 const & dt,
                                                    DomainPartition & domain,
                                                    MeshLevel & mesh,
-                                                   arrayView1d< string const > const & )
+                                                   arrayView1d< string const > const & regionNames)
 {
   NodeManager & nodeManager = mesh.getNodeManager();
 
-  arrayView2d< real32 > const p_n = nodeManager.getField< acousticfieldsdg::Pressure_n >();
-  arrayView2d< real32 > const p_np1 = nodeManager.getField< acousticfieldsdg::Pressure_np1 >();
-
-  // arrayView2d< real32 > const stiffnessVector = nodeManager.getField< acousticfieldsdg::StiffnessVector >();
-  //arrayView1d< real32 > const rhs = nodeManager.getField< acousticfieldsdg::ForcingRHS >();
-
-  /// synchronize pressure fields
-  FieldIdentifiers fieldsToBeSync;
-  fieldsToBeSync.addFields( FieldLocation::Node, { acousticfieldsdg::Pressure_np1::key() } );
+  mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const regionIndex,
+                                                                                          CellElementSubRegion & elementSubRegion )
+  {
+
+    FieldIdentifiers fieldsToBeSync;
+    fieldsToBeSync.addElementFields( {acousticfieldsdg::Pressure_nm1::key(), acousticfieldsdg::Pressure_n::key(), acousticfieldsdg::Pressure_np1::key()}, regionNames );
 
-  CommunicationTools & syncFields = CommunicationTools::getInstance();
-  syncFields.synchronizeFields( fieldsToBeSync,
+    CommunicationTools & syncFields = CommunicationTools::getInstance();
+    syncFields.synchronizeFields( fieldsToBeSync,
                                 mesh,
                                 domain.getNeighbors(),
                                 true );
+
+
+  } );
+
   /// compute the seismic traces since last step.
-  arrayView2d< real32 > const pReceivers = m_pressureNp1AtReceivers.toView();
+  //arrayView2d< real32 > const pReceivers = m_pressureNp1AtReceivers.toView();
 
   //computeAllSeismoTraces( time_n, dt, p_np1, p_n, pReceivers );
   //incrementIndexSeismoTrace( time_n );

diff --git a/...ers/wavePropagation/dg/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationDGKernel.hpp b/...ers/wavePropagation/dg/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationDGKernel.hpp
@@ -451,16 +451,14 @@ struct PressureComputationKernel
 
  {
 
+
     real64 const rickerValue = useSourceWaveletTables ? 0 : WaveSolverUtils::evaluateRicker( time_n, timeSourceFrequency, timeSourceDelay, rickerOrder );
 
-     printf("herebefore");
 
     //For now lots of comments with ideas  + needed array to add to the method prototype
     forAll< EXEC_POLICY >( 1, [=] GEOS_HOST_DEVICE ( localIndex const k )
     {
 
-      printf("here\n");
-
       real64 const dt2 = pow(dt,2);
 
       constexpr localIndex numNodesPerElem = FE_TYPE::numNodes;
@@ -483,6 +481,7 @@ struct PressureComputationKernel
       //Multiply by p_{n } by 2*Mass
       FE_TYPE::computeMassTerm(xLocal, [&] (const int i, const int j, const real64 val)
       {
+
          flowx[j] = 2.0*val*p_n[k][i];
          flowx[j] -= val*p_nm1[k][i];
       } );
@@ -493,6 +492,7 @@ struct PressureComputationKernel
       //First stiffness part (volume)
       FE_TYPE::computeStiffnessTerm(xLocal, [&] (const int i, const int j, real64 val)
       {
+
          flowx[j] -= dt2*val*p_n[k][i];
       } );
 
@@ -527,8 +527,8 @@ struct PressureComputationKernel
 
          //Flux computation
 
-         flowx[c1] -= 0.5*dt2*p_n[k][c2];
-         flowx[c1] += 0.5*dt2*p_n[elemNeigh][neighDof];
+         flowx[c1] -= 0.5*dt2*val*p_n[k][c2];
+         flowx[c1] += 0.5*dt2*val*p_n[elemNeigh][neighDof];
 
 
        },
@@ -562,8 +562,8 @@ struct PressureComputationKernel
 
          //Flux computation
 
-         flowx[c1] += 0.5*dt2*p_n[elemNeigh][neighDof];
-         flowx[c1] -= 0.5*dt2*p_n[k][c2];
+         flowx[c1] += 0.5*dt2*val*p_n[elemNeigh][neighDof];
+         flowx[c1] -= 0.5*dt2*val*p_n[k][c2];
 
          //Then we need a second time where we take the transpose of the previous values:
 
@@ -580,8 +580,8 @@ struct PressureComputationKernel
 
          //Flux computation
 
-         flowx[c2] -= 0.5*dt2*p_n[elemNeigh][neighDof2];
-         flowx[c2] += 0.5*dt2*p_n[k][c1];
+         flowx[c2] -= 0.5*dt2*val*p_n[elemNeigh][neighDof2];
+         flowx[c2] += 0.5*dt2*val*p_n[k][c1];
 
 
        } );