diff --git a/src/coreComponents/physicsSolvers/wavePropagation/dg/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationDG.cpp b/src/coreComponents/physicsSolvers/wavePropagation/dg/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationDG.cpp index 770fce42c2..cff2374886 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/dg/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationDG.cpp +++ b/src/coreComponents/physicsSolvers/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/src/coreComponents/physicsSolvers/wavePropagation/dg/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationDGKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/dg/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationDGKernel.hpp index 21ec550625..705e687565 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/dg/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationDGKernel.hpp +++ b/src/coreComponents/physicsSolvers/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]; } );