Skip to content

Commit

Permalink
wip on DG structure initialization
Browse files Browse the repository at this point in the history
  • Loading branch information
@sframba
sframba committed Dec 6, 2024
1 parent 2168ee4 commit 2358458
Show file tree
Hide file tree
Showing 6 changed files with 547 additions and 248 deletions.
4 changes: 2 additions & 2 deletions src/coreComponents/finiteElement/elementFormulations/BB_Tetrahedron.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -169,7 +169,7 @@ class BB_Tetrahedron final : public FiniteElementBase
static void calcN( localIndex const,
real64 (&)[numNodes] )
{
GEOS_ERROR( "Bernstein-Bézier basis is modal, not nodal." );
GEOS_ERROR( "Bernstein-Bézier basis is modal, not nodal. No quadrature points are defined." );
}

/**
Expand Down Expand Up @@ -491,7 +491,7 @@ class BB_Tetrahedron final : public FiniteElementBase
real64 const (&X)[numNodes][3],
real64 ( &gradN )[numNodes][3] )
{
GEOS_ERROR( "Bernstein-Bézier basis is modal, not nodal." );
GEOS_ERROR( "Bernstein-Bézier basis is modal, not nodal. No quadrature points are defined." );
return 0;
}
/**
Expand Down
103 changes: 57 additions & 46 deletions ...csSolvers/wavePropagation/dg/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationDG.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,7 @@
#include "physicsSolvers/wavePropagation/shared/WaveSolverUtils.hpp"
#include "events/EventManager.hpp"
#include "physicsSolvers/wavePropagation/shared/PrecomputeSourcesAndReceiversKernel.hpp"
#include "physicsSolvers/wavePropagation/dg/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationDGKernel.hpp"

namespace geos
{
Expand All @@ -53,7 +54,7 @@ AcousticWaveEquationDG::AcousticWaveEquationDG( const std::string & name,
setSizedFromParent( 0 ).
setDescription( "Element containing the sources" );

registerWrapper( viewKeyStruct::receiverElemString(), &m_rcvElem ).
registerWrapper( viewKeyStruct::receiverElemString(), &m_receiverElem ).
setInputFlag( InputFlags::FALSE ).
setSizedFromParent( 0 ).
setDescription( "Element containing the receivers" );
Expand Down Expand Up @@ -85,9 +86,6 @@ void AcousticWaveEquationDG::registerDataOnMesh( Group & meshBodies )
{
NodeManager & nodeManager = mesh.getNodeManager();

nodeManager.registerField< acousticfieldsdg::AcousticFreeSurfaceNodeIndicator >( this->getName() );

nodeManager.registerField<acousticfieldsdg::AcousticMassMatrix>( getName() );

FaceManager & faceManager = mesh.getFaceManager();
faceManager.registerField< acousticfieldsdg::AcousticFreeSurfaceFaceIndicator >( this->getName() );
Expand All @@ -102,8 +100,9 @@ void AcousticWaveEquationDG::registerDataOnMesh( Group & meshBodies )
subRegion.registerField< acousticfieldsdg::Pressure_nm1 >( this->getName() );
subRegion.registerField< acousticfieldsdg::Pressure_n >( this->getName() );
subRegion.registerField< acousticfieldsdg::Pressure_np1 >( this->getName() );
subRegion.registerField< acousticfieldsdg::StiffnessVector >( this->getName() );
subRegion.registerField< acousticfieldsdg::AcousticMassMatrix >( this->getName() );

subRegion.registerField< acousticfieldsdg::ElementToOpposite >( this->getName() );
subRegion.registerField< acousticfieldsdg::ElementToOppositePermutation >( this->getName() );

finiteElement::FiniteElementBase const & fe = subRegion.getReference< finiteElement::FiniteElementBase >( getDiscretizationName() );

Expand All @@ -117,8 +116,9 @@ void AcousticWaveEquationDG::registerDataOnMesh( Group & meshBodies )
subRegion.getField< acousticfieldsdg::Pressure_nm1 >().resizeDimension< 1 >( numNodesPerElem );
subRegion.getField< acousticfieldsdg::Pressure_n >().resizeDimension< 1 >( numNodesPerElem );
subRegion.getField< acousticfieldsdg::Pressure_np1 >().resizeDimension< 1 >( numNodesPerElem );
subRegion.getField< acousticfieldsdg::StiffnessVector >().resizeDimension< 1 >( numNodesPerElem );
subRegion.getField< acousticfieldsdg::AcousticMassMatrix >().resizeDimension< 1 >( numNodesPerElem );

subRegion.getField< acousticfieldsdg::ElementToOpposite >().resizeDimension< 1 >( 4 );
subRegion.getField< acousticfieldsdg::ElementToOppositePermutation >().resizeDimension< 1 >( 4 );

} );

Expand Down Expand Up @@ -149,7 +149,7 @@ void AcousticWaveEquationDG::precomputeSourceAndReceiverTerm( MeshLevel & baseMe
arrayView2d< localIndex > const receiverNodeIds = m_receiverNodeIds.toView();
arrayView2d< real64 > const receiverConstants = m_receiverConstants.toView();
arrayView1d< localIndex > const receiverIsLocal = m_receiverIsLocal.toView();
arrayView1d< localIndex > const rcvElem = m_rcvElem.toView();
arrayView1d< localIndex > const receiverElem = m_receiverElem.toView();
receiverNodeIds.setValues< EXEC_POLICY >( -1 );
receiverConstants.setValues< EXEC_POLICY >( -1 );
receiverIsLocal.zero();
Expand Down Expand Up @@ -177,36 +177,32 @@ void AcousticWaveEquationDG::precomputeSourceAndReceiverTerm( MeshLevel & baseMe
{
using FE_TYPE = TYPEOFREF( finiteElement );

constexpr localIndex numNodesPerElem = FE_TYPE::numNodes;
localIndex const numFacesPerElem = elementSubRegion.numFacesPerElement();

// AcousticWaveEquationDGKernels::
// PrecomputeSourceAndReceiverKernel::
// launch< EXEC_POLICY, FE_TYPE >
// ( elementSubRegion.size(),
// numNodesPerElem,
// numFacesPerElem,
// X,
// elemGhostRank,
// elemsToNodes,
// elemsToFaces,
// elemCenter,
// faceNormal,
// faceCenter,
// sourceCoordinates,
// sourceIsAccessible,
// sourceElem,
// sourceNodeIds,
// sourceConstants,
// receiverCoordinates,
// receiverIsLocal,
// rcvElem,
// receiverNodeIds,
// receiverConstants,
// sourceValue,
// dt,
// timeSourceFrequency,
// rickerOrder );
AcousticWaveEquationDGKernels::
PrecomputeSourceAndReceiverKernel::
launch< EXEC_POLICY, FE_TYPE >
( elementSubRegion.size(),
facesToNodes,
nodeCoords,
nodeLocalToGlobalMap,
elemLocalToGlobalMap,
nodesToElements,
baseElemsToNodes,
elemGhostRank,
elemsToNodes,
elemsToFaces,
elemCenter,
sourceCoordinates,
sourceIsAccessible,
sourceElem,
sourceNodeIds,
sourceConstants,
receiverCoordinates,
receiverIsLocal,
receiverElem,
receiverNodeIds,
receiverConstants );
} );
} );

Expand Down Expand Up @@ -237,6 +233,8 @@ void AcousticWaveEquationDG::initializePostInitialConditionsPreSubGroups()

/// get the array of indicators: 1 if the face is on the boundary; 0 otherwise
arrayView1d< integer const > const & facesDomainBoundaryIndicator = faceManager.getDomainBoundaryIndicator();
ArrayOfArraysView< localIndex const > const facesToNodes = faceManager.nodeList().toViewConst();
arrayView2d< localIndex const > const & facesToElems = faceManager.elementList();
arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords = nodeManager.getField< fields::referencePosition32 >().toViewConst();

/// get face to node map
Expand All @@ -252,6 +250,11 @@ void AcousticWaveEquationDG::initializePostInitialConditionsPreSubGroups()
elemManager.forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const,
CellElementSubRegion & elementSubRegion )
{
GEOS_THROW_IF( elementSubRegion.getElementType() != ElementType::Tetrahedron,
"Invalid type of element, the acoustic DG solver is designed for tetrahedral meshes only ",
InputError );

forAll< EXEC_POLICY >
finiteElement::FiniteElementBase const &
fe = elementSubRegion.getReference< finiteElement::FiniteElementBase >( getDiscretizationName() );

Expand All @@ -260,6 +263,19 @@ void AcousticWaveEquationDG::initializePostInitialConditionsPreSubGroups()

computeTargetNodeSet( elemsToNodes, elementSubRegion.size(), fe.getNumQuadraturePoints() );

// Compute auxiliary element to neighbor (opposite to vertex) maps

AcousticWaveEquationDGKernels::
PrecomputeNeighborhoodKernel::
launch< EXEC_POLICY, FE_TYPE >
( elemsToNodes,
elemsToFaces,
facesToElems,
facesToNodes,
freeSurfaceFaceIndicator,
elemsToOpposite,
elemsToOppositePermutation );

// arrayView1d< real32 const > const velocity = elementSubRegion.getField< acousticfieldsdgdgdg::AcousticVelocity >();
// arrayView1d< real32 const > const density = elementSubRegion.getField< acousticfieldsdgdgdg::AcousticDensity >();

Expand Down Expand Up @@ -328,12 +344,8 @@ void AcousticWaveEquationDG::applyFreeSurfaceBC( real64 const time, DomainPartit
/// array of indicators: 1 if a face is on on free surface; 0 otherwise
arrayView1d< localIndex > const freeSurfaceFaceIndicator = faceManager.getField< acousticfieldsdg::AcousticFreeSurfaceFaceIndicator >();

/// array of indicators: 1 if a node is on on free surface; 0 otherwise
arrayView1d< localIndex > const freeSurfaceNodeIndicator = nodeManager.getField< acousticfieldsdg::AcousticFreeSurfaceNodeIndicator >();


freeSurfaceFaceIndicator.zero();
freeSurfaceNodeIndicator.zero();

fsManager.apply< FaceManager >( time,
domain.getMeshBody( 0 ).getMeshLevel( m_discretizationName ),
Expand Down Expand Up @@ -405,7 +417,7 @@ void AcousticWaveEquationDG::prepareNextTimestep( MeshLevel & mesh )
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 stiffnessVector = nodeManager.getField< acousticfieldsdg::StiffnessVector >();
//arrayView2d< real32 > const rhs = nodeManager.getField< acousticfieldsdgdgdg::ForcingRHS >();

SortedArrayView< localIndex const > const solverTargetNodesSet = m_solverTargetNodesSet.toViewConst();
Expand Down Expand Up @@ -437,8 +449,8 @@ void AcousticWaveEquationDG::computeUnknowns( real64 const & time_n,
arrayView2d< real32 > const p_n = nodeManager.getField< acousticfieldsdg::Pressure_n >();
arrayView2d< real32 > const p_np1 = nodeManager.getField< acousticfieldsdg::Pressure_np1 >();

arrayView1d< localIndex const > const freeSurfaceNodeIndicator = nodeManager.getField< acousticfieldsdg::AcousticFreeSurfaceNodeIndicator >();
arrayView2d< real32 > const stiffnessVector = nodeManager.getField< acousticfieldsdg::StiffnessVector >();
// arrayView1d< localIndex const > const freeSurfaceNodeIndicator = nodeManager.getField< acousticfieldsdg::AcousticFreeSurfaceNodeIndicator >();
// arrayView2d< real32 > const stiffnessVector = nodeManager.getField< acousticfieldsdg::StiffnessVector >();
// arrayView1d< real32 > const rhs = nodeManager.getField< acousticfieldsdgdg::ForcingRHS >();

// auto kernelFactory = acousticWaveEquationSEMKernels::ExplicitAcousticSEMFactory( dt );
Expand Down Expand Up @@ -474,7 +486,7 @@ void AcousticWaveEquationDG::synchronizeUnknowns( real64 const & time_n,
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 stiffnessVector = nodeManager.getField< acousticfieldsdg::StiffnessVector >();
//arrayView1d< real32 > const rhs = nodeManager.getField< acousticfieldsdg::ForcingRHS >();

/// synchronize pressure fields
Expand Down Expand Up @@ -505,7 +517,6 @@ real64 AcousticWaveEquationDG::explicitStepInternal( real64 const & time_n,

GEOS_LOG_RANK_0_IF( dt < epsilonLoc, "Warning! Value for dt: " << dt << "s is smaller than local threshold: " << epsilonLoc );

real64 dtCompute;
forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
MeshLevel & mesh,
arrayView1d< string const > const & regionNames )
Expand Down
13 changes: 11 additions & 2 deletions ...csSolvers/wavePropagation/dg/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationDG.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -97,7 +97,7 @@ class AcousticWaveEquationDG : public WaveSolverBase
static constexpr char const * pressureNp1AtReceiversString() { return "pressureNp1AtReceivers"; }

static constexpr char const * sourceElemString() { return "sourceElem"; }
static constexpr char const * receiverElemString() { return "rcvElem"; }
static constexpr char const * receiverElemString() { return "receiverElem"; }
static constexpr char const * receiverRegionString() { return "receiverRegion"; }

} waveEquationViewKeys;
Expand Down Expand Up @@ -167,11 +167,20 @@ class AcousticWaveEquationDG : public WaveSolverBase
array1d< localIndex > m_sourceElem;

/// Array containing the elements which contain a receiver
array1d< localIndex > m_rcvElem;
array1d< localIndex > m_receiverElem;

/// Array containing the elements which contain the region which the receiver belongs
array1d< localIndex > m_receiverRegion;

/// Inverse of the mass matrix in the reference element for each subregion
array1d< array2d< real64 > > m_referenceInvMassMatrix

/// Inverse of the mass plus damping matrix in the reference element for each boundary element
array3d< real64 > m_boundaryInvMassPlusDamping

/// Index for each boundary element to m_boundaryInvMassPlusDamping
array3d< real64 > m_indexToBoundaryMatrix

};

} /* namespace geos */
Expand Down
Loading

0 comments on commit 2358458

Please sign in to comment.