From a294c95b70454d439f6bf6c1db124d043f2b9368 Mon Sep 17 00:00:00 2001 From: aurelien Date: Thu, 30 Jan 2025 12:02:26 +0100 Subject: [PATCH 1/3] Fix bug in map computation --- .../isotropic/AcousticWaveEquationDG.cpp | 14 +- .../AcousticWaveEquationDGKernel.hpp | 153 ++++++++++-------- .../dg/acoustic/shared/AcousticFieldsDG.hpp | 2 +- 3 files changed, 97 insertions(+), 72 deletions(-) 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 c9df7a16f37..fa4523afe06 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/dg/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationDG.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/dg/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationDG.cpp @@ -80,7 +80,7 @@ void AcousticWaveEquationDG::initializePreSubGroups() void AcousticWaveEquationDG::registerDataOnMesh( Group & meshBodies ) { - + WaveSolverBase::registerDataOnMesh( meshBodies ); forDiscretizationOnMeshTargets( meshBodies, [&] ( string const &, MeshLevel & mesh, arrayView1d< string const > const & ) @@ -222,7 +222,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, @@ -282,6 +282,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 ); @@ -300,7 +301,6 @@ void AcousticWaveEquationDG::initializePostInitialConditionsPreSubGroups() elemsToOpposite, elemsToOppositePermutation ); - // Precompute reference mass matrix for non-boundary elements m_referenceInvMassMatrix[ regionIndex ][ subRegionIndex ].resizeDimension< 0, 1 >( FE_TYPE::numNodes, FE_TYPE::numNodes ); array2d< real64 > massMatrix; @@ -378,7 +378,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(); @@ -500,9 +500,9 @@ 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() ); 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 aa399ec1934..16c43ce5174 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/dg/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationDGKernel.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/dg/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationDGKernel.hpp @@ -218,13 +218,23 @@ struct PrecomputeNeighborhoodKernel { forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k1 ) { + printf("sizeeTN=%d\n",elemsToNodes.size()); localIndex vertices[ 4 ] = { elemsToNodes( k1, 0 ), elemsToNodes( k1, 1 ), elemsToNodes( k1, 2 ), elemsToNodes( k1, 3 ) }; - + printf("vertices0=%d\n",elemsToNodes( k1, 0 )); + printf("vertices1=%d\n",elemsToNodes( k1, 1 )); + printf("vertices2=%d\n",elemsToNodes( k1, 2 )); + printf("vertices3=%d\n",elemsToNodes( k1, 3 )); for( int i = 0; i < 4; i++ ) { localIndex k1OrderedVertices[ 3 ]; localIndex f = elemsToFaces( k1, i ); localIndex faceVertices[ 3 ] = { facesToNodes( f, 0 ), facesToNodes( f, 1 ), facesToNodes( f, 2 ) }; + printf("index=%d\n",i); + printf("face=%d\n",f); + printf("sizefacetonodes=%d\n",facesToNodes.size()); + printf("faceVertices0=%d\n",facesToNodes( f, 0 )); + printf("faceVertices1=%d\n",facesToNodes( f, 1 )); + printf("faceVertices2=%d\n",facesToNodes( f, 2 )); // find neighboring element, if any localIndex k2 = facesToElems( f, 0 ); if( k2 == k1 ) @@ -233,8 +243,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++ ) { @@ -247,27 +260,37 @@ struct PrecomputeNeighborhoodKernel if( !found ) { o1 = vertex; + indexo1=k; + printf("o1=%d\n",o1); } 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; + printf("là\n"); + printf("k1=%d\n",k1); + printf("o1=%d\n",o1); + printf("k2=%d\n",k2); + elemsToOpposite( k1, indexo1 ) = k2; + printf("apres"); 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++ ) { @@ -280,7 +303,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 @@ -300,7 +325,7 @@ struct PrecomputeNeighborhoodKernel permutation = permutation + c * ( position + 1 ); c = c * 4; } - elemsToOppositePermutation( k1, o1 ) = permutation; + elemsToOppositePermutation( k1, indexo1 ) = permutation; } } } ); @@ -398,96 +423,96 @@ struct PressureComputation } ); - 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); + // 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); - // 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 p4 = (perm/64)-1; + // 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 + // // 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 = m_finiteElement.dofIndex(ii2,jj2,kk2); + // const int neighDof = m_finiteElement.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*p_n[k][c2]; + // flowx[c1] += 0.5*dt2*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 = m_finiteElement.dofIndex(ii2,jj2,kk2); + // const int neighDof = m_finiteElement.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*p_n[elemNeigh][neighDof]; + // flowx[c1] -= 0.5*dt2*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 = m_finiteElement.dofIndex(ii1,jj1,kk1); + // const int neighDof2 = m_finiteElement.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*p_n[elemNeigh][neighDof2]; + // flowx[c2] += 0.5*dt2*p_n[k][c1]; - } ); + // } ); //Source Injection for( localIndex isrc = 0; isrc < sourceConstants.size( 0 ); ++isrc ) diff --git a/src/coreComponents/physicsSolvers/wavePropagation/dg/acoustic/shared/AcousticFieldsDG.hpp b/src/coreComponents/physicsSolvers/wavePropagation/dg/acoustic/shared/AcousticFieldsDG.hpp index 0c5a4775b7e..7aa25134e43 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/dg/acoustic/shared/AcousticFieldsDG.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/dg/acoustic/shared/AcousticFieldsDG.hpp @@ -131,7 +131,7 @@ DECLARE_FIELD( ElementToOppositePermutation, WRITE_AND_READ, "Map from elements to the permutation of the neighboring element, opposite to each vertex." ); -DECLARE_FIELD( CharacteristicSize +DECLARE_FIELD( CharacteristicSize, "charactersticSize", array1d< real32 >, 0, From 252da892cbfa18413cebc5fe6419e795beeaeb00 Mon Sep 17 00:00:00 2001 From: aurelien Date: Thu, 30 Jan 2025 16:59:56 +0100 Subject: [PATCH 2/3] Bugfix --- .../AcousticWaveEquationDGKernel.hpp | 128 +++++++----------- 1 file changed, 51 insertions(+), 77 deletions(-) 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 33af3a757f4..21ec5506254 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/dg/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationDGKernel.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/dg/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationDGKernel.hpp @@ -219,23 +219,12 @@ struct PrecomputeNeighborhoodKernel { forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k1 ) { - printf("sizeeTN=%d\n",elemsToNodes.size()); localIndex vertices[ 4 ] = { elemsToNodes( k1, 0 ), elemsToNodes( k1, 1 ), elemsToNodes( k1, 2 ), elemsToNodes( k1, 3 ) }; - printf("vertices0=%d\n",elemsToNodes( k1, 0 )); - printf("vertices1=%d\n",elemsToNodes( k1, 1 )); - printf("vertices2=%d\n",elemsToNodes( k1, 2 )); - printf("vertices3=%d\n",elemsToNodes( k1, 3 )); for( int i = 0; i < 4; i++ ) { localIndex k1OrderedVertices[ 3 ]; localIndex f = elemsToFaces( k1, i ); localIndex faceVertices[ 3 ] = { facesToNodes( f, 0 ), facesToNodes( f, 1 ), facesToNodes( f, 2 ) }; - printf("index=%d\n",i); - printf("face=%d\n",f); - printf("sizefacetonodes=%d\n",facesToNodes.size()); - printf("faceVertices0=%d\n",facesToNodes( f, 0 )); - printf("faceVertices1=%d\n",facesToNodes( f, 1 )); - printf("faceVertices2=%d\n",facesToNodes( f, 2 )); // find neighboring element, if any localIndex k2 = facesToElems( f, 0 ); if( k2 == k1 ) @@ -262,7 +251,6 @@ struct PrecomputeNeighborhoodKernel { o1 = vertex; indexo1=k; - printf("o1=%d\n",o1); } else { @@ -279,12 +267,7 @@ struct PrecomputeNeighborhoodKernel } else { - printf("là\n"); - printf("k1=%d\n",k1); - printf("o1=%d\n",o1); - printf("k2=%d\n",k2); elemsToOpposite( k1, indexo1 ) = k2; - printf("apres"); localIndex oppositeElemVertices[ 4 ] = { elemsToNodes( k2, 0 ), elemsToNodes( k2, 1 ), elemsToNodes( k2, 2 ), elemsToNodes( k2, 3 ) }; // find opposite vertex in second element int o2 = -1; @@ -470,10 +453,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 >( size, [=] GEOS_HOST_DEVICE ( localIndex const k ) + 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; @@ -519,98 +506,85 @@ struct PressureComputationKernel //We take the neighbour element const localIndex 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 p1 = perm%4-1; - // const int p2 = (perm/4)%4-1; - // const int p3 = (perm/16)%4-1; - // const int p4 = (perm/64)-1; + 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 - // const int Indices[3] = {i2,j2,k2}; + // 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 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 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 = m_finiteElement.dofIndex(ii2,jj2,kk2); const int neighDof = FE_TYPE::dofIndex( ii2, jj2, kk2 ); - // //Flux computation - - // flowx[c1] -= 0.5*dt2*p_n[k][c2]; - // flowx[c1] += 0.5*dt2*p_n[elemNeigh][neighDof]; + //Flux computation + flowx[c1] -= 0.5*dt2*p_n[k][c2]; + flowx[c1] += 0.5*dt2*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); - // // 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 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 perm = elemsToOppositePermutation(elemNeigh,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 - // const int p1 = perm%4-1; - // const int p2 = (perm/4)%4-1; - // const int p3 = (perm/16)%4-1; + const int perm = elemsToOppositePermutation(elemNeigh,f1); - // // 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 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 - // 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 = m_finiteElement.dofIndex(ii2,jj2,kk2); 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*p_n[elemNeigh][neighDof]; + flowx[c1] -= 0.5*dt2*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 = m_finiteElement.dofIndex(ii1,jj1,kk1); 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*p_n[elemNeigh][neighDof2]; + flowx[c2] += 0.5*dt2*p_n[k][c1]; - // } ); + } ); //Source Injection for( localIndex isrc = 0; isrc < sourceConstants.size( 0 ); ++isrc ) From acfcb44eb83ea4c0e7fcc9b621244ad050e954fc Mon Sep 17 00:00:00 2001 From: aurelien Date: Thu, 30 Jan 2025 19:15:02 +0100 Subject: [PATCH 3/3] Bugfix --- .../isotropic/AcousticWaveEquationDG.cpp | 29 ++++++++++--------- .../AcousticWaveEquationDGKernel.hpp | 18 ++++++------ 2 files changed, 25 insertions(+), 22 deletions(-) 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 770fce42c2b..cff2374886f 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 21ec5506254..705e687565a 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]; } );