diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index 0ebe2652f98..97f13dd9062 100644 --- a/.integrated_tests.yaml +++ b/.integrated_tests.yaml @@ -1,6 +1,6 @@ baselines: bucket: geosx - baseline: integratedTests/baseline_integratedTests-pr3381-9063-399123c + baseline: integratedTests/baseline_integratedTests-pr3361-9139-2fc4131 allow_fail: all: '' streak: '' diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md index 1fbd4cbfa9d..e5e337912cb 100644 --- a/BASELINE_NOTES.md +++ b/BASELINE_NOTES.md @@ -6,11 +6,18 @@ This file is designed to track changes to the integrated test baselines. Any developer who updates the baseline ID in the .integrated_tests.yaml file is expected to create an entry in this file with the pull request number, date, and their justification for rebaselining. These notes should be in reverse-chronological order, and use the following time format: (YYYY-MM-DD). +PR #3361 (2024-12-03) +===================== +Baseline diffs after reimplementation of wave equation acoustic gradient for velocity and density parameters: new field "partialGradient2" and "pressureForward" field replacing "pressureDoubleDerivative". + +PR #3393 (2024-12-02) +===================== +Fix netToGross bug. + PR #3381 (2024-12-01) ===================== A few baseline diffs for order FaceElementSubRegion::m_toFacesRelation map. Not sure why this was changed by this PR, but the previous order seems incorrect for a couple of cases. - PR #2957 (2024-11-27) ===================== Added ExternalDataRepository. @@ -19,7 +26,6 @@ PR #3448 (2024-11-21) ===================== Switched the FaceElementSubRegion::m_toFacesRelation and FaceElementSubRegion::m_2dElemToElems back to array2d instead of ArrayOfArray. This results in a reordering m_toFacesRelation back to the "correct" assumed order of "original face first". This fixes a bug that failed to remove the CellStencil entry when a FaceElement splits two cells. - PR #2637 (2024-11-21) ===================== Added numberOfTargetProcesses. diff --git a/inputFiles/lagrangianContactMechanics/LagrangeContactBubbleStab_singleFracCompression_base.xml b/inputFiles/lagrangianContactMechanics/LagrangeContactBubbleStab_singleFracCompression_base.xml index 09ec110a624..854136ed032 100644 --- a/inputFiles/lagrangianContactMechanics/LagrangeContactBubbleStab_singleFracCompression_base.xml +++ b/inputFiles/lagrangianContactMechanics/LagrangeContactBubbleStab_singleFracCompression_base.xml @@ -130,8 +130,8 @@ plotLevel="2" format="binary" plotFileRoot="bubbleStab"/> - - + + diff --git a/inputFiles/singlePhaseFlow/FieldCaseTutorial3_composite_smoke.xml b/inputFiles/singlePhaseFlow/FieldCaseTutorial3_composite_smoke.xml index dd43e24c23a..9312ecaae37 100644 --- a/inputFiles/singlePhaseFlow/FieldCaseTutorial3_composite_smoke.xml +++ b/inputFiles/singlePhaseFlow/FieldCaseTutorial3_composite_smoke.xml @@ -15,7 +15,7 @@ + dataSourceName="synthetic"> string toMetricPrefixString( T const & value ) { + if( std::fpclassify( value ) == FP_ZERO ) + { + return " 0.0 "; + } + // These are the metric prefixes corrosponding to kilo, mega, giga...etc. char const prefixes[12] = { 'f', 'p', 'n', 'u', 'm', ' ', 'K', 'M', 'G', 'T', 'P', 'E'}; string rval; diff --git a/src/coreComponents/linearAlgebra/docs/DofManager.rst b/src/coreComponents/linearAlgebra/docs/DofManager.rst index d13e2ee7c9c..6948b36f1ae 100644 --- a/src/coreComponents/linearAlgebra/docs/DofManager.rst +++ b/src/coreComponents/linearAlgebra/docs/DofManager.rst @@ -2,8 +2,6 @@ DoF Manager ############################################################################### -This will contains a description of the DoF manager in GEOS. - Brief description ======================== diff --git a/src/coreComponents/linearAlgebra/docs/LinearSolvers.rst b/src/coreComponents/linearAlgebra/docs/LinearSolvers.rst index b4b68dfc682..6f68efac64d 100644 --- a/src/coreComponents/linearAlgebra/docs/LinearSolvers.rst +++ b/src/coreComponents/linearAlgebra/docs/LinearSolvers.rst @@ -14,14 +14,14 @@ Any physics solver relying on standard finite element and finite volume techniqu \mathsf{A} \mathsf{x} = \mathsf{b} -with a :math:`\mathsf{A}` a square sparse matrix, :math:`\mathsf{x}` the solution vector, and :math:`\mathsf{b}` the right-hand side. +where :math:`\mathsf{A}` is the square sparse matrix, :math:`\mathsf{x}` the solution vector, and :math:`\mathsf{b}` the right-hand side. For example, in a classical linear elastostatics problem :math:`\mathsf{A}` is the stiffness matrix, and :math:`\mathsf{x}` and :math:`\mathsf{b}` are the displacement and nodal force vectors, respectively. This solution stage represents the most computationally expensive portion of a typical simulation. Solution algorithms generally belong to two families of methods: direct methods and iterative methods. In GEOS both options are made available wrapping around well-established open-source linear algebra libraries, namely -`HYPRE `__, +`HYPRE `__, `PETSc `__, `SuperLU `__, and `Trilinos `__. @@ -38,7 +38,7 @@ Irrespective of the selected direct solver implementation, three stages can be i (#) **Solve Stage**: the solution to the linear systems involving the factorized matrix is computed (#) **Finalize Stage**: the systems involving the factorized matrix have been solved and the direct solver lifetime ends -The default option in GEOS relies on `SuperLU `__, a general purpose library for the direct solution of large, sparse, nonsymmetric systems of linear equations, that is called taking advantage of the interface provided in `HYPRE `__. +The default option in GEOS relies on `SuperLU `__, a general purpose library for the direct solution of large, sparse, nonsymmetric systems of linear equations, that is called taking advantage of the interface provided in `HYPRE `__. ****************** Iterative methods @@ -110,8 +110,8 @@ This section provides a brief description of the available preconditioners. - `PETSc documentation `__, - `Trilinos documentation `__. -* **MGR**: multigrid reduction. Available through *hypre* interface only. Specific documentation coming soon. - Further details can be found in `MGR documentation `__. +* **MGR**: multigrid reduction. Available through *hypre* interface only. + Further details can be found in `MGR documentation `__, also see section below. * **Block**: custom preconditioner designed for a 2 x 2 block matrix. @@ -119,7 +119,7 @@ This section provides a brief description of the available preconditioners. HYPRE MGR Preconditioner ************************ -MGR stands for multigrid reduction, a multigrid method that uses the interpolation, restriction operators, and the Galerkin triple product, to reduce a linear system to a smaller one, similar to a Schur complement approach. As such, it is designed to target block linear systems resulting from discretizations of multiphysics problems. GEOS uses MGR through an implementation in `HYPRE `__. More information regarding MGR can be found `here `__. Currently, MGR strategies are implemented for hydraulic fracturing, poroelastic, compositional flow with and without wells. More multiphysics solvers with MGR will be enabled in the future. +MGR stands for multigrid reduction, a multigrid method that uses the interpolation, restriction operators, and the Galerkin triple product, to reduce a linear system to a smaller one, similar to a Schur complement approach. As such, it is designed to target block linear systems resulting from discretizations of multiphysics problems. GEOS uses MGR through an implementation in `HYPRE `__. More information regarding MGR can be found `here `__. Currently, MGR strategies are implemented for hydraulic fracturing, singlephase and multiphase poromechanics, singlephase poromechanics with fractures, compositional flow with and without wells. More multiphysics solvers with MGR will be enabled in the future. To use MGR for a specific block system, several components need to be specified. @@ -157,3 +157,16 @@ Moreover, a block scaling is available. Feasible options are: * none: keep the original scaling; * Frobenius norm: equilibrate Frobenius norm of the diagonal blocks; * user provided. + +******************** +Adaptive tolerance +******************** + +This feature is available for iterative solvers and can be enabled using `krylovAdaptiveTol` flag in `LinearSolverParameters`. It follows the Eisenstat-Walker inexact Newton approach described in [Eisenstat and Walker 1996]. The key idea is to relax the linear solver tolerance at the beginning of the nonlinear iterations loop and tighten it when getting closer to the final solution. The initial tolerance is defined by `krylovWeakestTol` and starting from second nonlinear iteration the tolerance is chosen using the following steps: + +- compute the current to previous nonlinear norm ratio: :math:`\mathsf{nr} = \mathsf{min}( \mathsf{norm}^{curr} / \mathsf{norm}^{prev}, 1.0 )` +- estimate the new linear solver tolerance: :math:`\mathsf{tol}_{new} = \mathsf{\gamma} \cdot \mathsf{nr}^{ax}` +- compute a safeguard to avoid too sharp tolerance reduction: :math:`\mathsf{tol}_{alt} = \mathsf{tol}_{old}^{2}` (the bound is the quadratic reduction with respect to the previous tolerance value) +- apply safeguards and compute the final tolerance: :math:`\mathsf{tol} = \mathsf{max}( \mathsf{tol}_{new}, \mathsf{tol}_{alt} )`, :math:`\mathsf{tol} = \mathsf{min}( \mathsf{tol}_{max}, \mathsf{max}( \mathsf{tol}_{min}, \mathsf{tol} ) )` + +Here :math:`\mathsf{\gamma}` is the forcing term, :math:`ax` is the adaptivity exponent, :math:`\mathsf{tol}_{min}` and :math:`\mathsf{tol}_{max}` are prescribed tolerance bounds (defined by `krylovStrongestTol` and `krylovWeakestTol`, respectively). diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp index 08e7176fdbe..9c7d0108c7b 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp @@ -831,30 +831,25 @@ real64 CompositionalMultiphaseBase::updateFluidState( ElementSubRegionBase & sub } void CompositionalMultiphaseBase::initializeFluidState( MeshLevel & mesh, - DomainPartition & domain, arrayView1d< string const > const & regionNames ) { GEOS_MARK_FUNCTION; integer const numComp = m_numComponents; - // 1. Compute hydrostatic equilibrium in the regions for which corresponding field specification tag has been specified - computeHydrostaticEquilibrium(); - mesh.getElemManager().forElementSubRegions( regionNames, [&]( localIndex const, ElementSubRegionBase & subRegion ) { - // 2. Assume global component fractions have been prescribed. + // Assume global component fractions have been prescribed. // Initialize constitutive state to get fluid density. updateFluidModel( subRegion ); - // 3. Back-calculate global component densities from fractions and total fluid density + // Back-calculate global component densities from fractions and total fluid density // in order to initialize the primary solution variables - string const & fluidName = subRegion.getReference< string >( viewKeyStruct::fluidNamesString() ); + string const & fluidName = subRegion.template getReference< string >( viewKeyStruct::fluidNamesString() ); MultiFluidBase const & fluid = getConstitutiveModel< MultiFluidBase >( subRegion, fluidName ); arrayView2d< real64 const, multifluid::USD_FLUID > const totalDens = fluid.totalDensity(); - arrayView2d< real64 const, compflow::USD_COMP > const compFrac = subRegion.getField< fields::flow::globalCompFraction >(); arrayView2d< real64, compflow::USD_COMP > const compDens = @@ -867,10 +862,13 @@ void CompositionalMultiphaseBase::initializeFluidState( MeshLevel & mesh, compDens[ei][ic] = totalDens[ei][0] * compFrac[ei][ic]; } } ); - } ); - // with initial component densities defined - check if they need to be corrected to avoid zero diags etc - chopNegativeDensities( domain ); + // with initial component densities defined - check if they need to be corrected to avoid zero diags etc + if( m_allowCompDensChopping ) + { + chopNegativeDensities( subRegion ); + } + } ); // for some reason CUDA does not want the host_device lambda to be defined inside the generic lambda // I need the exact type of the subRegion for updateSolidflowProperties to work well. @@ -878,115 +876,65 @@ void CompositionalMultiphaseBase::initializeFluidState( MeshLevel & mesh, SurfaceElementSubRegion >( regionNames, [&]( localIndex const, auto & subRegion ) { - // 4. Initialize/update dependent state quantities + // Initialize/update dependent state quantities - // 4.1 Update the constitutive models that only depend on - // - the primary variables - // - the fluid constitutive quantities (as they have already been updated) - // We postpone the other constitutive models for now - // In addition, to avoid multiplying permeability/porosity bay netToGross in the assembly kernel, we do it once and for all here - arrayView1d< real64 const > const netToGross = subRegion.template getField< fields::flow::netToGross >(); - CoupledSolidBase const & porousSolid = - getConstitutiveModel< CoupledSolidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::solidNamesString() ) ); - PermeabilityBase const & permeabilityModel = - getConstitutiveModel< PermeabilityBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::permeabilityNamesString() ) ); - permeabilityModel.scaleHorizontalPermeability( netToGross ); - porousSolid.scaleReferencePorosity( netToGross ); - saveConvergedState( subRegion ); // necessary for a meaningful porosity update in sequential schemes - updatePorosityAndPermeability( subRegion ); updateCompAmount( subRegion ); updatePhaseVolumeFraction( subRegion ); - // Now, we initialize and update each constitutive model one by one + // Update the constitutive models that only depend on + // - the primary variables + // - the fluid constitutive quantities (as they have already been updated) + // We postpone the other constitutive models for now - // 4.2 Save the computed porosity into the old porosity - // - // Note: - // - This must be called after updatePorosityAndPermeability - // - This step depends on porosity - string const & solidName = subRegion.template getReference< string >( viewKeyStruct::solidNamesString() ); - CoupledSolidBase const & porousMaterial = getConstitutiveModel< CoupledSolidBase >( subRegion, solidName ); - porousMaterial.initializeState(); - - // 4.3 Initialize/update the relative permeability model using the initial phase volume fraction - // This is needed to handle relative permeability hysteresis - // Also, initialize the fluid model - // - // Note: - // - This must be called after updatePhaseVolumeFraction - // - This step depends on phaseVolFraction + // Now, we initialize and update each constitutive model one by one // initialized phase volume fraction arrayView2d< real64 const, compflow::USD_PHASE > const phaseVolFrac = subRegion.template getField< fields::flow::phaseVolumeFraction >(); - string const & relpermName = subRegion.template getReference< string >( viewKeyStruct::relPermNamesString() ); - RelativePermeabilityBase & relPermMaterial = - getConstitutiveModel< RelativePermeabilityBase >( subRegion, relpermName ); - relPermMaterial.saveConvergedPhaseVolFractionState( phaseVolFrac ); // this needs to happen before calling updateRelPermModel + // Initialize/update the relative permeability model using the initial phase volume fraction + // Note: + // - This must be called after updatePhaseVolumeFraction + // - This step depends on phaseVolFraction + RelativePermeabilityBase & relPerm = + getConstitutiveModel< RelativePermeabilityBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::relPermNamesString() ) ); + relPerm.saveConvergedPhaseVolFractionState( phaseVolFrac ); // this needs to happen before calling updateRelPermModel updateRelPermModel( subRegion ); - relPermMaterial.saveConvergedState(); // this needs to happen after calling updateRelPermModel + relPerm.saveConvergedState(); // this needs to happen after calling updateRelPermModel string const & fluidName = subRegion.template getReference< string >( viewKeyStruct::fluidNamesString() ); - MultiFluidBase & fluidMaterial = getConstitutiveModel< MultiFluidBase >( subRegion, fluidName ); - fluidMaterial.initializeState(); + MultiFluidBase & fluid = getConstitutiveModel< MultiFluidBase >( subRegion, fluidName ); + fluid.initializeState(); + + // Update the phase mobility + // Note: + // - This must be called after updateRelPermModel + // - This step depends phaseRelPerm + updatePhaseMobility( subRegion ); - // 4.4 Then, we initialize/update the capillary pressure model - // + // Initialize/update the capillary pressure model // Note: - // - This must be called after updatePorosityAndPermeability - // - This step depends on porosity and permeability + // - This must be called after updatePorosityAndPermeability and updatePhaseVolumeFraction + // - This step depends on porosity, permeability, and phaseVolFraction if( m_hasCapPressure ) { // initialized porosity - arrayView2d< real64 const > const porosity = porousMaterial.getPorosity(); + CoupledSolidBase const & porousSolid = + getConstitutiveModel< CoupledSolidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::solidNamesString() ) ); + arrayView2d< real64 const > const porosity = porousSolid.getPorosity(); - string const & permName = subRegion.template getReference< string >( viewKeyStruct::permeabilityNamesString() ); - PermeabilityBase const & permeabilityMaterial = - getConstitutiveModel< PermeabilityBase >( subRegion, permName ); // initialized permeability + PermeabilityBase const & permeabilityMaterial = + getConstitutiveModel< PermeabilityBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::permeabilityNamesString() ) ); arrayView3d< real64 const > const permeability = permeabilityMaterial.permeability(); - string const & capPressureName = subRegion.template getReference< string >( viewKeyStruct::capPressureNamesString() ); - CapillaryPressureBase const & capPressureMaterial = - getConstitutiveModel< CapillaryPressureBase >( subRegion, capPressureName ); - capPressureMaterial.initializeRockState( porosity, permeability ); // this needs to happen before calling updateCapPressureModel + CapillaryPressureBase const & capPressure = + getConstitutiveModel< CapillaryPressureBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::capPressureNamesString() ) ); + capPressure.initializeRockState( porosity, permeability ); // this needs to happen before calling updateCapPressureModel updateCapPressureModel( subRegion ); } - // 4.5 Update the phase mobility - // - // Note: - // - This must be called after updateRelPermModel - // - This step depends phaseRelPerm - updatePhaseMobility( subRegion ); - - // 4.6 We initialize the rock thermal quantities: conductivity and solid internal energy - // - // Note: - // - This must be called after updatePorosityAndPermeability and updatePhaseVolumeFraction - // - This step depends on porosity and phaseVolFraction - if( m_isThermal ) - { - // initialized porosity - arrayView2d< real64 const > const porosity = porousMaterial.getPorosity(); - - string const & thermalConductivityName = subRegion.template getReference< string >( viewKeyStruct::thermalConductivityNamesString() ); - MultiPhaseThermalConductivityBase const & conductivityMaterial = - getConstitutiveModel< MultiPhaseThermalConductivityBase >( subRegion, thermalConductivityName ); - conductivityMaterial.initializeRockFluidState( porosity, phaseVolFrac ); - // note that there is nothing to update here because thermal conductivity is explicit for now - - updateSolidInternalEnergyModel( subRegion ); - string const & solidInternalEnergyName = subRegion.template getReference< string >( viewKeyStruct::solidInternalEnergyNamesString() ); - SolidInternalEnergy const & solidInternalEnergyMaterial = - getConstitutiveModel< SolidInternalEnergy >( subRegion, solidInternalEnergyName ); - solidInternalEnergyMaterial.saveConvergedState(); - - updateEnergy( subRegion ); - } - - // Step 4.7: if the diffusion and/or dispersion is/are supported, initialize the two models + // If the diffusion and/or dispersion is/are supported, initialize the two models if( m_hasDiffusion ) { string const & diffusionName = subRegion.template getReference< string >( viewKeyStruct::diffusionNamesString() ); @@ -1004,24 +952,42 @@ void CompositionalMultiphaseBase::initializeFluidState( MeshLevel & mesh, } } ); +} - // 5. Save initial pressure - mesh.getElemManager().forElementSubRegions( regionNames, [&]( localIndex const, - ElementSubRegionBase & subRegion ) +void CompositionalMultiphaseBase::initializeThermalState( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) +{ + mesh.getElemManager().forElementSubRegions< CellElementSubRegion, + SurfaceElementSubRegion >( regionNames, [&]( localIndex const, + auto & subRegion ) { - arrayView1d< real64 const > const pres = subRegion.getField< fields::flow::pressure >(); - arrayView1d< real64 > const initPres = subRegion.getField< fields::flow::initialPressure >(); - arrayView1d< real64 const > const temp = subRegion.getField< fields::flow::temperature >(); - arrayView1d< real64 > const initTemp = subRegion.template getField< fields::flow::initialTemperature >(); - initPres.setValues< parallelDevicePolicy<> >( pres ); - initTemp.setValues< parallelDevicePolicy<> >( temp ); + // initialized porosity + CoupledSolidBase const & porousSolid = + getConstitutiveModel< CoupledSolidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::solidNamesString() ) ); + arrayView2d< real64 const > const porosity = porousSolid.getPorosity(); + + // initialized phase volume fraction + arrayView2d< real64 const, compflow::USD_PHASE > const phaseVolFrac = + subRegion.template getField< fields::flow::phaseVolumeFraction >(); + + string const & thermalConductivityName = subRegion.template getReference< string >( viewKeyStruct::thermalConductivityNamesString()); + MultiPhaseThermalConductivityBase const & conductivityMaterial = + getConstitutiveModel< MultiPhaseThermalConductivityBase >( subRegion, thermalConductivityName ); + conductivityMaterial.initializeRockFluidState( porosity, phaseVolFrac ); + // note that there is nothing to update here because thermal conductivity is explicit for now + + updateSolidInternalEnergyModel( subRegion ); + string const & solidInternalEnergyName = subRegion.template getReference< string >( viewKeyStruct::solidInternalEnergyNamesString()); + SolidInternalEnergy const & solidInternalEnergyMaterial = + getConstitutiveModel< SolidInternalEnergy >( subRegion, solidInternalEnergyName ); + solidInternalEnergyMaterial.saveConvergedState(); + + updateEnergy( subRegion ); } ); } -void CompositionalMultiphaseBase::computeHydrostaticEquilibrium() +void CompositionalMultiphaseBase::computeHydrostaticEquilibrium( DomainPartition & domain ) { FieldSpecificationManager & fsManager = FieldSpecificationManager::getInstance(); - DomainPartition & domain = this->getGroupByPath< DomainPartition >( "/Problem/domain" ); integer const numComps = m_numComponents; integer const numPhases = m_numPhases; @@ -1281,8 +1247,7 @@ void CompositionalMultiphaseBase::initializePostInitialConditionsPreSubGroups() arrayView1d< string const > const & regionNames ) { FieldIdentifiers fieldsToBeSync; - fieldsToBeSync.addElementFields( { fields::flow::pressure::key(), - fields::flow::globalCompDensity::key() }, + fieldsToBeSync.addElementFields( { fields::flow::globalCompDensity::key() }, regionNames ); CommunicationTools::getInstance().synchronizeFields( fieldsToBeSync, mesh, domain.getNeighbors(), false ); @@ -1295,35 +1260,10 @@ void CompositionalMultiphaseBase::initializePostInitialConditionsPreSubGroups() string const & fluidName = subRegion.template getReference< string >( viewKeyStruct::fluidNamesString() ); MultiFluidBase & fluid = getConstitutiveModel< MultiFluidBase >( subRegion, fluidName ); fluid.setMassFlag( m_useMass ); - - saveConvergedState( subRegion ); // necessary for a meaningful porosity update in sequential schemes - updatePorosityAndPermeability( subRegion ); - - CoupledSolidBase const & porousSolid = - getConstitutiveModel< CoupledSolidBase >( subRegion, - subRegion.template getReference< string >( viewKeyStruct::solidNamesString() ) ); - porousSolid.initializeState(); } ); - - // Initialize primary variables from applied initial conditions - initializeFluidState( mesh, domain, regionNames ); - - mesh.getElemManager().forElementRegions< SurfaceElementRegion >( regionNames, - [&]( localIndex const, - SurfaceElementRegion & region ) - { - region.forElementSubRegions< FaceElementSubRegion >( [&]( FaceElementSubRegion & subRegion ) - { - subRegion.getWrapper< real64_array >( fields::flow::hydraulicAperture::key() ). - setApplyDefaultValue( region.getDefaultAperture() ); - } ); - } ); - } ); - // report to the user if some pore volumes are very small - // note: this function is here because: 1) porosity has been initialized and 2) NTG has been applied - validatePoreVolumes( domain ); + initializeState( domain ); } void @@ -2071,9 +2011,6 @@ void CompositionalMultiphaseBase::chopNegativeDensities( DomainPartition & domai using namespace isothermalCompositionalMultiphaseBaseKernels; - integer const numComp = m_numComponents; - real64 const minCompDens = m_minCompDens; - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, MeshLevel & mesh, arrayView1d< string const > const & regionNames ) @@ -2082,25 +2019,33 @@ void CompositionalMultiphaseBase::chopNegativeDensities( DomainPartition & domai [&]( localIndex const, ElementSubRegionBase & subRegion ) { - arrayView1d< integer const > const ghostRank = subRegion.ghostRank(); + chopNegativeDensities( subRegion ); + } ); + } ); +} - arrayView2d< real64, compflow::USD_COMP > const compDens = - subRegion.getField< fields::flow::globalCompDensity >(); +void CompositionalMultiphaseBase::chopNegativeDensities( ElementSubRegionBase & subRegion ) +{ + integer const numComp = m_numComponents; + real64 const minCompDens = m_minCompDens; - forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei ) + arrayView1d< integer const > const ghostRank = subRegion.ghostRank(); + + arrayView2d< real64, compflow::USD_COMP > const compDens = + subRegion.getField< fields::flow::globalCompDensity >(); + + forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei ) + { + if( ghostRank[ei] < 0 ) + { + for( integer ic = 0; ic < numComp; ++ic ) { - if( ghostRank[ei] < 0 ) + if( compDens[ei][ic] < minCompDens ) { - for( integer ic = 0; ic < numComp; ++ic ) - { - if( compDens[ei][ic] < minCompDens ) - { - compDens[ei][ic] = minCompDens; - } - } + compDens[ei][ic] = minCompDens; } - } ); - } ); + } + } } ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp index cfd06428707..278d8d6b453 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp @@ -282,12 +282,14 @@ class CompositionalMultiphaseBase : public FlowSolverBase * from prescribed intermediate values (i.e. global densities from global fractions) * and any applicable hydrostatic equilibration of the domain */ - void initializeFluidState( MeshLevel & mesh, DomainPartition & domain, arrayView1d< string const > const & regionNames ); + virtual void initializeFluidState( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; + + virtual void initializeThermalState( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; /** * @brief Compute the hydrostatic equilibrium using the compositions and temperature input tables */ - void computeHydrostaticEquilibrium(); + virtual void computeHydrostaticEquilibrium( DomainPartition & domain ) override; /** * @brief Function to perform the Application of Dirichlet type BC's @@ -362,6 +364,8 @@ class CompositionalMultiphaseBase : public FlowSolverBase */ void chopNegativeDensities( DomainPartition & domain ); + void chopNegativeDensities( ElementSubRegionBase & subRegion ); + virtual real64 setNextDtBasedOnStateChange( real64 const & currentDt, DomainPartition & domain ) override; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp index 4c9dde33e1b..2e99cfe6f48 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp @@ -36,8 +36,6 @@ #include "mesh/mpiCommunications/CommunicationTools.hpp" #include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" #include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp" -#include "physicsSolvers/fluidFlow/kernels/compositional/AccumulationKernel.hpp" -#include "physicsSolvers/fluidFlow/kernels/compositional/ThermalAccumulationKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/ResidualNormKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/ThermalResidualNormKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/SolutionScalingKernel.hpp" @@ -50,10 +48,10 @@ #include "physicsSolvers/fluidFlow/kernels/compositional/ThermalDiffusionDispersionFluxComputeKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/StabilizedFluxComputeKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/DissipationFluxComputeKernel.hpp" -#include "physicsSolvers/fluidFlow/kernels/compositional/PhaseMobilityKernel.hpp" -#include "physicsSolvers/fluidFlow/kernels/compositional/ThermalPhaseMobilityKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/DirichletFluxComputeKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/ThermalDirichletFluxComputeKernel.hpp" +#include "physicsSolvers/fluidFlow/kernels/compositional/PhaseMobilityKernel.hpp" +#include "physicsSolvers/fluidFlow/kernels/compositional/ThermalPhaseMobilityKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/AquiferBCKernel.hpp" namespace geos diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp index bedfa0c3261..d6a0c9568d6 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp @@ -31,7 +31,6 @@ #include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" #include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/CompositionalMultiphaseHybridFVMKernels.hpp" -#include "physicsSolvers/fluidFlow/kernels/compositional/AccumulationKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/SolutionScalingKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/SolutionCheckKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/ResidualNormKernel.hpp" diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp index f803e7c7e1e..8cd3ff1b942 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp @@ -286,16 +286,6 @@ void FlowSolverBase::saveSequentialIterationState( DomainPartition & domain ) m_sequentialTempChange = m_isThermal ? MpiWrapper::max( maxTempChange ) : 0.0; } -void FlowSolverBase::enableFixedStressPoromechanicsUpdate() -{ - m_isFixedStressPoromechanicsUpdate = true; -} - -void FlowSolverBase::enableJumpStabilization() -{ - m_isJumpStabilized = true; -} - void FlowSolverBase::setConstitutiveNamesCallSuper( ElementSubRegionBase & subRegion ) const { PhysicsSolverBase::setConstitutiveNamesCallSuper( subRegion ); @@ -455,6 +445,12 @@ void FlowSolverBase::initializePostInitialConditionsPreSubGroups() arrayView1d< string const > const & regionNames ) { precomputeData( mesh, regionNames ); + + FieldIdentifiers fieldsToBeSync; + fieldsToBeSync.addElementFields( { fields::flow::pressure::key(), fields::flow::temperature::key() }, + regionNames ); + + CommunicationTools::getInstance().synchronizeFields( fieldsToBeSync, mesh, domain.getNeighbors(), false ); } ); } @@ -491,6 +487,97 @@ void FlowSolverBase::precomputeData( MeshLevel & mesh, } } +void FlowSolverBase::initializeState( DomainPartition & domain ) +{ + // Compute hydrostatic equilibrium in the regions for which corresponding field specification tag has been specified + computeHydrostaticEquilibrium( domain ); + + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) + { + initializePorosityAndPermeability( mesh, regionNames ); + initializeHydraulicAperture( mesh, regionNames ); + + // Initialize primary variables from applied initial conditions + initializeFluidState( mesh, regionNames ); + + // Initialize the rock thermal quantities: conductivity and solid internal energy + // Note: + // - This must be called after updatePorosityAndPermeability and updatePhaseVolumeFraction + // - This step depends on porosity and phaseVolFraction + if( m_isThermal ) + { + initializeThermalState( mesh, regionNames ); + } + + // Save initial pressure and temperature fields + saveInitialPressureAndTemperature( mesh, regionNames ); + } ); + + // report to the user if some pore volumes are very small + // note: this function is here because: 1) porosity has been initialized and 2) NTG has been applied + validatePoreVolumes( domain ); +} + +void FlowSolverBase::initializePorosityAndPermeability( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) +{ + // Update porosity and permeability + // In addition, to avoid multiplying permeability/porosity bay netToGross in the assembly kernel, we do it once and for all here + mesh.getElemManager().forElementSubRegions< CellElementSubRegion, SurfaceElementSubRegion >( regionNames, [&]( localIndex const, + auto & subRegion ) + { + // Apply netToGross to reference porosity and horizontal permeability + CoupledSolidBase const & porousSolid = + getConstitutiveModel< CoupledSolidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::solidNamesString() ) ); + PermeabilityBase const & permeability = + getConstitutiveModel< PermeabilityBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::permeabilityNamesString() ) ); + arrayView1d< real64 const > const netToGross = subRegion.template getField< fields::flow::netToGross >(); + porousSolid.scaleReferencePorosity( netToGross ); + permeability.scaleHorizontalPermeability( netToGross ); + + // in some initializeState versions it uses newPorosity, so let's run updatePorosityAndPermeability to compute something + saveConvergedState( subRegion ); // necessary for a meaningful porosity update in sequential schemes + updatePorosityAndPermeability( subRegion ); + porousSolid.initializeState(); + + // run final update + saveConvergedState( subRegion ); // necessary for a meaningful porosity update in sequential schemes + updatePorosityAndPermeability( subRegion ); + + // Save the computed porosity into the old porosity + // Note: + // - This must be called after updatePorosityAndPermeability + // - This step depends on porosity + porousSolid.saveConvergedState(); + } ); +} + +void FlowSolverBase::initializeHydraulicAperture( MeshLevel & mesh, const arrayView1d< const string > & regionNames ) +{ + mesh.getElemManager().forElementRegions< SurfaceElementRegion >( regionNames, + [&]( localIndex const, + SurfaceElementRegion & region ) + { + region.forElementSubRegions< SurfaceElementSubRegion >( [&]( SurfaceElementSubRegion & subRegion ) + { subRegion.getWrapper< real64_array >( fields::flow::hydraulicAperture::key()).setApplyDefaultValue( region.getDefaultAperture()); } ); + } ); +} + +void FlowSolverBase::saveInitialPressureAndTemperature( MeshLevel & mesh, const arrayView1d< const string > & regionNames ) +{ + mesh.getElemManager().forElementSubRegions( regionNames, [&]( localIndex const, + ElementSubRegionBase & subRegion ) + { + arrayView1d< real64 const > const pres = subRegion.getField< fields::flow::pressure >(); + arrayView1d< real64 > const initPres = subRegion.getField< fields::flow::initialPressure >(); + arrayView1d< real64 const > const temp = subRegion.getField< fields::flow::temperature >(); + arrayView1d< real64 > const initTemp = subRegion.template getField< fields::flow::initialTemperature >(); + initPres.setValues< parallelDevicePolicy<> >( pres ); + initTemp.setValues< parallelDevicePolicy<> >( temp ); + } ); +} + void FlowSolverBase::updatePorosityAndPermeability( CellElementSubRegion & subRegion ) const { GEOS_MARK_FUNCTION; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp index 6b65e91e979..23897ccb899 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp @@ -100,9 +100,9 @@ class FlowSolverBase : public PhysicsSolverBase */ void updateStencilWeights( DomainPartition & domain ) const; - void enableFixedStressPoromechanicsUpdate(); + void enableFixedStressPoromechanicsUpdate() { m_isFixedStressPoromechanicsUpdate = true; } - void enableJumpStabilization(); + void enableJumpStabilization() { m_isJumpStabilized = true; } void updatePorosityAndPermeability( CellElementSubRegion & subRegion ) const; @@ -114,39 +114,12 @@ class FlowSolverBase : public PhysicsSolverBase */ virtual void saveSequentialIterationState( DomainPartition & domain ) override; - /** - * @brief For each equilibrium initial condition, loop over all the target cells and compute the min/max elevation - * @param[in] domain the domain partition - * @param[in] equilNameToEquilId the map from the name of the initial condition to the initial condition index (used in min/maxElevation) - * @param[out] maxElevation the max elevation for each initial condition - * @param[out] minElevation the min elevation for each initial condition - */ - void findMinMaxElevationInEquilibriumTarget( DomainPartition & domain, // cannot be const... - std::map< string, localIndex > const & equilNameToEquilId, - arrayView1d< real64 > const & maxElevation, - arrayView1d< real64 > const & minElevation ) const; - - /** - * @brief For each source flux boundary condition, loop over all the target cells and sum the owned cells - * @param[in] time the time at the beginning of the time step - * @param[in] dt the time step size - * @param[in] domain the domain partition - * @param[in] bcNameToBcId the map from the name of the boundary condition to the boundary condition index - * @param[out] bcAllSetsSize the total number of owned cells for each source flux boundary condition - */ - void computeSourceFluxSizeScalingFactor( real64 const & time, - real64 const & dt, - DomainPartition & domain, // cannot be const... - std::map< string, localIndex > const & bcNameToBcId, - arrayView1d< globalIndex > const & bcAllSetsSize ) const; - integer & isThermal() { return m_isThermal; } /** * @return The unit in which we evaluate the amount of fluid per element (Mass or Mole). */ - virtual units::Unit getMassUnit() const - { return units::Unit::Mass; } + virtual units::Unit getMassUnit() const { return units::Unit::Mass; } /** * @brief Function to activate the flag allowing negative pressure @@ -173,6 +146,36 @@ class FlowSolverBase : public PhysicsSolverBase real64 const & timeAtBeginningOfStep, real64 const & dt ); + virtual void initializeFluidState( MeshLevel & mesh, const arrayView1d< const string > & regionNames ) { GEOS_UNUSED_VAR( mesh, regionNames ); } + + virtual void initializeThermalState( MeshLevel & mesh, const arrayView1d< const string > & regionNames ) { GEOS_UNUSED_VAR( mesh, regionNames ); } + + /** + * @brief For each equilibrium initial condition, loop over all the target cells and compute the min/max elevation + * @param[in] domain the domain partition + * @param[in] equilNameToEquilId the map from the name of the initial condition to the initial condition index (used in min/maxElevation) + * @param[out] maxElevation the max elevation for each initial condition + * @param[out] minElevation the min elevation for each initial condition + */ + void findMinMaxElevationInEquilibriumTarget( DomainPartition & domain, // cannot be const... + std::map< string, localIndex > const & equilNameToEquilId, + arrayView1d< real64 > const & maxElevation, + arrayView1d< real64 > const & minElevation ) const; + + /** + * @brief For each source flux boundary condition, loop over all the target cells and sum the owned cells + * @param[in] time the time at the beginning of the time step + * @param[in] dt the time step size + * @param[in] domain the domain partition + * @param[in] bcNameToBcId the map from the name of the boundary condition to the boundary condition index + * @param[out] bcAllSetsSize the total number of owned cells for each source flux boundary condition + */ + void computeSourceFluxSizeScalingFactor( real64 const & time, + real64 const & dt, + DomainPartition & domain, // cannot be const... + std::map< string, localIndex > const & bcNameToBcId, + arrayView1d< globalIndex > const & bcAllSetsSize ) const; + protected: /** @@ -207,6 +210,16 @@ class FlowSolverBase : public PhysicsSolverBase virtual void initializePostInitialConditionsPreSubGroups() override; + void initializeState( DomainPartition & domain ); + + virtual void computeHydrostaticEquilibrium( DomainPartition & domain ) { GEOS_UNUSED_VAR( domain ); } + + void initializePorosityAndPermeability( MeshLevel & mesh, arrayView1d< string const > const & regionNames ); + + void initializeHydraulicAperture( MeshLevel & mesh, const arrayView1d< const string > & regionNames ); + + void saveInitialPressureAndTemperature( MeshLevel & mesh, const arrayView1d< const string > & regionNames ); + virtual void setConstitutiveNamesCallSuper( ElementSubRegionBase & subRegion ) const override; /// the number of Degrees of Freedom per cell diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp index 429ea336ab7..812c342de0a 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp @@ -40,8 +40,6 @@ #include "mesh/mpiCommunications/CommunicationTools.hpp" #include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" #include "physicsSolvers/fluidFlow/SinglePhaseBaseFields.hpp" -#include "physicsSolvers/fluidFlow/kernels/singlePhase/AccumulationKernels.hpp" -#include "physicsSolvers/fluidFlow/kernels/singlePhase/ThermalAccumulationKernels.hpp" #include "physicsSolvers/fluidFlow/kernels/singlePhase/MobilityKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/singlePhase/SolutionCheckKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/singlePhase/SolutionScalingKernel.hpp" @@ -406,105 +404,12 @@ void SinglePhaseBase::initializePostInitialConditionsPreSubGroups() DomainPartition & domain = this->getGroupByPath< DomainPartition >( "/Problem/domain" ); - - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, - MeshLevel & mesh, - arrayView1d< string const > const & regionNames ) - { - FieldIdentifiers fieldsToBeSync; - fieldsToBeSync.addElementFields( { fields::flow::pressure::key() }, - regionNames ); - - CommunicationTools::getInstance().synchronizeFields( fieldsToBeSync, mesh, domain.getNeighbors(), false ); - - // Moved the following part from ImplicitStepSetup to here since it only needs to be initialized once - // They will be updated in applySystemSolution and ImplicitStepComplete, respectively - mesh.getElemManager().forElementSubRegions< CellElementSubRegion, SurfaceElementSubRegion >( regionNames, [&]( localIndex const, - auto & subRegion ) - { - // Compute hydrostatic equilibrium in the regions for which corresponding field specification tag has been specified - computeHydrostaticEquilibrium(); - - // 1. update porosity, permeability, and density/viscosity - // In addition, to avoid multiplying permeability/porosity bay netToGross in the assembly kernel, we do it once and for all here - arrayView1d< real64 const > const netToGross = subRegion.template getField< fields::flow::netToGross >(); - CoupledSolidBase const & porousSolid = - getConstitutiveModel< CoupledSolidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::solidNamesString() ) ); - PermeabilityBase const & permeabilityModel = - getConstitutiveModel< PermeabilityBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::permeabilityNamesString() ) ); - permeabilityModel.scaleHorizontalPermeability( netToGross ); - porousSolid.scaleReferencePorosity( netToGross ); - saveConvergedState( subRegion ); // necessary for a meaningful porosity update in sequential schemes - updatePorosityAndPermeability( subRegion ); - - SingleFluidBase const & fluid = - getConstitutiveModel< SingleFluidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::fluidNamesString() ) ); - updateFluidState( subRegion ); - - // 2. save the initial density (for use in the single-phase poromechanics solver to compute the deltaBodyForce) - fluid.initializeState(); - - // 3. save the initial/old porosity - porousSolid.initializeState(); - - // 4. initialize the rock thermal quantities: conductivity and solid internal energy - if( m_isThermal ) - { - // initialized porosity - arrayView2d< real64 const > const porosity = porousSolid.getPorosity(); - - string const & thermalConductivityName = subRegion.template getReference< string >( viewKeyStruct::thermalConductivityNamesString() ); - SinglePhaseThermalConductivityBase const & conductivityMaterial = - getConstitutiveModel< SinglePhaseThermalConductivityBase >( subRegion, thermalConductivityName ); - conductivityMaterial.initializeRockFluidState( porosity ); - // note that there is nothing to update here because thermal conductivity is explicit for now - - updateSolidInternalEnergyModel( subRegion ); - string const & solidInternalEnergyName = subRegion.template getReference< string >( viewKeyStruct::solidInternalEnergyNamesString() ); - SolidInternalEnergy const & solidInternalEnergyMaterial = - getConstitutiveModel< SolidInternalEnergy >( subRegion, solidInternalEnergyName ); - solidInternalEnergyMaterial.saveConvergedState(); - } - } ); - - mesh.getElemManager().forElementRegions< SurfaceElementRegion >( regionNames, - [&]( localIndex const, - SurfaceElementRegion & region ) - { - region.forElementSubRegions< SurfaceElementSubRegion >( [&]( SurfaceElementSubRegion & subRegion ) - { - subRegion.getWrapper< real64_array >( fields::flow::hydraulicAperture::key() ). - setApplyDefaultValue( region.getDefaultAperture() ); - } ); - } ); - - mesh.getElemManager().forElementSubRegions( regionNames, [&]( localIndex const, - ElementSubRegionBase & subRegion ) - { - // Save initial pressure field - arrayView1d< real64 const > const pres = subRegion.getField< fields::flow::pressure >(); - arrayView1d< real64 > const initPres = subRegion.getField< fields::flow::initialPressure >(); - arrayView1d< real64 const > const & temp = subRegion.template getField< fields::flow::temperature >(); - arrayView1d< real64 > const initTemp = subRegion.template getField< fields::flow::initialTemperature >(); - initPres.setValues< parallelDevicePolicy<> >( pres ); - initTemp.setValues< parallelDevicePolicy<> >( temp ); - - // finally update mass and energy - updateMass( subRegion ); - if( m_isThermal ) - updateEnergy( subRegion ); - } ); - } ); - - // report to the user if some pore volumes are very small - // note: this function is here because: 1) porosity has been initialized and 2) NTG has been applied - validatePoreVolumes( domain ); + FlowSolverBase::initializeState( domain ); } -void SinglePhaseBase::computeHydrostaticEquilibrium() +void SinglePhaseBase::computeHydrostaticEquilibrium( DomainPartition & domain ) { FieldSpecificationManager & fsManager = FieldSpecificationManager::getInstance(); - DomainPartition & domain = this->getGroupByPath< DomainPartition >( "/Problem/domain" ); real64 const gravVector[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3( gravityVector() ); @@ -678,6 +583,48 @@ void SinglePhaseBase::computeHydrostaticEquilibrium() } ); } +void SinglePhaseBase::initializeFluidState( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) +{ + mesh.getElemManager().forElementSubRegions< CellElementSubRegion, SurfaceElementSubRegion >( regionNames, [&]( localIndex const, + auto & subRegion ) + { + SingleFluidBase const & fluid = + getConstitutiveModel< SingleFluidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::fluidNamesString())); + updateFluidState( subRegion ); + + // 2. save the initial density (for use in the single-phase poromechanics solver to compute the deltaBodyForce) + fluid.initializeState(); + + updateMass( subRegion ); + } ); +} + +void SinglePhaseBase::initializeThermalState( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) +{ + mesh.getElemManager().forElementSubRegions< CellElementSubRegion, SurfaceElementSubRegion >( regionNames, [&]( localIndex const, + auto & subRegion ) + { + // initialized porosity + CoupledSolidBase const & porousSolid = + getConstitutiveModel< CoupledSolidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::solidNamesString() ) ); + arrayView2d< real64 const > const porosity = porousSolid.getPorosity(); + + string const & thermalConductivityName = subRegion.template getReference< string >( viewKeyStruct::thermalConductivityNamesString()); + SinglePhaseThermalConductivityBase const & conductivityMaterial = + getConstitutiveModel< SinglePhaseThermalConductivityBase >( subRegion, thermalConductivityName ); + conductivityMaterial.initializeRockFluidState( porosity ); + // note that there is nothing to update here because thermal conductivity is explicit for now + + updateSolidInternalEnergyModel( subRegion ); + string const & solidInternalEnergyName = subRegion.template getReference< string >( viewKeyStruct::solidInternalEnergyNamesString()); + SolidInternalEnergy const & solidInternalEnergyMaterial = + getConstitutiveModel< SolidInternalEnergy >( subRegion, solidInternalEnergyName ); + solidInternalEnergyMaterial.saveConvergedState(); + + updateEnergy( subRegion ); + } ); +} + void SinglePhaseBase::implicitStepSetup( real64 const & GEOS_UNUSED_PARAM( time_n ), real64 const & GEOS_UNUSED_PARAM( dt ), DomainPartition & domain ) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp index bf1cd19eadd..9e40ecda15f 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp @@ -335,10 +335,14 @@ class SinglePhaseBase : public FlowSolverBase virtual void initializePostInitialConditionsPreSubGroups() override; + virtual void initializeFluidState( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; + + virtual void initializeThermalState( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; + /** * @brief Compute the hydrostatic equilibrium using the compositions and temperature input tables */ - void computeHydrostaticEquilibrium(); + virtual void computeHydrostaticEquilibrium( DomainPartition & domain ) override; /** * @brief Update the cell-wise pressure gradient diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp index b87ce6bcf75..760394bb02b 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp @@ -33,8 +33,6 @@ #include "fieldSpecification/AquiferBoundaryCondition.hpp" #include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" #include "physicsSolvers/fluidFlow/SinglePhaseBaseFields.hpp" -#include "physicsSolvers/fluidFlow/kernels/singlePhase/AccumulationKernels.hpp" -#include "physicsSolvers/fluidFlow/kernels/singlePhase/ThermalAccumulationKernels.hpp" #include "physicsSolvers/fluidFlow/kernels/singlePhase/ResidualNormKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/singlePhase/FluxComputeKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/singlePhase/ThermalFluxComputeKernel.hpp" diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/CompositionalMultiphaseHybridFVMKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/CompositionalMultiphaseHybridFVMKernels.hpp index c388cad6da5..ea1f4ec0d12 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/CompositionalMultiphaseHybridFVMKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/CompositionalMultiphaseHybridFVMKernels.hpp @@ -33,7 +33,6 @@ #include "physicsSolvers/PhysicsSolverBaseKernels.hpp" #include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp" #include "physicsSolvers/fluidFlow/StencilAccessors.hpp" -#include "physicsSolvers/fluidFlow/kernels/compositional/AccumulationKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/PropertyKernelBase.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/KernelLaunchSelectors.hpp" diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/DiffusionDispersionFluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/DiffusionDispersionFluxComputeKernel.hpp index ac5ab4265da..e8ea185c993 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/DiffusionDispersionFluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/DiffusionDispersionFluxComputeKernel.hpp @@ -33,7 +33,6 @@ #include "mesh/ElementRegionManager.hpp" #include "physicsSolvers/fluidFlow/CompositionalMultiphaseUtilities.hpp" #include "physicsSolvers/fluidFlow/StencilAccessors.hpp" -#include "physicsSolvers/fluidFlow/kernels/compositional/KernelLaunchSelectors.hpp" namespace geos { diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/SolutionScalingKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/SolutionScalingKernel.hpp index 399ccaaa6bc..7cdbd3456cb 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/SolutionScalingKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/SolutionScalingKernel.hpp @@ -21,6 +21,7 @@ #define GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_SOLUTIONSCALINGKERNEL_HPP #include "physicsSolvers/fluidFlow/kernels/compositional/SolutionScalingAndCheckingKernelBase.hpp" +#include "physicsSolvers/fluidFlow/kernels/compositional/AccumulationKernel.hpp" #include "physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp" namespace geos diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/SinglePhaseHybridFVMKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/SinglePhaseHybridFVMKernels.hpp index 13e5eeabfcb..124e454d90f 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/SinglePhaseHybridFVMKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/SinglePhaseHybridFVMKernels.hpp @@ -37,7 +37,6 @@ #include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" #include "physicsSolvers/fluidFlow/SinglePhaseBaseFields.hpp" #include "physicsSolvers/fluidFlow/StencilAccessors.hpp" -#include "physicsSolvers/fluidFlow/kernels/singlePhase/AccumulationKernels.hpp" #include "physicsSolvers/fluidFlow/kernels/HybridFVMHelperKernels.hpp" #include "physicsSolvers/PhysicsSolverBaseKernels.hpp" #include "codingUtilities/Utilities.hpp" diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp index 95f59c179ec..dbd4637d01b 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp @@ -42,8 +42,6 @@ #include "physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.hpp" #include "physicsSolvers/fluidFlow/wells/kernels/ThermalCompositionalMultiphaseWellKernels.hpp" #include "physicsSolvers/fluidFlow/wells/kernels/PerforationFluxKernels.hpp" -#include "physicsSolvers/fluidFlow/kernels/compositional/AccumulationKernel.hpp" -#include "physicsSolvers/fluidFlow/kernels/compositional/ThermalAccumulationKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/SolutionScalingKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/ThermalSolutionScalingKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/SolutionCheckKernel.hpp" @@ -62,7 +60,6 @@ namespace geos using namespace dataRepository; using namespace constitutive; -using namespace compositionalMultiphaseWellKernels; CompositionalMultiphaseWell::CompositionalMultiphaseWell( const string & name, Group * const parent ) @@ -910,7 +907,8 @@ void CompositionalMultiphaseWell::updateTotalMassDensity( WellElementSubRegion & subRegion, fluid ) : - TotalMassDensityKernelFactory:: + compositionalMultiphaseWellKernels:: + TotalMassDensityKernelFactory:: createAndLaunch< parallelDevicePolicy<> >( m_numComponents, m_numPhases, subRegion, @@ -976,8 +974,10 @@ void CompositionalMultiphaseWell::initializeWells( DomainPartition & domain, rea { ElementRegionManager & elemManager = mesh.getElemManager(); - PresTempCompFracInitializationKernel::CompFlowAccessors resCompFlowAccessors( mesh.getElemManager(), flowSolver.getName() ); - PresTempCompFracInitializationKernel::MultiFluidAccessors resMultiFluidAccessors( mesh.getElemManager(), flowSolver.getName() ); + compositionalMultiphaseWellKernels::PresTempCompFracInitializationKernel::CompFlowAccessors + resCompFlowAccessors( mesh.getElemManager(), flowSolver.getName() ); + compositionalMultiphaseWellKernels::PresTempCompFracInitializationKernel::MultiFluidAccessors + resMultiFluidAccessors( mesh.getElemManager(), flowSolver.getName() ); elemManager.forElementSubRegions< WellElementSubRegion >( regionNames, [&]( localIndex const, @@ -1011,7 +1011,8 @@ void CompositionalMultiphaseWell::initializeWells( DomainPartition & domain, rea // 1) Loop over all perforations to compute an average mixture density and component fraction // 2) Initialize the reference pressure // 3) Estimate the pressures in the well elements using the average density - PresTempCompFracInitializationKernel:: + compositionalMultiphaseWellKernels:: + PresTempCompFracInitializationKernel:: launch( perforationData.size(), subRegion.size(), numComp, @@ -1053,11 +1054,12 @@ void CompositionalMultiphaseWell::initializeWells( DomainPartition & domain, rea wellElemCompFrac ); } ); - CompDensInitializationKernel::launch( subRegion.size(), - numComp, - wellElemCompFrac, - wellElemTotalDens, - wellElemCompDens ); + compositionalMultiphaseWellKernels:: + CompDensInitializationKernel::launch( subRegion.size(), + numComp, + wellElemCompFrac, + wellElemTotalDens, + wellElemCompDens ); // 5) Recompute the pressure-dependent properties updateSubRegionState( subRegion ); @@ -1310,7 +1312,7 @@ CompositionalMultiphaseWell::calculateResidualNorm( real64 const & time_n, else { real64 subRegionResidualNorm[1]{}; - ResidualNormKernelFactory:: + compositionalMultiphaseWellKernels::ResidualNormKernelFactory:: createAndLaunch< parallelDevicePolicy<> >( m_numComponents, numDofPerWellElement(), m_targetPhaseIndex, @@ -1914,24 +1916,25 @@ void CompositionalMultiphaseWell::assemblePressureRelations( real64 const & time bool controlHasSwitched = false; isothermalCompositionalMultiphaseBaseKernels:: - KernelLaunchSelectorCompTherm< PressureRelationKernel >( numFluidComponents(), - isThermal, - subRegion.size(), - dofManager.rankOffset(), - subRegion.isLocallyOwned(), - subRegion.getTopWellElementIndex(), - m_targetPhaseIndex, - wellControls, - time_n + dt, // controls evaluated with BHP/rate of the end of step - wellElemDofNumber, - wellElemGravCoef, - nextWellElemIndex, - wellElemPres, - wellElemTotalMassDens, - dWellElemTotalMassDens, - controlHasSwitched, - localMatrix, - localRhs ); + KernelLaunchSelectorCompTherm< compositionalMultiphaseWellKernels::PressureRelationKernel > + ( numFluidComponents(), + isThermal, + subRegion.size(), + dofManager.rankOffset(), + subRegion.isLocallyOwned(), + subRegion.getTopWellElementIndex(), + m_targetPhaseIndex, + wellControls, + time_n + dt, // controls evaluated with BHP/rate of the end of step + wellElemDofNumber, + wellElemGravCoef, + nextWellElemIndex, + wellElemPres, + wellElemTotalMassDens, + dWellElemTotalMassDens, + controlHasSwitched, + localMatrix, + localRhs ); if( controlHasSwitched ) { diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.hpp index 1ca6ebac997..1995119e9f4 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.hpp @@ -33,7 +33,6 @@ #include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp" #include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" #include "physicsSolvers/fluidFlow/StencilAccessors.hpp" -#include "physicsSolvers/fluidFlow/kernels/compositional/AccumulationKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/PropertyKernelBase.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/SolutionScalingKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/SolutionCheckKernel.hpp" diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.cpp index 059259cb12b..a3c8a6cc967 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.cpp @@ -77,7 +77,7 @@ void AcousticWaveEquationSEM::registerDataOnMesh( Group & meshBodies ) nodeManager.registerField< acousticfields::Pressure_nm1, acousticfields::Pressure_n, acousticfields::Pressure_np1, - acousticfields::PressureDoubleDerivative, + acousticfields::PressureForward, acousticfields::ForcingRHS, acousticfields::AcousticMassVector, acousticfields::DampingVector, @@ -106,6 +106,7 @@ void AcousticWaveEquationSEM::registerDataOnMesh( Group & meshBodies ) subRegion.registerField< acousticfields::AcousticVelocity >( getName() ); subRegion.registerField< acousticfields::AcousticDensity >( getName() ); subRegion.registerField< acousticfields::PartialGradient >( getName() ); + subRegion.registerField< acousticfields::PartialGradient2 >( getName() ); } ); } ); @@ -298,6 +299,8 @@ void AcousticWaveEquationSEM::initializePostInitialConditionsPreSubGroups() arrayView1d< real32 > grad = elementSubRegion.getField< acousticfields::PartialGradient >(); grad.zero(); + arrayView1d< real32 > grad2 = elementSubRegion.getField< acousticfields::PartialGradient2 >(); + grad2.zero(); finiteElement::FiniteElementDispatchHandler< SEM_FE_TYPES >::dispatch3D( fe, [&] ( auto const finiteElement ) { @@ -881,43 +884,35 @@ real64 AcousticWaveEquationSEM::explicitStepForward( real64 const & time_n, { NodeManager & nodeManager = mesh.getNodeManager(); - arrayView1d< real32 > const p_nm1 = nodeManager.getField< acousticfields::Pressure_nm1 >(); arrayView1d< real32 > const p_n = nodeManager.getField< acousticfields::Pressure_n >(); - arrayView1d< real32 > const p_np1 = nodeManager.getField< acousticfields::Pressure_np1 >(); if( computeGradient && cycleNumber >= 0 ) { - arrayView1d< real32 > const p_dt2 = nodeManager.getField< acousticfields::PressureDoubleDerivative >(); - if( m_enableLifo ) { if( !m_lifo ) { int const rank = MpiWrapper::commRank( MPI_COMM_GEOS ); - std::string lifoPrefix = GEOS_FMT( "lifo/rank_{:05}/pdt2_shot{:06}", rank, m_shotIndex ); - m_lifo = std::make_unique< LifoStorage< real32, localIndex > >( lifoPrefix, p_dt2, m_lifoOnDevice, m_lifoOnHost, m_lifoSize ); + std::string lifoPrefix = GEOS_FMT( "lifo/rank_{:05}/p_forward_shot{:06}", rank, m_shotIndex ); + m_lifo = std::make_unique< LifoStorage< real32, localIndex > >( lifoPrefix, p_n, m_lifoOnDevice, m_lifoOnHost, m_lifoSize ); } m_lifo->pushWait(); } - forAll< EXEC_POLICY >( nodeManager.size(), [=] GEOS_HOST_DEVICE ( localIndex const nodeIdx ) - { - p_dt2[nodeIdx] = (p_np1[nodeIdx] - 2*p_n[nodeIdx] + p_nm1[nodeIdx]) / pow( dt, 2 ); - } ); if( m_enableLifo ) { // Need to tell LvArray data is on GPU to avoir HtoD copy - p_dt2.move( LvArray::MemorySpace::cuda, false ); - m_lifo->pushAsync( p_dt2 ); + p_n.move( LvArray::MemorySpace::cuda, false ); + m_lifo->pushAsync( p_n ); } else { GEOS_MARK_SCOPE ( DirectWrite ); - p_dt2.move( LvArray::MemorySpace::host, false ); + p_n.move( LvArray::MemorySpace::host, false ); int const rank = MpiWrapper::commRank( MPI_COMM_GEOS ); - std::string fileName = GEOS_FMT( "lifo/rank_{:05}/pressuredt2_{:06}_{:08}.dat", rank, m_shotIndex, cycleNumber ); + std::string fileName = GEOS_FMT( "lifo/rank_{:05}/pressure_forward_{:06}_{:08}.dat", rank, m_shotIndex, cycleNumber ); int lastDirSeparator = fileName.find_last_of( "/\\" ); std::string dirName = fileName.substr( 0, lastDirSeparator ); if( string::npos != (size_t)lastDirSeparator && !directoryExists( dirName )) @@ -929,7 +924,7 @@ real64 AcousticWaveEquationSEM::explicitStepForward( real64 const & time_n, GEOS_THROW_IF( !wf, getDataContext() << ": Could not open file "<< fileName << " for writing", InputError ); - wf.write( (char *)&p_dt2[0], p_dt2.size()*sizeof( real32 ) ); + wf.write( (char *)&p_n[0], p_n.size()*sizeof( real32 ) ); wf.close( ); GEOS_THROW_IF( !wf.good(), getDataContext() << ": An error occured while writing "<< fileName, @@ -959,12 +954,18 @@ real64 AcousticWaveEquationSEM::explicitStepBackward( real64 const & time_n, { NodeManager & nodeManager = mesh.getNodeManager(); - arrayView1d< real32 const > const mass = nodeManager.getField< acousticfields::AcousticMassVector >(); - arrayView1d< real32 > const p_nm1 = nodeManager.getField< acousticfields::Pressure_nm1 >(); arrayView1d< real32 > const p_n = nodeManager.getField< acousticfields::Pressure_n >(); arrayView1d< real32 > const p_np1 = nodeManager.getField< acousticfields::Pressure_np1 >(); + //// Compute q_dt2 and store it in p_nm1 + 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_np1[a] - 2*p_n[a] + p_nm1[a]) / pow( dt, 2 ); + } ); + EventManager const & event = getGroupByPath< EventManager >( "/Problem/Events" ); real64 const & maxTime = event.getReference< real64 >( EventManager::viewKeyStruct::maxTimeString() ); int const maxCycle = int(round( maxTime / dt )); @@ -973,11 +974,11 @@ real64 AcousticWaveEquationSEM::explicitStepBackward( real64 const & time_n, { ElementRegionManager & elemManager = mesh.getElemManager(); - arrayView1d< real32 > const p_dt2 = nodeManager.getField< acousticfields::PressureDoubleDerivative >(); + arrayView1d< real32 > const p_forward = nodeManager.getField< acousticfields::PressureForward >(); if( m_enableLifo ) { - m_lifo->pop( p_dt2 ); + m_lifo->pop( p_forward ); if( m_lifo->empty() ) delete m_lifo.release(); } @@ -986,37 +987,60 @@ real64 AcousticWaveEquationSEM::explicitStepBackward( real64 const & time_n, GEOS_MARK_SCOPE ( DirectRead ); int const rank = MpiWrapper::commRank( MPI_COMM_GEOS ); - std::string fileName = GEOS_FMT( "lifo/rank_{:05}/pressuredt2_{:06}_{:08}.dat", rank, m_shotIndex, cycleNumber ); + std::string fileName = GEOS_FMT( "lifo/rank_{:05}/pressure_forward_{:06}_{:08}.dat", rank, m_shotIndex, cycleNumber ); std::ifstream wf( fileName, std::ios::in | std::ios::binary ); GEOS_THROW_IF( !wf, getDataContext() << ": Could not open file "<< fileName << " for reading", InputError ); - p_dt2.move( LvArray::MemorySpace::host, true ); - wf.read( (char *)&p_dt2[0], p_dt2.size()*sizeof( real32 ) ); + p_forward.move( LvArray::MemorySpace::host, true ); + wf.read( (char *)&p_forward[0], p_forward.size()*sizeof( real32 ) ); wf.close( ); remove( fileName.c_str() ); } elemManager.forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, CellElementSubRegion & elementSubRegion ) { - arrayView1d< real32 const > const velocity = elementSubRegion.getField< acousticfields::AcousticVelocity >(); + arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords = nodeManager.getField< fields::referencePosition32 >().toViewConst(); arrayView1d< real32 > grad = elementSubRegion.getField< acousticfields::PartialGradient >(); + arrayView1d< real32 > grad2 = elementSubRegion.getField< acousticfields::PartialGradient2 >(); arrayView2d< localIndex const, cells::NODE_MAP_USD > const & elemsToNodes = elementSubRegion.nodeList(); - constexpr localIndex numNodesPerElem = 8; arrayView1d< integer const > const elemGhostRank = elementSubRegion.ghostRank(); GEOS_MARK_SCOPE ( updatePartialGradient ); - forAll< EXEC_POLICY >( elementSubRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const eltIdx ) + + //COMPUTE GRADIENTS with respect to K=1/rho*c2 (grad) and b=1/rho (grad2) + finiteElement::FiniteElementBase const & + fe = elementSubRegion.getReference< finiteElement::FiniteElementBase >( getDiscretizationName() ); + + finiteElement::FiniteElementDispatchHandler< SEM_FE_TYPES >::dispatch3D( fe, [&] ( auto const finiteElement ) { - if( elemGhostRank[eltIdx]<0 ) - { - for( localIndex i = 0; i < numNodesPerElem; ++i ) - { - localIndex nodeIdx = elemsToNodes[eltIdx][i]; - grad[eltIdx] += (-2/velocity[eltIdx]) * mass[nodeIdx]/8.0 * (p_dt2[nodeIdx] * p_n[nodeIdx]); - } - } + using FE_TYPE = TYPEOFREF( finiteElement ); + + AcousticMatricesSEM::GradientKappaBuoyancy< FE_TYPE > kernelG( finiteElement ); + kernelG.template computeGradient< EXEC_POLICY, ATOMIC_POLICY >( elementSubRegion.size(), + nodeCoords, + elemsToNodes, + elemGhostRank, + p_nm1, + p_n, + p_forward, + grad, + grad2 ); + + } ); + + // // Change of variables to return grad with respect to c and rho + // arrayView1d< real32 const > const velocity = elementSubRegion.getField< acousticfields::AcousticVelocity >(); + // arrayView1d< real32 const > const density = elementSubRegion.getField< acousticfields::AcousticDensity >(); + // forAll< EXEC_POLICY >( elementSubRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const eltIdx ) + // { + // if( elemGhostRank[eltIdx]<0 ) + // { + // grad2[eltIdx] = -1/(pow(density[eltIdx]*velocity[eltIdx],2)) * grad[eltIdx] - 1/pow(density[eltIdx],2) * grad2[eltIdx]; + // grad[eltIdx]= -2/(density[eltIdx]*pow(velocity[eltIdx],3)) * grad[eltIdx]; + // } + // } ); } ); } diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticFields.hpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticFields.hpp index dd067a11a1f..de61245599c 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticFields.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticFields.hpp @@ -59,13 +59,13 @@ DECLARE_FIELD( Pressure_np1, WRITE_AND_READ, "Scalar pressure at time n+1." ); -DECLARE_FIELD( PressureDoubleDerivative, - "pressureDoubleDerivative", +DECLARE_FIELD( PressureForward, + "pressureForward", array1d< real32 >, 0, NOPLOT, WRITE_AND_READ, - "Double derivative of the pressure for each node to compute the gradient" ); + "Pressure field from forward pass on each node to compute the gradient" ); DECLARE_FIELD( Velocity_x, "velocity_x", @@ -99,6 +99,14 @@ DECLARE_FIELD( PartialGradient, WRITE_AND_READ, "Partiel gradient computed during backward propagation" ); +DECLARE_FIELD( PartialGradient2, + "partialGradient2", + array1d< real32 >, + 0, + NOPLOT, + WRITE_AND_READ, + "Partial gradient for density/velocity computed during backward propagation" ); + DECLARE_FIELD( ForcingRHS, "rhs", array1d< real32 >, diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticMatricesSEMKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticMatricesSEMKernel.hpp index b311d04d45d..9d3d14fa31c 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticMatricesSEMKernel.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticMatricesSEMKernel.hpp @@ -152,7 +152,70 @@ struct AcousticMatricesSEM }; + template< typename FE_TYPE > + struct GradientKappaBuoyancy + { + + GradientKappaBuoyancy( FE_TYPE const & finiteElement ) + : m_finiteElement( finiteElement ) + {} + + /** + * @brief Launch the computation of the 2 gradients relative to the coeff of the wave equation K=1/rho*c2 and b=1/rho + * @tparam EXEC_POLICY the execution policy + * @tparam ATOMIC_POLICY the atomic policy + * @param[in] size the number of cells in the subRegion + * @param[in] nodeCoords coordinates of the nodes + * @param[in] elemsToNodes map from element to nodes + * @param[in] q_dt2 second order derivative in time of backward + * @param[in] q_n current time step of backward + * @param[in] p_n current time step of forward + * @param[out] grad first part of gradient vector with respect to K=1/rho*c2 + * @param[out] grad2 second part of gradient vector with respact to b=1/rho + */ + template< typename EXEC_POLICY, typename ATOMIC_POLICY > + void + computeGradient( localIndex const size, + arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords, + arrayView2d< localIndex const, cells::NODE_MAP_USD > const elemsToNodes, + arrayView1d< integer const > const elemGhostRank, + arrayView1d< real32 const > const q_dt2, + arrayView1d< real32 const > const q_n, + arrayView1d< real32 const > const p_n, + arrayView1d< real32 > const grad, + arrayView1d< real32 > const grad2 ) + { + forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const e ) + { + if( elemGhostRank[e]<0 ) + { + // only the eight corners of the mesh cell are needed to compute the Jacobian + real64 xLocal[ 8 ][ 3 ]; + for( localIndex a = 0; a < 8; ++a ) + { + localIndex const nodeIndex = elemsToNodes( e, FE_TYPE::meshIndexToLinearIndex3D( a ) ); + for( localIndex i = 0; i < 3; ++i ) + { + xLocal[a][i] = nodeCoords( nodeIndex, i ); + } + } + constexpr localIndex numQuadraturePointsPerElem = FE_TYPE::numQuadraturePoints; + for( localIndex q = 0; q < numQuadraturePointsPerElem; ++q ) + { + localIndex nodeIdx = elemsToNodes( e, q ); + grad[e] += q_dt2[nodeIdx] * p_n[nodeIdx] * m_finiteElement.computeMassTerm( q, xLocal ); + m_finiteElement.template computeStiffnessTerm( q, xLocal, [&] ( const int i, const int j, const real64 val ) + { + grad2[e] += val* q_n[elemsToNodes( e, j )] * p_n[elemsToNodes( e, i )]; + } ); + } + } + } ); // end loop over element + } + /// The finite element space/discretization object for the element type in the subRegion + FE_TYPE const & m_finiteElement; + }; }; diff --git a/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.cpp b/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.cpp index c5baf9c7983..f1ac14c24ea 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.cpp @@ -484,9 +484,19 @@ void WaveSolverBase::computeTargetNodeSet( arrayView2d< localIndex const, cells: void WaveSolverBase::incrementIndexSeismoTrace( real64 const time_n ) { - while( (m_dtSeismoTrace * m_indexSeismoTrace) <= (time_n + epsilonLoc) && m_indexSeismoTrace < m_nsamplesSeismoTrace ) + if( m_forward ) + { + while( (m_dtSeismoTrace * m_indexSeismoTrace) <= (time_n + epsilonLoc) && m_indexSeismoTrace < m_nsamplesSeismoTrace ) + { + m_indexSeismoTrace++; + } + } + else { - m_indexSeismoTrace++; + while( (m_dtSeismoTrace * m_indexSeismoTrace) >= (time_n - epsilonLoc) && m_indexSeismoTrace > 0 ) + { + m_indexSeismoTrace--; + } } } @@ -515,12 +525,14 @@ void WaveSolverBase::computeAllSeismoTraces( real64 const time_n, if( m_nsamplesSeismoTrace == 0 ) return; integer const dir = m_forward ? +1 : -1; - for( localIndex iSeismo = m_indexSeismoTrace; iSeismo < m_nsamplesSeismoTrace; iSeismo++ ) + integer const beginIndex = m_forward ? m_indexSeismoTrace : m_nsamplesSeismoTrace-m_indexSeismoTrace; + for( localIndex iSeismo = beginIndex; iSeismo < m_nsamplesSeismoTrace; iSeismo++ ) { - real64 const timeSeismo = m_dtSeismoTrace * (m_forward ? iSeismo : (m_nsamplesSeismoTrace - 1) - iSeismo); - if( dir * timeSeismo > dir * (time_n + epsilonLoc) ) + localIndex seismoIndex = m_forward ? iSeismo : m_nsamplesSeismoTrace-iSeismo; + real64 const timeSeismo = m_dtSeismoTrace * seismoIndex; + if( dir * timeSeismo > dir * time_n + epsilonLoc ) break; - WaveSolverUtils::computeSeismoTrace( time_n, dir * dt, timeSeismo, iSeismo, m_receiverNodeIds, + WaveSolverUtils::computeSeismoTrace( time_n, dir * dt, timeSeismo, seismoIndex, m_receiverNodeIds, m_receiverConstants, m_receiverIsLocal, var_np1, var_n, varAtReceivers, coeffs, add ); } } @@ -535,12 +547,14 @@ void WaveSolverBase::compute2dVariableAllSeismoTraces( localIndex const regionIn if( m_nsamplesSeismoTrace == 0 ) return; integer const dir = m_forward ? +1 : -1; - for( localIndex iSeismo = m_indexSeismoTrace; iSeismo < m_nsamplesSeismoTrace; iSeismo++ ) + integer const beginIndex = m_forward ? m_indexSeismoTrace : m_nsamplesSeismoTrace-m_indexSeismoTrace; + for( localIndex iSeismo = beginIndex; iSeismo < m_nsamplesSeismoTrace; iSeismo++ ) { - real64 const timeSeismo = m_dtSeismoTrace * (m_forward ? iSeismo : (m_nsamplesSeismoTrace - 1) - iSeismo); - if( dir * timeSeismo > dir * (time_n + epsilonLoc)) + localIndex seismoIndex = m_forward ? iSeismo : m_nsamplesSeismoTrace-iSeismo; + real64 const timeSeismo = m_dtSeismoTrace * seismoIndex; + if( dir * timeSeismo > dir * time_n + epsilonLoc ) break; - WaveSolverUtils::compute2dVariableSeismoTrace( time_n, dir * dt, regionIndex, m_receiverRegion, timeSeismo, iSeismo, m_receiverElem, + WaveSolverUtils::compute2dVariableSeismoTrace( time_n, dir * dt, regionIndex, m_receiverRegion, timeSeismo, seismoIndex, m_receiverElem, m_receiverConstants, m_receiverIsLocal, var_np1, var_n, varAtReceivers ); } } diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd index 394f021a29e..7a6e7967c03 100644 --- a/src/coreComponents/schema/schema.xsd +++ b/src/coreComponents/schema/schema.xsd @@ -1761,6 +1761,8 @@ stress - traction is applied to the faces as specified by the inner product of i + + @@ -1781,8 +1783,6 @@ stress - traction is applied to the faces as specified by the inner product of i - - diff --git a/src/coreComponents/schema/schema.xsd.other b/src/coreComponents/schema/schema.xsd.other index 6488245af14..33386f19956 100644 --- a/src/coreComponents/schema/schema.xsd.other +++ b/src/coreComponents/schema/schema.xsd.other @@ -1027,7 +1027,7 @@ - + @@ -1282,23 +1282,6 @@ - - - - - - - - - - - - - - - - - diff --git a/src/coreComponents/unitTests/wavePropagationTests/CMakeLists.txt b/src/coreComponents/unitTests/wavePropagationTests/CMakeLists.txt index fab0a741e73..c1296f06b10 100644 --- a/src/coreComponents/unitTests/wavePropagationTests/CMakeLists.txt +++ b/src/coreComponents/unitTests/wavePropagationTests/CMakeLists.txt @@ -6,7 +6,9 @@ set( gtest_geosx_tests testWavePropagationDAS.cpp testWavePropagationElasticVTI.cpp testWavePropagationAttenuation.cpp - testWavePropagationAcousticFirstOrder.cpp ) + testWavePropagationAcousticFirstOrder.cpp + testWavePropagationAdjoint1.cpp + ) set( tplDependencyList ${parallelDeps} gtest ) diff --git a/src/coreComponents/unitTests/wavePropagationTests/testWavePropagationAdjoint1.cpp b/src/coreComponents/unitTests/wavePropagationTests/testWavePropagationAdjoint1.cpp new file mode 100644 index 00000000000..68c7618fecc --- /dev/null +++ b/src/coreComponents/unitTests/wavePropagationTests/testWavePropagationAdjoint1.cpp @@ -0,0 +1,404 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 Total, S.A + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +// using some utility classes from the following unit test +#include "unitTests/fluidFlowTests/testCompFlowUtils.hpp" + +#include "common/DataTypes.hpp" +#include "mainInterface/initialization.hpp" +#include "mainInterface/ProblemManager.hpp" +#include "mesh/DomainPartition.hpp" +#include "mainInterface/GeosxState.hpp" +#include "physicsSolvers/PhysicsSolverManager.hpp" +#include "physicsSolvers/wavePropagation/shared/WaveSolverBase.hpp" +#include "physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.hpp" + +#include + +using namespace geos; +using namespace geos::dataRepository; +using namespace geos::testing; + +CommandLineOptions g_commandLineOptions; + +// This unit test checks the interpolation done to extract seismic traces from a wavefield. +// It computes a seismogram at a receiver co-located with the source and compares it to the surrounding receivers. +char const * xmlInput = + R"xml( + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + )xml"; + +class AcousticWaveEquationSEMTest : public ::testing::Test +{ +public: + + AcousticWaveEquationSEMTest(): + state( std::make_unique< CommandLineOptions >( g_commandLineOptions ) ) + {} + +protected: + + void SetUp() override + { + setupProblemFromXML( state.getProblemManager(), xmlInput ); + } + + static real64 constexpr time = 0.0; + static real64 constexpr dt = 5e-3; + static real64 constexpr eps = std::numeric_limits< real64 >::epsilon(); + + GeosxState state; + AcousticWaveEquationSEM * propagator; +}; + +real64 constexpr AcousticWaveEquationSEMTest::time; +real64 constexpr AcousticWaveEquationSEMTest::dt; +real64 constexpr AcousticWaveEquationSEMTest::eps; + +TEST_F( AcousticWaveEquationSEMTest, SeismoTrace ) +{ + + DomainPartition & domain = state.getProblemManager().getDomainPartition(); + propagator = &state.getProblemManager().getPhysicsSolverManager().getGroup< AcousticWaveEquationSEM >( "acousticSolver" ); + + // Check source term (sourceCoordinates and sourceValue) + array2d< real32 > rhsForward; + rhsForward.resize( 51, 1 ); + real32 * ptrTimeSourceFrequency = &propagator->getReference< real32 >( AcousticWaveEquationSEM::viewKeyStruct::timeSourceFrequencyString() ); + real32 * ptrTimeSourceDelay = &propagator->getReference< real32 >( AcousticWaveEquationSEM::viewKeyStruct::timeSourceDelayString() ); + localIndex * ptrRickerOrder = &propagator->getReference< localIndex >( AcousticWaveEquationSEM::viewKeyStruct::rickerOrderString() ); + + real64 time_n = time; + std::cout << "Begin forward:" << time_n << std::endl; + // run for 0.25s (100 steps) + for( int i=0; i<50; i++ ) + { + rhsForward[i][0]=WaveSolverUtils::evaluateRicker( time_n, *ptrTimeSourceFrequency, *ptrTimeSourceDelay, *ptrRickerOrder ); + propagator->explicitStepForward( time_n, dt, i, domain, false ); + time_n += dt; + } + // cleanup (triggers calculation of the remaining seismograms data points) + propagator->cleanup( 1.0, 50, 0, 0, domain ); + + // retrieve seismo + arrayView2d< real32 > const pReceivers = propagator->getReference< array2d< real32 > >( AcousticWaveEquationSEM::viewKeyStruct::pressureNp1AtReceiversString() ).toView(); + + // move it to CPU, if needed + pReceivers.move( hostMemorySpace, false ); + + // check number of seismos and trace length + ASSERT_EQ( pReceivers.size( 1 ), 5 ); + ASSERT_EQ( pReceivers.size( 0 ), 51 ); + + /*----------Save receiver forward----------------------*/ + array2d< real32 > uForward; + uForward.resize( 51, 1 ); + + // save receiver value forward on uForward. + for( int i = 0; i < 51; i++ ) + { + /*std::cout << "time: " << i*dt << std::endl; + std::cout << "pReceivers1 " << i << ":" << pReceivers[i][0] << std::endl; + std::cout << "pReceivers2 " << i << ":" << pReceivers[i][1] << std::endl; + std::cout << "pReceivers3 " << i << ":" << pReceivers[i][2] << std::endl; + std::cout << "pReceivers4 " << i << ":" << pReceivers[i][3] << std::endl; + std::cout << "rhsForward " << i << ":" << rhsForward[i][0] << std::endl;*/ + uForward[i][0] = pReceivers[i][0]; + pReceivers[i][0] = 0.; + pReceivers[i][1] = 0.; + pReceivers[i][2] = 0.; + pReceivers[i][3] = 0.; + } + + ASSERT_EQ( rhsForward.size( 1 ), 1 ); + ASSERT_EQ( rhsForward.size( 0 ), 51 ); + + arrayView2d< localIndex > const rNodeIds = propagator->getReference< array2d< localIndex > >( AcousticWaveEquationSEM::viewKeyStruct::receiverNodeIdsString() ).toView(); + rNodeIds.move( hostMemorySpace, false ); + localIndex sNodesIdsAfterModif=rNodeIds[0][0]; + std::cout << "ref back sNodeIds[0][0]:" << sNodesIdsAfterModif << std::endl; + + /*---------------------------------------------------*/ + + std::cout << "Begin backward:" << time_n << std::endl; + + //----------Switch source and receiver1 position for backward----------------------// + arrayView2d< real64 > const sCoord = propagator->getReference< array2d< real64 > >( AcousticWaveEquationSEM::viewKeyStruct::sourceCoordinatesString() ).toView(); + arrayView2d< real64 > const rCoord = propagator->getReference< array2d< real64 > >( AcousticWaveEquationSEM::viewKeyStruct::receiverCoordinatesString() ).toView(); + + for( int i = 0; i < 3; i++ ) + { + real64 tmp_double; + tmp_double=rCoord[0][i]; + rCoord[0][i]=sCoord[0][i]; + sCoord[0][i]=tmp_double; + } + + sCoord.registerTouch( hostMemorySpace ); + rCoord.registerTouch( hostMemorySpace ); + + std::cout << "sCoord :" << sCoord[0][0] <<" "<< sCoord[0][1] <<" "<< sCoord[0][2] << std::endl; + std::cout << "rCoord1 :" << rCoord[0][0] <<" "<< rCoord[0][1] <<" "<< rCoord[0][2] << std::endl; + std::cout << "rCoord2 :" << rCoord[1][0] <<" "<< rCoord[1][1] <<" "<< rCoord[1][2] << std::endl; + std::cout << "rCoord3 :" << rCoord[2][0] <<" "<< rCoord[2][1] <<" "<< rCoord[2][2] << std::endl; + std::cout << "rCoord4 :" << rCoord[3][0] <<" "<< rCoord[3][1] <<" "<< rCoord[3][2] << std::endl; + + //change timeSourceFrequency + std::cout << "timeSourceFrequency forward:" << *ptrTimeSourceFrequency << std::endl; + real32 newTimeFreq=2; + *ptrTimeSourceFrequency = newTimeFreq; + std::cout << "timeSourceFrequency backward:" << *ptrTimeSourceFrequency << std::endl; + + //reinit m_indexSeismoTrace + localIndex * ptrISeismo = &propagator->getReference< localIndex >( AcousticWaveEquationSEM::viewKeyStruct::indexSeismoTraceString() ); + *ptrISeismo = pReceivers.size( 0 )-1; + //reinit m_forward + localIndex * ptrForward = &propagator->getReference< localIndex >( AcousticWaveEquationSEM::viewKeyStruct::forwardString() ); + *ptrForward = 0; + + //"propagator->reinit()" not enough because state field not reinit to zero + //propagator->reinit(); + state.getProblemManager().applyInitialConditions(); + + array2d< real32 > rhsBackward; + rhsBackward.resize( 51, 1 ); + + arrayView2d< localIndex > const sNodeIds_new2 = propagator->getReference< array2d< localIndex > >( AcousticWaveEquationSEM::viewKeyStruct::sourceNodeIdsString() ).toView(); + sNodeIds_new2.move( hostMemorySpace, false ); + std::cout << "sNodeIds[0][0] second get2:" << sNodeIds_new2[0][0] << std::endl; + ASSERT_TRUE( sNodeIds_new2[0][0] == sNodesIdsAfterModif ); + + /*---------------------------------------------------*/ + // run backward solver + for( int i = 50; i > 0; i-- ) + { + rhsBackward[i][0]=WaveSolverUtils::evaluateRicker( time_n, *ptrTimeSourceFrequency, *ptrTimeSourceDelay, *ptrRickerOrder ); + propagator->explicitStepBackward( time_n, dt, i, domain, false ); + time_n -= dt; + //check source node in backward loop + arrayView2d< localIndex > const sNodeIds_loop = propagator->getReference< array2d< localIndex > >( AcousticWaveEquationSEM::viewKeyStruct::sourceNodeIdsString() ).toView(); + sNodeIds_loop.move( hostMemorySpace, false ); + ASSERT_TRUE( sNodeIds_loop[0][0] == sNodesIdsAfterModif ); + } + + // move it to CPU, if needed + pReceivers.move( hostMemorySpace, false ); + + localIndex mForward2 = propagator->getReference< localIndex >( AcousticWaveEquationSEM::viewKeyStruct::forwardString() ); + std::cout << "m_forward second get:" << mForward2 << std::endl; + ASSERT_TRUE( mForward2 == 0 ); + + arrayView2d< localIndex > const sNodeIds_new3 = propagator->getReference< array2d< localIndex > >( AcousticWaveEquationSEM::viewKeyStruct::sourceNodeIdsString() ).toView(); + sNodeIds_new3.move( hostMemorySpace, false ); + std::cout << "sNodeIds[0][0] get3:" << sNodeIds_new3[0][0] << std::endl; + ASSERT_TRUE( sNodeIds_new3[0][0] == sNodesIdsAfterModif ); + + real32 const timeSourceFrequency_new = propagator->getReference< real32 >( AcousticWaveEquationSEM::viewKeyStruct::timeSourceFrequencyString() ); + ASSERT_TRUE( std::abs( timeSourceFrequency_new - newTimeFreq ) < 1.e-8 ); + + /*std::cout << "pReceiver size(0):" << pReceivers.size(0) << std::endl; + std::cout << "pReceiver size(1):" << pReceivers.size(1) << std::endl;*/ + + + /*----------Save receiver backward----------------------*/ + array2d< real32 > qBackward; + qBackward.resize( 51, 1 ); + + real32 sum_ufb=0.; + real32 sum_qff=0.; + real32 sum_u2=0.; + real32 sum_q2=0.; + real32 sum_ff2=0.; + real32 sum_fb2=0.; + + // fill backward field at receiver. + for( int i=50; i > 0; i-- ) + { + /*std::cout << "back time: " << i*dt << std::endl; + std::cout << "back pReceivers1 " << i << ":" << pReceivers[i][0] << std::endl; + std::cout << "back pReceivers2 " << i << ":" << pReceivers[i][1] << std::endl; + std::cout << "back pReceivers3 " << i << ":" << pReceivers[i][2] << std::endl; + std::cout << "back pReceivers4 " << i << ":" << pReceivers[i][3] << std::endl; + std::cout << "back rhsBackward " << i << ":" << rhsBackward[i][0] << std::endl;*/ + qBackward[i][0] = pReceivers[i][0]; + } + + //check transitivity with sum + for( int i=0; i<51; i++ ) + { + sum_ufb += uForward[i][0]*rhsBackward[i][0]; + sum_qff += qBackward[i][0]*rhsForward[i][0]; + + sum_u2 += uForward[i][0]*uForward[i][0]; + sum_q2 += qBackward[i][0]*qBackward[i][0]; + sum_ff2 += rhsForward[i][0]*rhsForward[i][0]; + sum_fb2 += rhsBackward[i][0]*rhsBackward[i][0]; + /*std::cout << "sum evol sum_ufb:" << sum_ufb << " / sum_qff:" << sum_qff << std::endl; + std::cout << "uForward:" << uForward[i][0] << " / qBackward:" << qBackward[i][0] << std::endl; + std::cout << "ufb:" << uForward[i][0]*rhsBackward[i][0] << " / qff:" << qBackward[i][0]*rhsForward[i][0] << std::endl;*/ + } + + // check scalar products and are non null + ASSERT_TRUE( sum_ufb > 1.e-8 ); + ASSERT_TRUE( sum_qff > 1.e-8 ); + + // check || - ||/max(||f||.||q||,||f'||.||u||) < 10^1or2 x epsilon_machine with f rhs direct and f' rhs backward + std::cout << ": " << sum_ufb << " / : " << sum_qff << std::endl; + std::cout << "|| - ||=" << std::abs( sum_ufb-sum_qff ) << " / ||f||.||q||=" << std::sqrt( sum_q2*sum_ff2 ); + std::cout << " / ||f'||.||u||=" << std::sqrt( sum_fb2*sum_u2 ) << " / ||f||.||f'||=" << std::sqrt( sum_ff2*sum_fb2 ) << std::endl; + real32 diffToCheck; + diffToCheck=std::abs( sum_ufb-sum_qff ) / std::max( std::sqrt( sum_fb2*sum_u2 ), std::sqrt( sum_q2*sum_ff2 )); + std::cout << " Diff to compare with 2.e-4: " << diffToCheck << std::endl; + ASSERT_TRUE( diffToCheck < 2.e-4 ); +} + +int main( int argc, char * * argv ) +{ + ::testing::InitGoogleTest( &argc, argv ); + g_commandLineOptions = *geos::basicSetup( argc, argv ); + int const result = RUN_ALL_TESTS(); + geos::basicCleanup(); + return result; +}