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 986b1ec279..cff2374886 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/dg/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationDG.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/dg/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationDG.cpp @@ -79,7 +79,7 @@ void AcousticWaveEquationDG::initializePreSubGroups() void AcousticWaveEquationDG::registerDataOnMesh( Group & meshBodies ) { - + WaveSolverBase::registerDataOnMesh( meshBodies ); forDiscretizationOnMeshTargets( meshBodies, [&] ( string const &, MeshLevel & mesh, arrayView1d< string const > const & ) @@ -221,7 +221,7 @@ void AcousticWaveEquationDG::initializePostInitialConditionsPreSubGroups() DomainPartition & domain = getGroupByPath< DomainPartition >( "/Problem/domain" ); - applyFreeSurfaceBC( 0.0, domain ); + //applyFreeSurfaceBC( 0.0, domain ); forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const & meshBodyName, MeshLevel & mesh, @@ -246,7 +246,7 @@ void AcousticWaveEquationDG::initializePostInitialConditionsPreSubGroups() /// get array of indicators: 1 if face is on the free surface; 0 otherwise arrayView1d< localIndex const > const freeSurfaceFaceIndicator = faceManager.getField< acousticfieldsdg::AcousticFreeSurfaceFaceIndicator >(); - + m_referenceInvMassMatrix.resize( elemManager.numRegions() ); m_boundaryInvMassPlusDamping.resize( elemManager.numRegions() ); @@ -269,7 +269,7 @@ void AcousticWaveEquationDG::initializePostInitialConditionsPreSubGroups() "Invalid type of element, the acoustic DG solver is designed for tetrahedral meshes only ", InputError ); - + finiteElement::FiniteElementBase const & fe = elementSubRegion.getReference< finiteElement::FiniteElementBase >( getDiscretizationName() ); @@ -289,6 +289,7 @@ void AcousticWaveEquationDG::initializePostInitialConditionsPreSubGroups() /// Partial gradient if gradient as to be computed + finiteElement::FiniteElementDispatchHandler< DG_FE_TYPES >::dispatch3D( fe, [&] ( auto const finiteElement ) { using FE_TYPE = TYPEOFREF( finiteElement ); @@ -306,26 +307,32 @@ void AcousticWaveEquationDG::initializePostInitialConditionsPreSubGroups() freeSurfaceFaceIndicator, elemsToOpposite, elemsToOppositePermutation ); - - AcousticWaveEquationDGKernels:: - PrecomputePenaltyGeomKernel:: - launch< EXEC_POLICY, FE_TYPE > - ( elementSubRegion.size(), - nodeCoords, - elemsToNodes, - characteristicSize ); - // Pre-compute inverse of mass + damping matrix for each boundary element and penalty coefficient parameter - AcousticWaveEquationDGKernels:: - PrecomputeMassDampingKernel:: - launch< EXEC_POLICY, ATOMIC_POLICY, FE_TYPE > - ( elementSubRegion.size(), - nodeCoords, - elemsToNodes, - elemsToOpposite, - m_referenceInvMassMatrix[ regionIndex ][ subRegionIndex ], - m_boundaryInvMassPlusDamping[ regionIndex ][ subRegionIndex ], - m_timeStep ); + // Precompute reference mass matrix for non-boundary elements + m_referenceInvMassMatrix[ regionIndex ][ subRegionIndex ].resizeDimension< 0, 1 >( FE_TYPE::numNodes, FE_TYPE::numNodes ); + array2d< real64 > massMatrix; + massMatrix.resize( FE_TYPE::numNodes, FE_TYPE::numNodes ); + massMatrix.zero(); + FE_TYPE::computeReferenceMassMatrix( massMatrix ); + BlasLapackLA::matrixInverse( massMatrix, m_referenceInvMassMatrix[ regionIndex ][ subRegionIndex ] ); + + // Pre-compute inverse of mass + damping matrix for each boundary element + // localIndex nAbsBdryElems = 0; + // forAll< EXEC_POLICY >( elementSubRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) + // { + // characteristicSize[ k ] = WaveSolverUtils::computeReferenceLengthForPenalty( elemsToNodes, nodeCoords, k ); + // bool bdry = false; + // for( int i = 0; i < 4; i++ ) + // { + // if( elemsToOpposite( k, i ) == -1 ) + // { + // bdry = true; + // break; + // } + // } + // RAJA::atomicInc< ATOMIC_POLICY >( &m_indexToBoundaryMatrix[ k ]) + // } ); + // m_boundaryInvMassPlusDamping[ regionIndex ][ subRegionIndex ].resizeDimension< 0, 1, 2 >( nAbsBdryElems, FE_TYPE::numNodes, FE_TYPE::numNodes ); // AcousticMatricesSEM::MassMatrix< FE_TYPE > kernelM( finiteElement ); // kernelM.template computeMassMatrix< EXEC_POLICY, ATOMIC_POLICY >( elementSubRegion.size(), @@ -378,7 +385,7 @@ void AcousticWaveEquationDG::applyFreeSurfaceBC( real64 const time, DomainPartit FaceManager & faceManager = domain.getMeshBody( 0 ).getMeshLevel( m_discretizationName ).getFaceManager(); NodeManager & nodeManager = domain.getMeshBody( 0 ).getMeshLevel( m_discretizationName ).getNodeManager(); - arrayView2d< real32 > const p_np1 = nodeManager.getField< acousticfieldsdg::Pressure_np1 >(); + // arrayView2d< real32 > const p_np1 = nodeManager.getField< acousticfieldsdg::Pressure_np1 >(); ArrayOfArraysView< localIndex const > const faceToNodeMap = faceManager.nodeList().toViewConst(); @@ -450,30 +457,30 @@ real64 AcousticWaveEquationDG::explicitStepBackward( real64 const & time_n, } -void AcousticWaveEquationDG::prepareNextTimestep( MeshLevel & mesh ) -{ - NodeManager & nodeManager = mesh.getNodeManager(); - +void AcousticWaveEquationDG::prepareNextTimestep( MeshLevel & mesh ) +{ + NodeManager & nodeManager = mesh.getNodeManager(); + arrayView2d< real32 > const p_nm1 = nodeManager.getField< acousticfieldsdg::Pressure_nm1 >(); 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 >(); //arrayView2d< real32 > const rhs = nodeManager.getField< acousticfieldsdgdgdg::ForcingRHS >(); - + SortedArrayView< localIndex const > const solverTargetNodesSet = m_solverTargetNodesSet.toViewConst(); - + // forAll< EXEC_POLICY >( solverTargetNodesSet.size(), [=] GEOS_HOST_DEVICE ( localIndex const n ) - // { - // localIndex const a = solverTargetNodesSet[n]; - - // p_nm1[a] = p_n[a]; - // p_n[a] = p_np1[a]; - - // stiffnessVector[a] = rhs[a] = 0.0; - // } ); + // { + // localIndex const a = solverTargetNodesSet[n]; + + // p_nm1[a] = p_n[a]; + // p_n[a] = p_np1[a]; -} + // stiffnessVector[a] = rhs[a] = 0.0; + // } ); + +} void AcousticWaveEquationDG::computeUnknowns( real64 const & time_n, real64 const & dt, @@ -500,17 +507,19 @@ void AcousticWaveEquationDG::computeUnknowns( real64 const & time_n, arrayView2d< localIndex > const elemsToOpposite = elementSubRegion.getField< acousticfieldsdg::ElementToOpposite >(); arrayView2d< integer > const elemsToOppositePermutation = elementSubRegion.getField< acousticfieldsdg::ElementToOppositePermutation >(); - arrayView2d< real32 const > const p_nm1 = nodeManager.getField< acousticfieldsdg::Pressure_nm1 >(); - arrayView2d< real32 const > const p_n = nodeManager.getField< acousticfieldsdg::Pressure_n >(); - arrayView2d< real32 > const p_np1 = nodeManager.getField< acousticfieldsdg::Pressure_np1 >(); + arrayView2d< real32 const > const p_nm1 = elementSubRegion.getField< acousticfieldsdg::Pressure_nm1 >(); + arrayView2d< real32 const > const p_n = elementSubRegion.getField< acousticfieldsdg::Pressure_n >(); + arrayView2d< real32 > const p_np1 = elementSubRegion.getField< acousticfieldsdg::Pressure_np1 >(); finiteElement::FiniteElementBase const & fe = elementSubRegion.getReference< finiteElement::FiniteElementBase >( getDiscretizationName() ); finiteElement::FiniteElementDispatchHandler< DG_FE_TYPES >::dispatch3D( fe, [&] ( auto const finiteElement ) { using FE_TYPE = TYPEOFREF( finiteElement ); - - + + + printf("before call"); + AcousticWaveEquationDGKernels:: PressureComputationKernel:: pressureComputation< FE_TYPE, EXEC_POLICY, ATOMIC_POLICY > @@ -536,7 +545,7 @@ void AcousticWaveEquationDG::computeUnknowns( real64 const & time_n, p_np1 ); } ); - + } ); // /// calculate your time integrators @@ -552,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 0a8ae3b56e..705e687565 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/dg/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationDGKernel.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/dg/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationDGKernel.hpp @@ -220,7 +220,6 @@ struct PrecomputeNeighborhoodKernel forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k1 ) { localIndex vertices[ 4 ] = { elemsToNodes( k1, 0 ), elemsToNodes( k1, 1 ), elemsToNodes( k1, 2 ), elemsToNodes( k1, 3 ) }; - for( int i = 0; i < 4; i++ ) { localIndex k1OrderedVertices[ 3 ]; @@ -234,8 +233,11 @@ struct PrecomputeNeighborhoodKernel } // find opposite vertex in first element int o1 = -1; + int indexo1= -1; + int vertex = -1; int count = 0; - for ( localIndex vertex : vertices) { + for ( localIndex k=0; k< 4; ++k) { + vertex = vertices[k]; bool found = false; for ( int j = 0; j < 3; j++ ) { @@ -248,27 +250,31 @@ struct PrecomputeNeighborhoodKernel if( !found ) { o1 = vertex; + indexo1=k; } else { k1OrderedVertices[ count++ ] = vertex; } } + GEOS_ERROR_IF( o1 < 0, "Topological error in mesh: a face and its adjacent element share all vertices."); if( k2 < 0 ) { // boundary element, either free surface, or absorbing boundary - elemsToOpposite( k1, o1 ) = freeSurfaceFaceIndicator( f ) == 1 ? -2 : -1; - elemsToOppositePermutation( k1, o1 ) = 0; + elemsToOpposite( k1, indexo1 ) = freeSurfaceFaceIndicator( f ) == 1 ? -2 : -1; + elemsToOppositePermutation( k1, indexo1 ) = 0; } else { - elemsToOpposite( k1, o1 ) = k2; + elemsToOpposite( k1, indexo1 ) = k2; localIndex oppositeElemVertices[ 4 ] = { elemsToNodes( k2, 0 ), elemsToNodes( k2, 1 ), elemsToNodes( k2, 2 ), elemsToNodes( k2, 3 ) }; // find opposite vertex in second element int o2 = -1; + int indexo2 = -1; count = 0; - for ( localIndex vertex : oppositeElemVertices) { + for ( localIndex k=0; k<4; ++k) { + vertex = vertices[k]; bool found = false; for ( int j = 0; j < 3; j++ ) { @@ -281,7 +287,9 @@ struct PrecomputeNeighborhoodKernel if( !found ) { o2 = vertex; + indexo2 = k; } + } GEOS_ERROR_IF( o2 < 0, "Topological error in mesh: a face and its adjacent element share all vertices."); // compute permutation @@ -301,7 +309,7 @@ struct PrecomputeNeighborhoodKernel permutation = permutation + c * ( position + 1 ); c = c * 4; } - elemsToOppositePermutation( k1, o1 ) = permutation; + elemsToOppositePermutation( k1, indexo1 ) = permutation; } } } ); @@ -364,7 +372,7 @@ struct PrecomputeMassDampingKernel { // Precompute reference mass matrix for non-boundary elements referenceInvMassMatrix.resizeDimension< 0, 1 >( FE_TYPE::numNodes, FE_TYPE::numNodes ); - array2d< real64 > massMatrix; + array2d< real64 > massMatrix; massMatrix.resize( FE_TYPE::numNodes, FE_TYPE::numNodes ); massMatrix.zero(); FE_TYPE::computeReferenceMassMatrix( massMatrix ); @@ -387,7 +395,7 @@ struct PrecomputeMassDampingKernel RAJA::atomicInc< ATOMIC_POLICY >( &nAbsBdryElems ); } } ); - + boundaryInvMassPlusDamping.resizeDimension< 0, 1, 2 >( nAbsBdryElems, FE_TYPE::numNodes, FE_TYPE::numNodes ); forAll< EXEC_POLICY >( size, [&] ( localIndex const k ) { @@ -443,10 +451,12 @@ struct PressureComputationKernel { + real64 const rickerValue = useSourceWaveletTables ? 0 : WaveSolverUtils::evaluateRicker( time_n, timeSourceFrequency, timeSourceDelay, rickerOrder ); + //For now lots of comments with ideas + needed array to add to the method prototype - forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k ) + forAll< EXEC_POLICY >( 1, [=] GEOS_HOST_DEVICE ( localIndex const k ) { real64 const dt2 = pow(dt,2); @@ -471,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]; } ); @@ -481,10 +492,15 @@ 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]; } ); + // m_finiteElement.template computeSurfaceTerms(xLocal, [&] (const int c1, const int c2, const int f1, const int , const int , const int ,const int i2, const int j2, const int k2, real64 val) + // { + // //We take the neighbour element + // const localIndex elemNeigh = elemsToOpposite(k,f1); FE_TYPE::computeSurfaceTerms(xLocal, [&] (const int c1, const int c2, const int f1, const int , const int , const int ,const int i2, const int j2, const int k2, real64 val) { //We take the neighbour element @@ -495,86 +511,80 @@ struct PressureComputationKernel // permutation value contained inside elemsToOppositePermutation permutation const int perm = elemsToOppositePermutation(elemNeigh,f1); - const int p1 = perm%4-1; const int p2 = (perm/4)%4-1; const int p3 = (perm/16)%4-1; - // const int p4 = (perm/64)-1; // Then we transform the 3 indices returned by the callback (i2,j2,k2) using the permutations. One of this permutation, will be 0 (depending on which - // degree of freedom is the one at the opposite of the face shared with the neighbour element) and will correspond to the one where p* will be negative - + // degree of freedom is the one at the opposite of the face shared with the neighbour element) and will correspond to the one where p* will be negati const int Indices[3] = {i2,j2,k2}; - const int ii2 = p1 < 0 ? 0 : Indices[p1]; const int jj2 = p2 < 0 ? 0 : Indices[p2]; const int kk2 = p3 < 0 ? 0 : Indices[p3]; - // Finally, using the dofIndex function, we compute the number of the global degree of freedom on the element const int neighDof = FE_TYPE::dofIndex( ii2, jj2, kk2 ); - //Flux computation + //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]; - }, - [&] (const int c1, const int c2, const int f1, const int i1, const int j1, const int k1, const int i2, const int j2, const int k2, real64 val) - { - //We take the neighbour element - const int elemNeigh = elemsToOpposite(k,f1); + }, + [&] (const int c1, const int c2, const int f1, const int i1, const int j1, const int k1, const int i2, const int j2, const int k2, real64 val) + { + //We take the neighbour element + const int elemNeigh = elemsToOpposite(k,f1); - // Now we seek the degree of freedom on the neighbour element to use for the computation of the flux (or the penalty) - // First, we compute the four possible values of the permutation of the degrees of freedom depending on the the fixed - // permutation value contained inside elemsToOppositePermutation permutation + // Now we seek the degree of freedom on the neighbour element to use for the computation of the flux (or the penalty) + // First, we compute the four possible values of the permutation of the degrees of freedom depending on the the fixed + // permutation value contained inside elemsToOppositePermutation permutation - const int perm = elemsToOppositePermutation(elemNeigh,f1); + const int perm = elemsToOppositePermutation(elemNeigh,f1); - const int p1 = perm%4-1; - const int p2 = (perm/4)%4-1; - const int p3 = (perm/16)%4-1; + const int p1 = perm%4-1; + const int p2 = (perm/4)%4-1; + const int p3 = (perm/16)%4-1; + // Then we transform the 3 indices returned by the callback (i2,j2,k2) using the permutations. One of this permutation, will be 0 (depending on which + // degree of freedom is the one at the opposite of the face shared with the neighbour element) and will correspond to the one where p* will be negative - // Then we transform the 3 indices returned by the callback (i2,j2,k2) using the permutations. One of this permutation, will be 0 (depending on which - // degree of freedom is the one at the opposite of the face shared with the neighbour element) and will correspond to the one where p* will be negative + const int Indices[3] = {i2,j2,k2}; - const int Indices[3] = {i2,j2,k2}; + const int ii2 = p1 < 0 ? 0 : Indices[p1]; + const int jj2 = p2 < 0 ? 0 : Indices[p2]; + const int kk2 = p3 < 0 ? 0 : Indices[p3]; - const int ii2 = p1 < 0 ? 0 : Indices[p1]; - const int jj2 = p2 < 0 ? 0 : Indices[p2]; - const int kk2 = p3 < 0 ? 0 : Indices[p3]; - - // Finally, using the dofIndex function, we compute the number of the global degree of freedom on the element + // Finally, using the dofIndex function, we compute the number of the global degree of freedom on the element const int neighDof = FE_TYPE::dofIndex( ii2, jj2, kk2 ); - //Flux computation + //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: + //Then we need a second time where we take the transpose of the previous values: - const int IndicesTranspose[3] = {i1,j1,k1}; + const int IndicesTranspose[3] = {i1,j1,k1}; - const int ii1 = p1 < 0 ? 0 : IndicesTranspose[p1]; - const int jj1 = p2 < 0 ? 0 : IndicesTranspose[p2]; - const int kk1 = p3 < 0 ? 0 : IndicesTranspose[p3]; + const int ii1 = p1 < 0 ? 0 : IndicesTranspose[p1]; + const int jj1 = p2 < 0 ? 0 : IndicesTranspose[p2]; + const int kk1 = p3 < 0 ? 0 : IndicesTranspose[p3]; - // Finally, using the dofIndex function, we compute the number of the global degree of freedom on the element + // Finally, using the dofIndex function, we compute the number of the global degree of freedom on the element const int neighDof2 = FE_TYPE::dofIndex( ii1, jj1, kk1 ); - //Flux computation + //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]; - } ); + } ); //Source Injection for( localIndex isrc = 0; isrc < sourceConstants.size( 0 ); ++isrc )