diff --git a/.github/workflows/ci_tests.yml b/.github/workflows/ci_tests.yml index 0a327822c9b..0606ecf9a96 100644 --- a/.github/workflows/ci_tests.yml +++ b/.github/workflows/ci_tests.yml @@ -7,7 +7,7 @@ concurrency: cancel-in-progress: true env: - GEOSX_TPL_TAG: 238-63 + GEOSX_TPL_TAG: 228-73 jobs: @@ -23,9 +23,12 @@ jobs: steps: - name: Check that the PR is not a draft (cancel rest of jobs otherwise) run: | - echo "Is PR a draft?" - echo ${{ toJSON(github.event.pull_request.draft) }} - if [[ ${{ toJSON(github.event.pull_request.draft) }} == true ]]; then "false" ; else "true"; fi + draft_status=$(curl -H "Accept: application/vnd.github+json" \ + https://api.github.com/repos/${{ github.repository }}/pulls/${{ github.event.number }} \ + | jq '.draft') + echo "Is PR a draft? (Check PR's draft status )" + echo $draft_status + if [[ $draft_status == false ]]; then "true" ; else "false"; fi # PR must be assigned to be merged. # This job will fail if this is not the case. @@ -43,6 +46,23 @@ jobs: echo $id if [[ $id != null ]]; then "true" ; else "false"; fi + # PR must be labeled "ready to be merged" to run CUDA jobs. + # This job will fail if this is not the case. + # Note: Label should only be added if PR is in merge queue. + check_ready_to_be_merged: + needs: [check_pull_request_is_not_a_draft] + runs-on: ubuntu-22.04 + steps: + - name: Check that the PR is ready to be merged + run: | + labels=$(curl -H "Accept: application/vnd.github+json" \ + https://api.github.com/repos/${{ github.repository }}/pulls/${{ github.event.number }} \ + | jq '.labels') + ready_label="ready to be merged" + echo "Is PR Ready to be Merged? (Check for $ready_label label)" + echo "Labels found are: $labels" + if [[ $labels == *"$ready_label"* ]]; then "true" ; else "false"; fi + check_submodules: needs: [check_pull_request_is_not_a_draft] runs-on: ubuntu-22.04 @@ -164,8 +184,9 @@ jobs: # CUDA jobs should only be run if PR is ready to merge. cuda_builds_merge_ready: name: ${{ matrix.name }} - needs: [check_pull_request_is_not_a_draft] - if: contains( github.event.pull_request.labels.*.name, format('flag{0} ready to be merged', ':')) + needs: + - check_pull_request_is_not_a_draft + - check_ready_to_be_merged strategy: # In-progress jobs will not be cancelled if there is a failure @@ -226,6 +247,7 @@ jobs: needs: - check_pull_request_is_not_a_draft - check_pull_request_is_assigned + - check_ready_to_be_merged - check_submodules - code_style - documentation diff --git a/.readthedocs.yml b/.readthedocs.yml index 760835abafc..cac04680e9d 100644 --- a/.readthedocs.yml +++ b/.readthedocs.yml @@ -8,9 +8,10 @@ build: apt_packages: - npm - plantuml + - texlive tools: nodejs: "16" - python: "3.7" + python: "3.10" jobs: post_install: - npm install -g npm@9.8.1 diff --git a/host-configs/TOTAL/pangea3-gcc8.4.1-openmpi-4.1.2.cmake b/host-configs/TOTAL/pangea3-gcc8.4.1-openmpi-4.1.2.cmake index 5f3c90ce11f..ee75f9079b7 100644 --- a/host-configs/TOTAL/pangea3-gcc8.4.1-openmpi-4.1.2.cmake +++ b/host-configs/TOTAL/pangea3-gcc8.4.1-openmpi-4.1.2.cmake @@ -108,4 +108,7 @@ set(ENABLE_PETSC OFF CACHE BOOL "") set(ENABLE_HYPRE ON CACHE BOOL "") set(ENABLE_HYPRE_DEVICE "CUDA" CACHE BOOL "") +# activate workaround for fmt formatter +set(ENABLE_FMT_CONST_FORMATTER_WORKAROUND ON CACHE BOOL "") + include( ${CMAKE_CURRENT_LIST_DIR}/../tpls.cmake ) diff --git a/inputFiles/multiphaseFlowFractures/deadoil_3ph_corey_2d_edfm.xml b/inputFiles/multiphaseFlowFractures/deadoil_3ph_corey_2d_edfm.xml index 6b46a7ceb90..410fa9f286c 100644 --- a/inputFiles/multiphaseFlowFractures/deadoil_3ph_corey_2d_edfm.xml +++ b/inputFiles/multiphaseFlowFractures/deadoil_3ph_corey_2d_edfm.xml @@ -28,7 +28,7 @@ + + + diff --git a/inputFiles/multiphaseFlowFractures/deadoil_3ph_corey_2d_impermeableFault.xml b/inputFiles/multiphaseFlowFractures/deadoil_3ph_corey_2d_impermeableFault.xml index 2118a7beb18..3b037407799 100644 --- a/inputFiles/multiphaseFlowFractures/deadoil_3ph_corey_2d_impermeableFault.xml +++ b/inputFiles/multiphaseFlowFractures/deadoil_3ph_corey_2d_impermeableFault.xml @@ -28,7 +28,7 @@ + + + diff --git a/inputFiles/multiphaseFlowFractures/deadoil_3ph_corey_edfm_1d.xml b/inputFiles/multiphaseFlowFractures/deadoil_3ph_corey_edfm_1d.xml index 6f1addf105c..6d979155181 100644 --- a/inputFiles/multiphaseFlowFractures/deadoil_3ph_corey_edfm_1d.xml +++ b/inputFiles/multiphaseFlowFractures/deadoil_3ph_corey_edfm_1d.xml @@ -10,7 +10,7 @@ useMass="0" targetRegions="{ Domain, Fracture }"> @@ -23,7 +23,7 @@ @@ -99,6 +99,11 @@ + + + diff --git a/inputFiles/multiphaseFlowFractures/deadoil_3ph_corey_pedfm_impermeableFault_smoke.xml b/inputFiles/multiphaseFlowFractures/deadoil_3ph_corey_pedfm_impermeableFault_smoke.xml index 2118a7beb18..3b037407799 100644 --- a/inputFiles/multiphaseFlowFractures/deadoil_3ph_corey_pedfm_impermeableFault_smoke.xml +++ b/inputFiles/multiphaseFlowFractures/deadoil_3ph_corey_pedfm_impermeableFault_smoke.xml @@ -28,7 +28,7 @@ + + + diff --git a/inputFiles/wavePropagation/acous3D_vti_smoke.xml b/inputFiles/wavePropagation/acous3D_vti_smoke.xml new file mode 100644 index 00000000000..c3da486c9b3 --- /dev/null +++ b/inputFiles/wavePropagation/acous3D_vti_smoke.xml @@ -0,0 +1,236 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/integratedTests b/integratedTests index 7f2d88633c7..f682bc89094 160000 --- a/integratedTests +++ b/integratedTests @@ -1 +1 @@ -Subproject commit 7f2d88633c7f35528a49257f89bb11ae475c6e4c +Subproject commit f682bc89094a79ca4d0f059378f69c3ad0a7f82b diff --git a/scripts/SchemaToRSTDocumentation.py b/scripts/SchemaToRSTDocumentation.py index 2dcd03d6f95..0710515381c 100644 --- a/scripts/SchemaToRSTDocumentation.py +++ b/scripts/SchemaToRSTDocumentation.py @@ -258,13 +258,14 @@ def buildTableValues(type_map, link_string='XML', include_defaults=True): # Build documentation tables with open('%s.rst' % (complete_output), 'w') as output_handle: # Write the file header - output_handle.write('======================\n') + output_handle.write('######################\n') output_handle.write('Datastructure Index\n') - output_handle.write('======================\n\n') + output_handle.write('######################\n\n') # Parse the input schema definitions + output_handle.write('**************************\n\n') output_handle.write('Input Schema Definitions\n') - output_handle.write('========================\n\n') + output_handle.write('**************************\n\n') output_handle.write(':download:`XML Schema <%s/../schema.xsd>`\n\n' % (sphinx_path)) diff --git a/src/cmake/GeosxConfig.cmake b/src/cmake/GeosxConfig.cmake index c6760667f37..124b5da4b8d 100644 --- a/src/cmake/GeosxConfig.cmake +++ b/src/cmake/GeosxConfig.cmake @@ -4,6 +4,7 @@ set( PREPROCESSOR_DEFINES ARRAY_BOUNDS_CHECK CUDA CUDA_NVTOOLSEXT HIP + FMT_CONST_FORMATTER_WORKAROUND FORTRAN_MANGLE_NO_UNDERSCORE FPE HYPRE diff --git a/src/cmake/blt b/src/cmake/blt index 5a792c1775e..a7f0a6ecc4f 160000 --- a/src/cmake/blt +++ b/src/cmake/blt @@ -1 +1 @@ -Subproject commit 5a792c1775e7a7628d84dcde31652a689f1df7b5 +Subproject commit a7f0a6ecc4fdfa1724399b1454c3909b9ee02e81 diff --git a/src/conf.py b/src/conf.py index 5ec6a685023..cfcb471952f 100644 --- a/src/conf.py +++ b/src/conf.py @@ -99,7 +99,7 @@ # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom # ones. extensions = [ - 'sphinx.ext.imgmath', + 'sphinx_design', 'sphinx.ext.todo', 'sphinx.ext.autodoc', 'sphinx.ext.doctest', @@ -149,31 +149,19 @@ pygments_style = 'sphinx' -# -- Options for HTML output ------------------------------------------------- - -try: - if read_the_docs_build: - import sphinx_rtd_theme - else: - import pydata_sphinx_theme -except: - html_theme = 'classic' - html_theme_options = { - 'codebgcolor': 'lightgrey', - 'stickysidebar': 'true' - } - html_theme_path = [] -else: - html_theme = "pydata_sphinx_theme" - html_theme_options = {'navigation_depth': -1, 'collapse_navigation': False} - html_theme_path = [ - "_themes", - ] - - if read_the_docs_build: - html_theme = 'sphinx_rtd_theme' - html_theme_options = {} - html_theme_path = [sphinx_rtd_theme.get_html_theme_path()] +# -- Theme options ---------------------------------------------- +extensions += [ + 'sphinx_rtd_theme', +] + +html_theme = "sphinx_rtd_theme" +# html_theme = "pydata_sphinx_theme" + +html_theme_options = { + 'navigation_depth': -1, + 'collapse_navigation': False +} + # The name for this set of Sphinx documents. If None, it defaults to # " v documentation". @@ -194,58 +182,15 @@ # Add any paths that contain custom static files (such as style sheets) here, # relative to this directory. They are copied after the builtin static files, # so a file named "default.css" will overwrite the builtin "default.css". -html_static_path = ['./docs/sphinx/_static'] - -html_context = { - 'css_files': [ - '_static/theme_overrides.css', # override wide tables in RTD theme - ], - } - -# If not '', a 'Last updated on:' timestamp is inserted at every page bottom, -# using the given strftime format. -#html_last_updated_fmt = '%b %d, %Y' - -# If true, SmartyPants will be used to convert quotes and dashes to -# typographically correct entities. -#html_use_smartypants = True - -# Custom sidebar templates, maps document names to template names. -#html_sidebars = {} - -# Additional templates that should be rendered to pages, maps page names to -# template names. -#html_additional_pages = {} - -# If false, no module index is generated. -#html_domain_indices = True -# If false, no index is generated. -#html_use_index = True - -# If true, the index is split into individual pages for each letter. -#html_split_index = False - -# If true, links to the reST sources are added to the pages. -#html_show_sourcelink = True - -# If true, "Created using Sphinx" is shown in the HTML footer. Default is True. -#html_show_sphinx = True - -# If true, "(C) Copyright ..." is shown in the HTML footer. Default is True. -#html_show_copyright = True - -# If true, an OpenSearch description file will be output, and all pages will -# contain a tag referring to it. The value of this option must be the -# base URL from which the finished HTML is served. -#html_use_opensearch = '' +html_static_path = ['./docs/sphinx/_static'] -# This is the file name suffix for HTML files (e.g. ".xhtml"). -#html_file_suffix = None +html_css_files = [ + 'theme_overrides.css', +] # -- Options for HTMLHelp output --------------------------------------------- - # Output file base name for HTML help builder. htmlhelp_basename = 'GEOSXdoc' @@ -309,21 +254,17 @@ # Additional stuff for the LaTeX preamble. latex_elements['preamble'] = '\\usepackage{amsmath}\n\\usepackage{amssymb}\n\\usepackage[retainorgcmds]{IEEEtrantools}\n' -imgmath_image_format='svg' -imgmath_font_size=14 + ##################################################### # add LaTeX macros f = open('docs/sphinx/latex_macros.sty') - -try: - imgmath_latex_preamble # check whether this is already defined -except NameError: - imgmath_latex_preamble = "" +imgmath_latex_preamble = "" for macro in f: # used when building latex and pdf versions latex_elements['preamble'] += macro + '\n' # used when building html version imgmath_latex_preamble += macro + '\n' + ##################################################### diff --git a/src/coreComponents/common/CMakeLists.txt b/src/coreComponents/common/CMakeLists.txt index 823195f0ab0..47e7076a7be 100644 --- a/src/coreComponents/common/CMakeLists.txt +++ b/src/coreComponents/common/CMakeLists.txt @@ -15,6 +15,7 @@ set( common_headers Path.hpp Span.hpp Stopwatch.hpp + Timer.hpp Tensor.hpp TimingMacros.hpp TypeDispatch.hpp diff --git a/src/coreComponents/common/GeosxConfig.hpp.in b/src/coreComponents/common/GeosxConfig.hpp.in index 6e366aa3020..2e0bad1a7f9 100644 --- a/src/coreComponents/common/GeosxConfig.hpp.in +++ b/src/coreComponents/common/GeosxConfig.hpp.in @@ -41,6 +41,9 @@ /// Enables use of HIP (CMake option ENABLE_HIP) #cmakedefine GEOS_USE_HIP +/// Workaround for FMT compilation issue on some NVCC/PowerPC machines (CMake option ENABLE_FMT_CONST_FORMATTER_WORKAROUND) +#cmakedefine GEOS_USE_FMT_CONST_FORMATTER_WORKAROUND + /// Enables use of PVTPackage (CMake option ENABLE_PVTPackage) #cmakedefine GEOSX_USE_PVTPackage diff --git a/src/coreComponents/common/Timer.hpp b/src/coreComponents/common/Timer.hpp new file mode 100644 index 00000000000..909340b7e5e --- /dev/null +++ b/src/coreComponents/common/Timer.hpp @@ -0,0 +1,57 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 TotalEnergies + * Copyright (c) 2019- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file Timer.hpp + */ + +#ifndef GEOS_COMMON_TIMER_HPP +#define GEOS_COMMON_TIMER_HPP + +#include + +namespace geos +{ + +/** + * @class Timer + * @brief Object that times the duration of its existence. + */ +class Timer +{ +public: + + /** + * @brief Constructor. The time the object is alive is added to @p duration. + * @param duration A reference to the duration to add to. + */ + Timer( std::chrono::system_clock::duration & duration ): + m_start( std::chrono::system_clock::now() ), + m_duration( duration ) + {} + + /// Destructor. Adds to the referenced duration. + ~Timer() + { m_duration += std::chrono::system_clock::now() - m_start; } + +private: + /// The time at which this object was constructed. + std::chrono::system_clock::time_point const m_start; + /// A reference to the duration to add to. + std::chrono::system_clock::duration & m_duration; +}; + +} + +#endif // GEOS_COMMON_TIMER_HPP diff --git a/src/coreComponents/constitutive/contact/ContactBase.cpp b/src/coreComponents/constitutive/contact/ContactBase.cpp index 01d9c5738c8..e8006e27f8f 100644 --- a/src/coreComponents/constitutive/contact/ContactBase.cpp +++ b/src/coreComponents/constitutive/contact/ContactBase.cpp @@ -104,6 +104,12 @@ void ContactBase::allocateConstitutiveData( Group & parent, // this check is necessary to ensure that the coordinates are strictly increasing if( apertureTransition > apertureValues[apertureValues.size()-1] ) { + GEOS_LOG( GEOS_FMT ( "Adding aperture transition for table {}:", m_apertureTableName ) ); + std::ostringstream s_orig; + for( localIndex i = 0; i < apertureValues.size(); i++ ) + s_orig << "[ " << apertureValues[i] << ", " << hydraulicApertureValues[i] << " ] "; + GEOS_LOG( GEOS_FMT ( " Original table = {}", s_orig.str())); + coords.emplaceBack( 0, apertureTransition ); hydraulicApertureValues.emplace_back( apertureTransition ); // if the aperture transition is larger than 0, we keep enlarging the table @@ -114,6 +120,11 @@ void ContactBase::allocateConstitutiveData( Group & parent, hydraulicApertureValues.emplace_back( apertureTransition*10e9 ); apertureTable.reInitializeFunction(); } + + std::ostringstream s_mod; + for( localIndex i = 0; i < apertureValues.size(); i++ ) + s_mod << "[ " << apertureValues[i] << ", " << hydraulicApertureValues[i] << " ] "; + GEOS_LOG( GEOS_FMT ( " Modified table = {}", s_mod.str())); } m_apertureTable = &apertureTable; diff --git a/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/BrineEnthalpy.cpp b/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/BrineEnthalpy.cpp index 3059fdfaef4..ecb661cf1f7 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/BrineEnthalpy.cpp +++ b/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/BrineEnthalpy.cpp @@ -37,7 +37,7 @@ namespace PVTProps namespace { -real64 MichaelidesBrineEnthalpy( real64 const & T, +real64 michaelidesBrineEnthalpy( real64 const & T, real64 const & m ) { @@ -51,26 +51,21 @@ real64 MichaelidesBrineEnthalpy( real64 const & T, real64 x1, x2, h1, h2, dh; x1 = 1000.0 / (1000.0 + 58.44 * m); - x2 = 58.44 * m / (1000.0 + 58.44 * m); - - dh = 0.0; for( localIndex i = 0; i < 4; ++i ) + { for( localIndex j = 0; j < 3; ++j ) { - dh += a[i][j] * pow( T, real64( i )) * pow( m, real64( j )); - } + } dh *= 4.184 / (1000.0 + 58.44 * m); - h1 = 0.12453e-4 * pow( T, 3.0 ) - 0.45137e-2 * pow( T, 2.0 ) + 4.81155 * T - 29.578; - h2 = (-0.83624e-3 * pow( T, 3.0 ) + 0.16792 * pow( T, 2.0 ) - 25.9293 * T) * 4.184 / 58.44; return ( (x1 * h1 + x2 * h2 + m * dh) * 1000.0 ); @@ -86,7 +81,7 @@ void calculateBrineEnthalpy( PTTableCoordinates const & tableCoords, for( localIndex i = 0; i < nTemperatures; ++i ) { real64 const TC = tableCoords.getTemperature( i ); - enthalpies[i] = MichaelidesBrineEnthalpy( TC, m ); + enthalpies[i] = michaelidesBrineEnthalpy( TC, m ); } } diff --git a/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/CO2Enthalpy.cpp b/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/CO2Enthalpy.cpp index dbe77c24909..b71f44fadbf 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/CO2Enthalpy.cpp +++ b/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/CO2Enthalpy.cpp @@ -35,7 +35,7 @@ namespace PVTProps namespace { -real64 HelmholtzCO2Enthalpy( real64 const & T, +real64 helmholtzCO2Enthalpy( real64 const & T, real64 const & rho ) { static const real64 dc = 467.6; @@ -269,12 +269,16 @@ CO2Enthalpy::calculateCO2Enthalpy( PTTableCoordinates const & tableCoords, localIndex const nPressures = tableCoords.nPressures(); localIndex const nTemperatures = tableCoords.nTemperatures(); + // Note that the enthalpy values given in Span and Wagner assume a reference enthalphy defined as: h_0 = 0 J/kg at T_0 = 298.15 K + // Therefore, the enthalpy computed using the Span and Wagner methid must be shifted by the enthalpy at 298.15 K + real64 const referenceEnthalpy = 5.0584e5; // J/kg + for( localIndex i = 0; i < nPressures; ++i ) { for( localIndex j = 0; j < nTemperatures; ++j ) { real64 const TC = tableCoords.getTemperature( j ); - enthalpies[j*nPressures+i] = HelmholtzCO2Enthalpy( TC, densities[j*nPressures+i] ); + enthalpies[j*nPressures+i] = helmholtzCO2Enthalpy( TC, densities[j*nPressures+i] ) + referenceEnthalpy; } } } diff --git a/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/PhillipsBrineDensity.cpp b/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/PhillipsBrineDensity.cpp index 2f89d6fb734..b0c4c35eaf2 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/PhillipsBrineDensity.cpp +++ b/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/PhillipsBrineDensity.cpp @@ -17,7 +17,7 @@ */ #include "constitutive/fluid/multifluid/CO2Brine/functions/PhillipsBrineDensity.hpp" - +#include "constitutive/fluid/multifluid/CO2Brine/functions/PureWaterProperties.hpp" #include "functions/FunctionManager.hpp" namespace geos @@ -70,6 +70,44 @@ void calculateBrineDensity( PTTableCoordinates const & tableCoords, } } +void calculatePureWaterDensity( PTTableCoordinates const & tableCoords, + string const & functionName, + array1d< real64 > const & densities ) +{ + // if no salinity, we fall back to the standard approach in three steps + // 1- Get the saturation density as a function of temperature + // 2- Get the saturation pressure as a function of temperature + // 3- Get the pure water density + + TableFunction const * waterSatDensityTable = + PureWaterProperties::makeSaturationDensityTable( functionName, FunctionManager::getInstance() ); + TableFunction const * waterSatPressureTable = + PureWaterProperties::makeSaturationPressureTable( functionName, FunctionManager::getInstance() ); + + localIndex const nPressures = tableCoords.nPressures(); + localIndex const nTemperatures = tableCoords.nTemperatures(); + + for( localIndex i = 0; i < nPressures; ++i ) + { + real64 const P = tableCoords.getPressure( i ) / 1e5; + + for( localIndex j = 0; j < nTemperatures; ++j ) + { + real64 const T = tableCoords.getTemperature( j ); + + // Step 1: get the saturation density + real64 const waterSatDensity = waterSatDensityTable->evaluate( &T ); + // Step 2: get the saturation pressure + real64 const waterSatPressure = waterSatPressureTable->evaluate( &T ); + // Step 3: get the pure water density + // Note: for now, we keep a constant water compressibility for consistency with the Ezrokhi model + // In the future, we should query the water compressibility as a function of pressure and temperature in a table + real64 const waterCompressibility = 4.5e-10; // Pa-1 // TODO: consolidate to a unique file as Dick started doing + densities[j*nPressures+i] = waterSatDensity * exp( waterCompressibility * ( P - waterSatPressure ) ); + } + } +} + TableFunction const * makeDensityTable( string_array const & inputParams, string const & functionName, FunctionManager & functionManager ) @@ -93,7 +131,21 @@ TableFunction const * makeDensityTable( string_array const & inputParams, } array1d< real64 > densities( tableCoords.nPressures() * tableCoords.nTemperatures() ); - calculateBrineDensity( tableCoords, salinity, densities ); + if( !isZero( salinity ) ) + { + // if we are in the range of validity of the Phillips method, everything is good + // if we are not, we issue a warning message + calculateBrineDensity( tableCoords, salinity, densities ); + GEOS_LOG_RANK_0_IF( salinity < 0.25, + GEOS_FMT( "{}: Warning! The salinity value of {} is below the range of validity of the Phillips model, results may be inaccurate", + functionName, salinity ) ); + } + else + { + // the Phillips correlation is inaccurate in the absence of salinity. + // since this is a very frequent case, we implement an alternate (more accurate) method below + calculatePureWaterDensity( tableCoords, functionName, densities ); + } string const tableName = functionName + "_table"; if( functionManager.hasGroup< TableFunction >( tableName ) ) diff --git a/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/PureWaterProperties.cpp b/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/PureWaterProperties.cpp index 906b0e6860d..6305d9588c1 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/PureWaterProperties.cpp +++ b/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/PureWaterProperties.cpp @@ -152,7 +152,7 @@ PureWaterProperties::makeSaturationDensityTable( string const & functionName, densities[2] = 998.21; densities[3] = 997.05; densities[4] = 995.65; - densities[5] = 995.65; + densities[5] = 992.25; densities[6] = 988.04; densities[7] = 983.2; densities[8] = 977.76; diff --git a/src/coreComponents/dataRepository/DataContext.hpp b/src/coreComponents/dataRepository/DataContext.hpp index 6b30face341..f04de33375c 100644 --- a/src/coreComponents/dataRepository/DataContext.hpp +++ b/src/coreComponents/dataRepository/DataContext.hpp @@ -223,4 +223,26 @@ struct GEOS_FMT_NS::formatter< geos::dataRepository::DataContext > : GEOS_FMT_NS } }; +// The following workaround is needed to fix compilation with NVCC on some PowerPC machines. +// The issue causes the following assertion error message: +// "Cannot format an argument. To make type T formattable provide a formatter specialization" +// The standard definition of the has_const_formatter check of fmt fails due to a compiler bug, see the issue below: +// https://github.com/fmtlib/fmt/issues/2746 +// The workaround was originally implemented in fmt: +// https://github.com/fmtlib/fmt/commit/70de324aa801eaf52e94c402d526a3849797c620 +// but later removed: +// https://github.com/fmtlib/fmt/commit/466e0650ec2d153d255a40ec230eb77d7f1c3334 +// This workaround provides a specialization of the const formatter check for the DataContext object. +// The formatter is defined within this file, and therefore the check is not needed. +// The scope of the check override is as small as possible to solve the current issue. +#ifdef GEOS_USE_FMT_CONST_FORMATTER_WORKAROUND +template<> +constexpr auto GEOS_FMT_NS::detail::has_const_formatter< geos::dataRepository::DataContext, GEOS_FMT_NS::format_context >() -> bool +{ + return true; +} +#endif +// End of the workaround for fmt compilation issues + + #endif /* GEOS_DATAREPOSITORY_DATACONTEXT_HPP_ */ diff --git a/src/coreComponents/events/EventManager.cpp b/src/coreComponents/events/EventManager.cpp index 510a3e25dca..d06353fae8d 100644 --- a/src/coreComponents/events/EventManager.cpp +++ b/src/coreComponents/events/EventManager.cpp @@ -179,9 +179,8 @@ bool EventManager::run( DomainPartition & domain ) subEvent->checkEvents( m_time, m_dt, m_cycle, domain ); // Print debug information for logLevel >= 1 - GEOS_LOG_LEVEL_RANK_0( 1, - " Event: " << m_currentSubEvent << " (" << subEvent->getName() << "), dt_request=" << subEvent->getCurrentEventDtRequest() << ", forecast=" << - subEvent->getForecast() ); + GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( "Event: {} ({}), dt_request={}, forecast={}", + m_currentSubEvent, subEvent->getName(), subEvent->getCurrentEventDtRequest(), subEvent->getForecast() ) ); // Execute, signal events bool earlyReturn = false; diff --git a/src/coreComponents/fileIO/vtk/VTKPolyDataWriterInterface.cpp b/src/coreComponents/fileIO/vtk/VTKPolyDataWriterInterface.cpp index b762d23a51f..d2f1b2ac3b7 100644 --- a/src/coreComponents/fileIO/vtk/VTKPolyDataWriterInterface.cpp +++ b/src/coreComponents/fileIO/vtk/VTKPolyDataWriterInterface.cpp @@ -58,7 +58,7 @@ VTKPolyDataWriterInterface::VTKPolyDataWriterInterface( string name ): {} static int -toVTKCellType( ElementType const elementType ) +toVTKCellType( ElementType const elementType, localIndex const numNodes ) { switch( elementType ) { @@ -70,7 +70,16 @@ toVTKCellType( ElementType const elementType ) case ElementType::Tetrahedron: return VTK_TETRA; case ElementType::Pyramid: return VTK_PYRAMID; case ElementType::Wedge: return VTK_WEDGE; - case ElementType::Hexahedron: return VTK_HEXAHEDRON; + case ElementType::Hexahedron: + switch( numNodes ) + { + case 8: + return VTK_HEXAHEDRON; + case 27: + return VTK_QUADRATIC_HEXAHEDRON; + default: + return VTK_HEXAHEDRON; + } case ElementType::Prism5: return VTK_PENTAGONAL_PRISM; case ElementType::Prism6: return VTK_HEXAGONAL_PRISM; case ElementType::Prism7: return VTK_POLYHEDRON; @@ -101,7 +110,7 @@ toVTKCellType( ElementType const elementType ) * nodes for mapping purpose while keeping a face streams data structure. The negative values are * converted to positives when generating the VTK_POLYHEDRON. Check getVtkCells() for more details. */ -static std::vector< int > getVtkConnectivity( ElementType const elementType ) +static std::vector< int > getVtkConnectivity( ElementType const elementType, localIndex const numNodes ) { switch( elementType ) { @@ -113,7 +122,21 @@ static std::vector< int > getVtkConnectivity( ElementType const elementType ) case ElementType::Tetrahedron: return { 0, 1, 2, 3 }; case ElementType::Pyramid: return { 0, 1, 3, 2, 4 }; case ElementType::Wedge: return { 0, 4, 2, 1, 5, 3 }; - case ElementType::Hexahedron: return { 0, 1, 3, 2, 4, 5, 7, 6 }; + case ElementType::Hexahedron: + switch( numNodes ) + { + case 8: + return { 0, 1, 3, 2, 4, 5, 7, 6 }; + break; + case 27: + // Numbering convention changed between VTK writer 1.0 and 2.2, see + // https://discourse.julialang.org/t/writevtk-node-numbering-for-27-node-lagrange-hexahedron/93698/8 + // GEOS uses 1.0 API + return { 0, 2, 8, 6, 18, 20, 26, 24, 1, 5, 7, 3, 19, 23, 25, 21, 9, 11, 17, 15, 12, 14, 10, 16, 4, 22, 17 }; + break; + default: + return { }; // TODO + } case ElementType::Prism5: return { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 }; case ElementType::Prism6: return { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 }; case ElementType::Prism7: return {-9, @@ -289,6 +312,7 @@ getSurface( FaceElementSubRegion const & subRegion, for( localIndex ei = 0; ei < subRegion.size(); ei++ ) { auto const & nodes = elemToNodes[ei]; + auto const numNodes = nodes.size(); ElementType const elementType = subRegion.getElementType( ei ); std::vector< int > vtkOrdering; @@ -299,7 +323,7 @@ getSurface( FaceElementSubRegion const & subRegion, } else { - vtkOrdering = getVtkConnectivity( elementType ); + vtkOrdering = getVtkConnectivity( elementType, numNodes ); } connectivity.clear(); @@ -317,7 +341,7 @@ getSurface( FaceElementSubRegion const & subRegion, } cellArray->InsertNextCell( vtkOrdering.size(), connectivity.data() ); - cellTypes.emplace_back( toVTKCellType( elementType ) ); + cellTypes.emplace_back( toVTKCellType( elementType, numNodes ) ); } auto points = vtkSmartPointer< vtkPoints >::New(); @@ -436,7 +460,7 @@ getVtkCells( CellElementRegion const & region, localIndex numConn = 0; region.forElementSubRegions< CellElementSubRegion >( [&]( CellElementSubRegion const & subRegion ) { - numConn += subRegion.size() * getVtkConnectivity( subRegion.getElementType() ).size(); + numConn += subRegion.size() * getVtkConnectivity( subRegion.getElementType(), subRegion.nodeList().size( 1 ) ).size(); } ); return numConn; }(); @@ -455,18 +479,19 @@ getVtkCells( CellElementRegion const & region, localIndex connOffset = 0; region.forElementSubRegions< CellElementSubRegion >( [&]( CellElementSubRegion const & subRegion ) { - cellTypes.insert( cellTypes.end(), subRegion.size(), toVTKCellType( subRegion.getElementType() ) ); - std::vector< int > const vtkOrdering = getVtkConnectivity( subRegion.getElementType() ); - localIndex const numVtkData = vtkOrdering.size(); auto const nodeList = subRegion.nodeList().toViewConst(); + auto subRegionNumNodes = nodeList.size( 1 ); + cellTypes.insert( cellTypes.end(), subRegion.size(), toVTKCellType( subRegion.getElementType(), subRegionNumNodes ) ); + std::vector< int > const vtkOrdering = getVtkConnectivity( subRegion.getElementType(), subRegionNumNodes ); + localIndex const numVtkData = vtkOrdering.size(); -// For all geosx element, the corresponding VTK data are copied in "connectivity". -// Local nodes are mapped to global indices. Any negative value in "vtkOrdering" -// corresponds to the number of faces or the number of nodes per faces, and they -// are copied as positive values. -// Here we privilege code simplicity. This can be more efficient (less tests) if the code is -// specialized for each type of subregion. -// This is not a time sensitive part of the code. Can be optimized later if needed. + // For all geosx element, the corresponding VTK data are copied in "connectivity". + // Local nodes are mapped to global indices. Any negative value in "vtkOrdering" + // corresponds to the number of faces or the number of nodes per faces, and they + // are copied as positive values. + // Here we privilege code simplicity. This can be more efficient (less tests) if the code is + // specialized for each type of subregion. + // This is not a time sensitive part of the code. Can be optimized later if needed. forAll< parallelHostPolicy >( subRegion.size(), [&]( localIndex const c ) { localIndex const elemConnOffset = connOffset + c * numVtkData; diff --git a/src/coreComponents/finiteElement/elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp b/src/coreComponents/finiteElement/elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp index 6d613207f6b..c6fddf9ad6a 100644 --- a/src/coreComponents/finiteElement/elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp +++ b/src/coreComponents/finiteElement/elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp @@ -415,7 +415,7 @@ class Qk_Hexahedron_Lagrange_GaussLobatto final : public FiniteElementBase real64 ( &B )[6] ); /** - * @brief computes the non-zero contributions of the d.o.f. indexd by q to the + * @brief computes the non-zero contributions of the d.o.f. indexed by q to the * stiffness matrix R, i.e., the superposition matrix of first derivatives * of the shape functions. * @param q The quadrature point index @@ -429,6 +429,92 @@ class Qk_Hexahedron_Lagrange_GaussLobatto final : public FiniteElementBase real64 const (&X)[numNodes][3], FUNC && func ); + /** + * @brief computes the matrix B in the case of quasi-stiffness (e.g. for pseudo-acoustic case), defined as J^{-T}A_xy J^{-1}/det(J), where + * J is the Jacobian matrix, and A_xy is a zero matrix except on A_xy(1,1) = 1 and A_xy(2,2) = 1. + * @param qa The 1d quadrature point index in xi0 direction (0,1) + * @param qb The 1d quadrature point index in xi1 direction (0,1) + * @param qc The 1d quadrature point index in xi2 direction (0,1) + * @param X Array containing the coordinates of the support points. + * @param J Array to store the Jacobian + * @param B Array to store the matrix B, in Voigt notation + */ + GEOS_HOST_DEVICE + static void + computeBxyMatrix( int const qa, + int const qb, + int const qc, + real64 const (&X)[numNodes][3], + real64 ( &J )[3][3], + real64 ( &B )[6] ); + + /** + * @brief computes the non-zero contributions of the d.o.f. indexed by q to the + * partial-stiffness matrix R, i.e., the superposition matrix of first derivatives in x and y + * of the shape functions. Warning, the matrix B is obtained by computeBxyMatrix instead of usual one. + * @param q The quadrature point index + * @param X Array containing the coordinates of the support points. + * @param func Callback function accepting three parameters: i, j and R_ij + */ + template< typename FUNC > + GEOS_HOST_DEVICE + static void + computeStiffnessxyTerm( int q, + real64 const (&X)[numNodes][3], + FUNC && func ); + + /** + * @brief computes the matrix B in the case of quasi-stiffness (e.g. for pseudo-acoustic case), defined as J^{-T}A_z J^{-1}/det(J), where + * J is the Jacobian matrix, and A_z is a zero matrix except on A_z(3,3) = 1. + * @param qa The 1d quadrature point index in xi0 direction (0,1) + * @param qb The 1d quadrature point index in xi1 direction (0,1) + * @param qc The 1d quadrature point index in xi2 direction (0,1) + * @param X Array containing the coordinates of the support points. + * @param J Array to store the Jacobian + * @param B Array to store the matrix B, in Voigt notation + */ + GEOS_HOST_DEVICE + static void + computeBzMatrix( int const qa, + int const qb, + int const qc, + real64 const (&X)[numNodes][3], + real64 ( &J )[3][3], + real64 ( &B )[6] ); + + /** + * @brief computes the non-zero contributions of the d.o.f. indexed by q to the + * partial-stiffness matrix R, i.e., the superposition matrix of first derivatives in z only + * of the shape functions. Warning, the matrix B is obtained by computeBzMatrix instead of usual one. + * @param q The quadrature point index + * @param X Array containing the coordinates of the support points. + * @param func Callback function accepting three parameters: i, j and R_ij + */ + template< typename FUNC > + GEOS_HOST_DEVICE + static void + computeStiffnesszTerm( int q, + real64 const (&X)[numNodes][3], + FUNC && func ); + +/** + * @brief Computes the "Grad(Phi)*B*Grad(Phi)" coefficient of the stiffness term. The matrix B must be provided and Phi denotes a basis + * function. + * @param qa The 1d quadrature point index in xi0 direction (0,1) + * @param qb The 1d quadrature point index in xi1 direction (0,1) + * @param qc The 1d quadrature point index in xi2 direction (0,1) + * @param B Array of the B matrix, in Voigt notation + * @param func Callback function accepting three parameters: i, j and R_ij + */ + template< typename FUNC > + GEOS_HOST_DEVICE + static void + computeGradPhiBGradPhi( int qa, + int qb, + int qc, + real64 const (&B)[6], + FUNC && func ); + /** * @brief computes the non-zero contributions of the d.o.f. indexd by q to the * x-part of the first order stiffness matrix R, i.e., the matrix composed of the @@ -885,16 +971,6 @@ computeDampingTerm( int q, return sqrt( LvArray::math::abs( LvArray::tensorOps::symDeterminant< 2 >( B ) ) )*GL_BASIS::weight( qa )*GL_BASIS::weight( qb ); } -/** - * @brief computes the matrix B, defined as J^{-T}J^{-1}/det(J), where J is the Jacobian matrix, - * at the given Gauss-Lobatto point. - * @param qa The 1d quadrature point index in xi0 direction (0,1) - * @param qb The 1d quadrature point index in xi1 direction (0,1) - * @param qc The 1d quadrature point index in xi2 direction (0,1) - * @param X Array containing the coordinates of the support points. - * @param J Array to store the Jacobian - * @param B Array to store the matrix B, in Voigt notation - */ template< typename GL_BASIS > GEOS_HOST_DEVICE inline @@ -922,7 +998,183 @@ computeBMatrix( int const qa, LvArray::tensorOps::symInvert< 3 >( B ); } +template< typename GL_BASIS > +GEOS_HOST_DEVICE +GEOS_FORCE_INLINE +void +Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >:: +computeBzMatrix( int const qa, + int const qb, + int const qc, + real64 const (&X)[numNodes][3], + real64 (& J)[3][3], + real64 (& B)[6] ) +{ + jacobianTransformation( qa, qb, qc, X, J ); + real64 const detJ = LvArray::tensorOps::determinant< 3 >( J ); + + real64 Jinv[3][3] = {{0}}; + LvArray::tensorOps::invert< 3 >( Jinv, J ); + + // compute det(J)*J^{-1}Az*J^{-T}, using Voigt notation for B + B[0] = detJ*(Jinv[0][2]*Jinv[0][2]); + B[1] = detJ*(Jinv[1][2]*Jinv[1][2]); + B[2] = detJ*(Jinv[2][2]*Jinv[2][2]); + B[3] = detJ*(Jinv[1][2]*Jinv[2][2]); + B[4] = detJ*(Jinv[0][2]*Jinv[2][2]); + B[5] = detJ*(Jinv[0][2]*Jinv[1][2]); +} + +template< typename GL_BASIS > +GEOS_HOST_DEVICE +GEOS_FORCE_INLINE +void +Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >:: +computeBxyMatrix( int const qa, + int const qb, + int const qc, + real64 const (&X)[numNodes][3], + real64 (& J)[3][3], + real64 (& B)[6] ) +{ + jacobianTransformation( qa, qb, qc, X, J ); + real64 const detJ = LvArray::tensorOps::determinant< 3 >( J ); + + real64 Jinv[3][3] = {{0}}; + LvArray::tensorOps::invert< 3 >( Jinv, J ); + + // compute det(J)*J^{-1}Axy*J^{-T}, using Voigt notation for B + B[0] = detJ*(Jinv[0][0]*Jinv[0][0] + Jinv[0][1]*Jinv[0][1]); + B[1] = detJ*(Jinv[1][1]*Jinv[1][1] + Jinv[1][0]*Jinv[1][0]); + B[2] = detJ*(Jinv[2][0]*Jinv[2][0] + Jinv[2][1]*Jinv[2][1]); + B[3] = detJ*(Jinv[1][0]*Jinv[2][0] + Jinv[1][1]*Jinv[2][1]); + B[4] = detJ*(Jinv[0][0]*Jinv[2][0] + Jinv[0][1]*Jinv[2][1]); + B[5] = detJ*(Jinv[0][0]*Jinv[1][0] + Jinv[0][1]*Jinv[1][1]); +} + +template< typename GL_BASIS > +template< typename FUNC > +GEOS_HOST_DEVICE +GEOS_FORCE_INLINE +void +Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >:: +computeGradPhiBGradPhi( int qa, + int qb, + int qc, + real64 const (&B)[6], + FUNC && func ) +{ + // diagonal terms + for( int i=0; i +template< typename FUNC > +GEOS_HOST_DEVICE +GEOS_FORCE_INLINE +void +Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >:: +computeStiffnessxyTerm( int q, + real64 const (&X)[numNodes][3], + FUNC && func ) +{ + real64 B[6] = {0}; + real64 J[3][3] = {{0}}; + int qa, qb, qc; + GL_BASIS::TensorProduct3D::multiIndex( q, qa, qb, qc ); + computeBxyMatrix( qa, qb, qc, X, J, B ); // The only change! + computeGradPhiBGradPhi( qa, qb, qc, B, func ); +} +template< typename GL_BASIS > +template< typename FUNC > +GEOS_HOST_DEVICE +GEOS_FORCE_INLINE +void +Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >:: +computeStiffnesszTerm( int q, + real64 const (&X)[numNodes][3], + FUNC && func ) +{ + real64 B[6] = {0}; + real64 J[3][3] = {{0}}; + int qa, qb, qc; + GL_BASIS::TensorProduct3D::multiIndex( q, qa, qb, qc ); + computeBzMatrix( qa, qb, qc, X, J, B ); // The only change! + computeGradPhiBGradPhi( qa, qb, qc, B, func ); +} template< typename GL_BASIS > template< typename FUNC > @@ -1104,8 +1356,6 @@ computeFirstOrderStiffnessTerm( int q, } } - - for( int i=0; i::solve( Vector const & rhs, if( m_params.logLevel >= 1 ) { - GEOS_LOG_RANK_0( "\t\tLinear Solver | " << m_result.status << + GEOS_LOG_RANK_0( " Linear Solver | " << m_result.status << " | Iterations: " << m_result.numIterations << " | Final Rel Res: " << m_result.residualReduction << " | Setup Time: " << m_result.setupTime << " s" << diff --git a/src/coreComponents/linearAlgebra/interfaces/direct/SuperLUDist.cpp b/src/coreComponents/linearAlgebra/interfaces/direct/SuperLUDist.cpp index 3e510657c59..52a62dfa2b3 100644 --- a/src/coreComponents/linearAlgebra/interfaces/direct/SuperLUDist.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/direct/SuperLUDist.cpp @@ -298,7 +298,7 @@ void SuperLUDist< LAI >::solve( Vector const & rhs, if( m_params.logLevel >= 1 ) { - GEOS_LOG_RANK_0( "\t\tLinear Solver | " << m_result.status << + GEOS_LOG_RANK_0( " Linear Solver | " << m_result.status << " | Iterations: " << m_result.numIterations << " | Final Rel Res: " << m_result.residualReduction << " | Setup Time: " << m_result.setupTime << " s" << diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMGR.cpp b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMGR.cpp index a595f7f6c07..f82594e87ef 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMGR.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMGR.cpp @@ -65,11 +65,11 @@ void hypre::mgr::createMGR( LinearSolverParameters const & params, if( params.logLevel >= 1 ) { - GEOS_LOG_RANK_0( numComponentsPerField ); + GEOS_LOG_RANK_0( GEOS_FMT( " MGR preconditioner: numComponentsPerField = {}", numComponentsPerField ) ); } if( params.logLevel >= 4 ) { - GEOS_LOG_RANK_VAR( mgrData.pointMarkers ); + GEOS_LOG_RANK( GEOS_FMT( " MGR preconditioner: pointMarkers = {}", mgrData.pointMarkers ) ); } switch( params.mgr.strategy ) diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreSolver.cpp b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreSolver.cpp index 0780fb9d01c..04516475be2 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreSolver.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreSolver.cpp @@ -275,7 +275,7 @@ void HypreSolver::solve( HypreVector const & rhs, if( m_params.logLevel >= 1 ) { - GEOS_LOG_RANK_0( "\t\tLinear Solver | " << m_result.status << + GEOS_LOG_RANK_0( " Linear Solver | " << m_result.status << " | Iterations: " << m_result.numIterations << " | Final Rel Res: " << m_result.residualReduction << " | Make Restrictor Time: " << m_makeRestrictorTime << diff --git a/src/coreComponents/linearAlgebra/interfaces/petsc/PetscSolver.cpp b/src/coreComponents/linearAlgebra/interfaces/petsc/PetscSolver.cpp index b877ae3ce17..19e5b1a4eae 100644 --- a/src/coreComponents/linearAlgebra/interfaces/petsc/PetscSolver.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/petsc/PetscSolver.cpp @@ -148,7 +148,7 @@ void PetscSolver::solve( PetscVector const & rhs, if( m_params.logLevel >= 1 ) { - GEOS_LOG_RANK_0( "\t\tLinear Solver | " << m_result.status << + GEOS_LOG_RANK_0( " Linear Solver | " << m_result.status << " | Iterations: " << m_result.numIterations << " | Final Rel Res: " << m_result.residualReduction << " | Setup Time: " << m_result.setupTime << " s" << diff --git a/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosSolver.cpp b/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosSolver.cpp index bf38e8b901e..01b97d058c3 100644 --- a/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosSolver.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosSolver.cpp @@ -156,7 +156,7 @@ void TrilinosSolver::solve( EpetraVector const & rhs, if( m_params.logLevel >= 1 ) { - GEOS_LOG_RANK_0( "\t\tLinear Solver | " << m_result.status << + GEOS_LOG_RANK_0( " Linear Solver | " << m_result.status << " | Iterations: " << m_result.numIterations << " | Final Rel Res: " << m_result.residualReduction << " | Setup Time: " << m_result.setupTime << " s" << diff --git a/src/coreComponents/mainInterface/GeosxState.cpp b/src/coreComponents/mainInterface/GeosxState.cpp index 1a2c29928ba..e965076b89c 100644 --- a/src/coreComponents/mainInterface/GeosxState.cpp +++ b/src/coreComponents/mainInterface/GeosxState.cpp @@ -18,6 +18,7 @@ #include "mainInterface/ProblemManager.hpp" #include "mainInterface/initialization.hpp" #include "mesh/mpiCommunications/CommunicationTools.hpp" +#include "common/Timer.hpp" // TPL includes #include @@ -43,35 +44,6 @@ GeosxState & getGlobalState() return *currentGlobalState; } - -/** - * @class Timer - * @brief Object that times the duration of its existence. - */ -class Timer -{ -public: - - /** - * @brief Constructor. The time the object is alive is added to @p duration. - * @param duration A reference to the duration to add to. - */ - Timer( std::chrono::system_clock::duration & duration ): - m_start( std::chrono::system_clock::now() ), - m_duration( duration ) - {} - - /// Destructor. Adds to the referenced duration. - ~Timer() - { m_duration += std::chrono::system_clock::now() - m_start; } - -private: - /// The time at which this object was constructed. - std::chrono::system_clock::time_point const m_start; - /// A reference to the duration to add to. - std::chrono::system_clock::duration & m_duration; -}; - ////////////////////////////////////////////////////////////////////////////////////////////////////////////////// string durationToString( std::chrono::system_clock::duration const duration ) { diff --git a/src/coreComponents/mainInterface/GeosxState.hpp b/src/coreComponents/mainInterface/GeosxState.hpp index e32a484f27d..8f3bb7402f5 100644 --- a/src/coreComponents/mainInterface/GeosxState.hpp +++ b/src/coreComponents/mainInterface/GeosxState.hpp @@ -24,7 +24,6 @@ // System includes #include -#include #include // Forward declaration of conduit::Node. diff --git a/src/coreComponents/mesh/ElementRegionBase.cpp b/src/coreComponents/mesh/ElementRegionBase.cpp index 707b5470a6a..349e06234f3 100644 --- a/src/coreComponents/mesh/ElementRegionBase.cpp +++ b/src/coreComponents/mesh/ElementRegionBase.cpp @@ -54,6 +54,9 @@ string ElementRegionBase::verifyMeshBodyName( Group const & meshBodies, string meshBodyName = meshBodyBlockName; localIndex const numberOfMeshBodies = meshBodies.numSubGroups(); + GEOS_THROW_IF( numberOfMeshBodies == 0, + "No MeshBodies found in this problem, please check if correct input file is provided", InputError ); + if( numberOfMeshBodies == 1 ) { string const & onlyMeshBodyName = meshBodies.getGroup( 0 ).getName(); diff --git a/src/coreComponents/mesh/generators/InternalMeshGenerator.cpp b/src/coreComponents/mesh/generators/InternalMeshGenerator.cpp index 1850338ffe3..f5b1a04e9c0 100644 --- a/src/coreComponents/mesh/generators/InternalMeshGenerator.cpp +++ b/src/coreComponents/mesh/generators/InternalMeshGenerator.cpp @@ -534,6 +534,7 @@ static void getElemToNodesRelationInBox( ElementType const elementType, void InternalMeshGenerator::fillCellBlockManager( CellBlockManager & cellBlockManager, array1d< int > const & partition ) { GEOS_MARK_FUNCTION; + m_spatialPartition.setPartitions( partition[0], partition[1], partition[2] ); // Partition based on even spacing to get load balance // Partition geometrical boundaries will be corrected in the end. @@ -583,6 +584,7 @@ void InternalMeshGenerator::fillCellBlockManager( CellBlockManager & cellBlockMa { m_numElemsTotal[i] += m_nElems[i][block]; } + GEOS_ERROR_IF( partition[i] > m_numElemsTotal[i], "Number of partitions in a direction should not exceed the number of elements in that direction" ); elemCenterCoords[i].resize( m_numElemsTotal[i] ); array1d< real64 > elemCenterCoordsLocal( m_numElemsTotal[i] ); diff --git a/src/coreComponents/mesh/generators/VTKUtilities.cpp b/src/coreComponents/mesh/generators/VTKUtilities.cpp index 4546860d608..7cbc37f0a83 100644 --- a/src/coreComponents/mesh/generators/VTKUtilities.cpp +++ b/src/coreComponents/mesh/generators/VTKUtilities.cpp @@ -802,13 +802,12 @@ ensureNoEmptyRank( vtkSmartPointer< vtkDataSet > mesh, // First we sort the donor in order of the number of elems they contain std::stable_sort( donorRanks.begin(), donorRanks.end(), [&elemCounts] ( auto i1, auto i2 ) - { return elemCounts[i1] < elemCounts[i2]; } ); + { return elemCounts[i1] > elemCounts[i2]; } ); // Then, if my position is "i" in the donorRanks array, I will send half of my elems to the i-th recipient integer const myRank = MpiWrapper::commRank(); - auto const myPosition = - LvArray::sortedArrayManipulation::find( donorRanks.begin(), donorRanks.size(), myRank ); - bool const isDonor = myPosition != donorRanks.size(); + auto const pos = std::find( donorRanks.begin(), donorRanks.end(), myRank ); + bool const isDonor = ( pos != donorRanks.end() ); // step 3: my rank was selected to donate cells, let's proceed // we need to make a distinction between two configurations @@ -820,6 +819,7 @@ ensureNoEmptyRank( vtkSmartPointer< vtkDataSet > mesh, // we use a strategy that preserves load balancing if( isDonor && donorRanks.size() >= recipientRanks.size() ) { + auto const myPosition = std::distance( donorRanks.begin(), pos ); if( myPosition < recipientRanks.size() ) { integer const recipientRank = recipientRanks[myPosition]; @@ -834,10 +834,11 @@ ensureNoEmptyRank( vtkSmartPointer< vtkDataSet > mesh, // we just want the simulation to run and count on ParMetis/PTScotch to restore load balancing else if( isDonor && donorRanks.size() < recipientRanks.size() ) { + auto const myPosition = std::distance( donorRanks.begin(), pos ); localIndex firstRecipientPosition = 0; for( integer iRank = 0; iRank < myPosition; ++iRank ) { - firstRecipientPosition += elemCounts[iRank] - 1; + firstRecipientPosition += elemCounts[donorRanks[iRank]] - 1; } if( firstRecipientPosition < recipientRanks.size() ) { diff --git a/src/coreComponents/mesh/mpiCommunications/NeighborCommunicator.cpp b/src/coreComponents/mesh/mpiCommunications/NeighborCommunicator.cpp index 7672a1b7ea3..62d8ac5de66 100644 --- a/src/coreComponents/mesh/mpiCommunications/NeighborCommunicator.cpp +++ b/src/coreComponents/mesh/mpiCommunications/NeighborCommunicator.cpp @@ -501,7 +501,7 @@ int NeighborCommunicator::packCommSizeForSync( FieldIdentifiers const & fieldsTo for( auto const & iter : fieldsToBeSync.getFields() ) { - FieldLocation location; + FieldLocation location{}; fieldsToBeSync.getLocation( iter.first, location ); switch( location ) { @@ -559,7 +559,7 @@ void NeighborCommunicator::packCommBufferForSync( FieldIdentifiers const & field for( auto const & iter : fieldsToBeSync.getFields() ) { - FieldLocation location; + FieldLocation location{}; fieldsToBeSync.getLocation( iter.first, location ); switch( location ) { diff --git a/src/coreComponents/physicsSolvers/CMakeLists.txt b/src/coreComponents/physicsSolvers/CMakeLists.txt index e834480cd08..3468462b8b6 100644 --- a/src/coreComponents/physicsSolvers/CMakeLists.txt +++ b/src/coreComponents/physicsSolvers/CMakeLists.txt @@ -135,6 +135,8 @@ set( physicsSolvers_headers wavePropagation/ElasticFirstOrderWaveEquationSEMKernel.hpp wavePropagation/AcousticFirstOrderWaveEquationSEM.hpp wavePropagation/AcousticFirstOrderWaveEquationSEMKernel.hpp + wavePropagation/AcousticVTIWaveEquationSEM.hpp + wavePropagation/AcousticVTIWaveEquationSEMKernel.hpp ) # @@ -198,7 +200,7 @@ set( physicsSolvers_sources wavePropagation/ElasticWaveEquationSEM.cpp wavePropagation/ElasticFirstOrderWaveEquationSEM.cpp wavePropagation/AcousticFirstOrderWaveEquationSEM.cpp - + wavePropagation/AcousticVTIWaveEquationSEM.cpp ) include( solidMechanics/kernels/SolidMechanicsKernels.cmake) diff --git a/src/coreComponents/physicsSolvers/SolverBase.cpp b/src/coreComponents/physicsSolvers/SolverBase.cpp index d3908f86b33..1d3add01579 100644 --- a/src/coreComponents/physicsSolvers/SolverBase.cpp +++ b/src/coreComponents/physicsSolvers/SolverBase.cpp @@ -20,6 +20,7 @@ #include "linearAlgebra/solvers/KrylovSolver.hpp" #include "mesh/DomainPartition.hpp" #include "math/interpolation/Interpolation.hpp" +#include "common/Timer.hpp" #if defined(GEOSX_USE_PYGEOSX) #include "python/PySolverType.hpp" @@ -97,7 +98,6 @@ SolverBase::SolverBase( string const & name, SolverBase::~SolverBase() = default; - void SolverBase::initialize_postMeshGeneration() { ExecutableGroup::initialize_postMeshGeneration(); @@ -362,11 +362,13 @@ real64 SolverBase::linearImplicitStep( real64 const & time_n, // TODO: Nonlinear step does not call its own setup, need to decide on consistent behavior implicitStepSetup( time_n, dt, domain ); - // zero out matrix/rhs before assembly - m_localMatrix.zero(); - m_rhs.zero(); - { + Timer timer( m_timers["assemble"] ); + + // zero out matrix/rhs before assembly + m_localMatrix.zero(); + m_rhs.zero(); + arrayView1d< real64 > const localRhs = m_rhs.open(); // call assemble to fill the matrix and the rhs @@ -386,30 +388,38 @@ real64 SolverBase::linearImplicitStep( real64 const & time_n, localRhs ); m_rhs.close(); - } - if( m_assemblyCallback ) - { - // Make a copy of LA objects and ship off to the callback - array1d< real64 > localRhsCopy( m_rhs.localSize() ); - localRhsCopy.setValues< parallelDevicePolicy<> >( m_rhs.values() ); - m_assemblyCallback( m_localMatrix, std::move( localRhsCopy ) ); + if( m_assemblyCallback ) + { + // Make a copy of LA objects and ship off to the callback + array1d< real64 > localRhsCopy( m_rhs.localSize() ); + localRhsCopy.setValues< parallelDevicePolicy<> >( m_rhs.values() ); + m_assemblyCallback( m_localMatrix, std::move( localRhsCopy ) ); + } } - // TODO: Trilinos currently requires this, re-evaluate after moving to Tpetra-based solvers - if( m_precond ) { - m_precond->clear(); - } + Timer timer( m_timers["linear solver total"] ); - // Compose parallel LA matrix out of local matrix - m_matrix.create( m_localMatrix.toViewConst(), m_dofManager.numLocalDofs(), MPI_COMM_GEOSX ); + // TODO: Trilinos currently requires this, re-evaluate after moving to Tpetra-based solvers + if( m_precond ) + { + m_precond->clear(); + } - // Output the linear system matrix/rhs for debugging purposes - debugOutputSystem( 0.0, 0, 0, m_matrix, m_rhs ); + { + Timer timer_create( m_timers["linear solver create"] ); - // Solve the linear system - solveLinearSystem( m_dofManager, m_matrix, m_rhs, m_solution ); + // Compose parallel LA matrix out of local matrix + m_matrix.create( m_localMatrix.toViewConst(), m_dofManager.numLocalDofs(), MPI_COMM_GEOSX ); + } + + // Output the linear system matrix/rhs for debugging purposes + debugOutputSystem( 0.0, 0, 0, m_matrix, m_rhs ); + + // Solve the linear system + solveLinearSystem( m_dofManager, m_matrix, m_rhs, m_solution ); + } // Increment the solver statistics for reporting purposes m_solverStatistics.logNonlinearIteration( m_linearSolverResult.numIterations ); @@ -417,11 +427,19 @@ real64 SolverBase::linearImplicitStep( real64 const & time_n, // Output the linear system solution for debugging purposes debugOutputSolution( 0.0, 0, 0, m_solution ); - // apply the system solution to the fields/variables - applySystemSolution( m_dofManager, m_solution.values(), 1.0, dt, domain ); + { + Timer timer( m_timers["apply solution"] ); + + // apply the system solution to the fields/variables + applySystemSolution( m_dofManager, m_solution.values(), 1.0, dt, domain ); + } + + { + Timer timer( m_timers["update state"] ); - // update non-primary variables (constitutive models) - updateState( domain ); + // update non-primary variables (constitutive models) + updateState( domain ); + } // final step for completion of timestep. typically secondary variable updates and cleanup. implicitStepComplete( time_n, dt, domain ); @@ -459,27 +477,37 @@ bool SolverBase::lineSearch( real64 const & time_n, // main loop for the line search. for( integer lineSearchIteration = 0; lineSearchIteration < maxNumberLineSearchCuts; ++lineSearchIteration ) { - // cut the scale factor by half. This means that the scale factors will - // have values of -0.5, -0.25, -0.125, ... - localScaleFactor *= lineSearchCutFactor; - cumulativeScale += localScaleFactor; - - if( !checkSystemSolution( domain, dofManager, solution.values(), localScaleFactor ) ) { - GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( " Line search {}, solution check failed", lineSearchIteration ) ); - continue; - } + Timer timer( m_timers["apply solution"] ); - applySystemSolution( dofManager, solution.values(), localScaleFactor, dt, domain ); + // cut the scale factor by half. This means that the scale factors will + // have values of -0.5, -0.25, -0.125, ... + localScaleFactor *= lineSearchCutFactor; + cumulativeScale += localScaleFactor; - // update non-primary variables (constitutive models) - updateState( domain ); + if( !checkSystemSolution( domain, dofManager, solution.values(), localScaleFactor ) ) + { + GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( " Line search {}, solution check failed", lineSearchIteration ) ); + continue; + } - // re-assemble system - localMatrix.zero(); - rhs.zero(); + applySystemSolution( dofManager, solution.values(), localScaleFactor, dt, domain ); + } { + Timer timer( m_timers["update state"] ); + + // update non-primary variables (constitutive models) + updateState( domain ); + } + + { + Timer timer( m_timers["assemble"] ); + + // re-assemble system + localMatrix.zero(); + rhs.zero(); + arrayView1d< real64 > const localRhs = rhs.open(); assembleSystem( time_n, dt, domain, dofManager, localMatrix, localRhs ); applyBoundaryConditions( time_n, dt, domain, dofManager, localMatrix, localRhs ); @@ -491,12 +519,12 @@ bool SolverBase::lineSearch( real64 const & time_n, std::cout << GEOS_FMT( " Line search @ {:0.3f}: ", cumulativeScale ); } - // get residual norm - residualNorm = calculateResidualNorm( time_n, dt, domain, dofManager, rhs.values() ); - - if( getLogLevel() >= 1 && logger::internal::rank==0 ) { - std::cout << std::endl; + Timer timer( m_timers["convergence check"] ); + + // get residual norm + residualNorm = calculateResidualNorm( time_n, dt, domain, dofManager, rhs.values() ); + GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( " ( R ) = ( {:4.2e} )", residualNorm ) ); } // if the residual norm is less than the last residual, we can proceed to the @@ -545,41 +573,52 @@ bool SolverBase::lineSearchWithParabolicInterpolation( real64 const & time_n, while( residualNormT >= (1.0 - alpha*localScaleFactor)*residualNorm0 ) { - real64 const previousLocalScaleFactor = localScaleFactor; - // Apply the three point parabolic model - if( lineSearchIteration == 0 ) - { - localScaleFactor *= sigma1; - } - else { - localScaleFactor = interpolation::parabolicInterpolationThreePoints( lamc, lamm, ff0, ffT, ffm ); - } + Timer timer( m_timers["apply solution"] ); - // Update x; keep the books on lambda - real64 const deltaLocalScaleFactor = ( localScaleFactor - previousLocalScaleFactor ); - cumulativeScale += deltaLocalScaleFactor; + real64 const previousLocalScaleFactor = localScaleFactor; + // Apply the three point parabolic model + if( lineSearchIteration == 0 ) + { + localScaleFactor *= sigma1; + } + else + { + localScaleFactor = interpolation::parabolicInterpolationThreePoints( lamc, lamm, ff0, ffT, ffm ); + } - if( !checkSystemSolution( domain, dofManager, solution.values(), deltaLocalScaleFactor ) ) - { - GEOS_LOG_LEVEL_RANK_0( 1, " Line search " << lineSearchIteration << ", solution check failed" ); - continue; + // Update x; keep the books on lambda + real64 const deltaLocalScaleFactor = ( localScaleFactor - previousLocalScaleFactor ); + cumulativeScale += deltaLocalScaleFactor; + + if( !checkSystemSolution( domain, dofManager, solution.values(), deltaLocalScaleFactor ) ) + { + GEOS_LOG_LEVEL_RANK_0( 1, " Line search " << lineSearchIteration << ", solution check failed" ); + continue; + } + + applySystemSolution( dofManager, solution.values(), deltaLocalScaleFactor, dt, domain ); } - applySystemSolution( dofManager, solution.values(), deltaLocalScaleFactor, dt, domain ); + { + Timer timer( m_timers["update state"] ); - updateState( domain ); + updateState( domain ); + } lamm = lamc; lamc = localScaleFactor; // Keep the books on the function norms - // re-assemble system - // TODO: add a flag to avoid a completely useless Jacobian computation: rhs is enough - localMatrix.zero(); - rhs.zero(); { + Timer timer( m_timers["assemble"] ); + + // re-assemble system + // TODO: add a flag to avoid a completely useless Jacobian computation: rhs is enough + localMatrix.zero(); + rhs.zero(); + arrayView1d< real64 > const localRhs = rhs.open(); assembleSystem( time_n, dt, domain, dofManager, localMatrix, localRhs ); applyBoundaryConditions( time_n, dt, domain, dofManager, localMatrix, localRhs ); @@ -591,8 +630,13 @@ bool SolverBase::lineSearchWithParabolicInterpolation( real64 const & time_n, std::cout << GEOS_FMT( " Line search @ {:0.3f}: ", cumulativeScale ); } - // get residual norm - residualNormT = calculateResidualNorm( time_n, dt, domain, dofManager, rhs.values() ); + { + Timer timer( m_timers["convergence check"] ); + + // get residual norm + residualNormT = calculateResidualNorm( time_n, dt, domain, dofManager, rhs.values() ); + GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( " ( R ) = ( {:4.2e} )", residualNormT ) ); + } ffm = ffT; ffT = residualNormT*residualNormT; lineSearchIteration += 1; @@ -717,9 +761,7 @@ real64 SolverBase::nonlinearImplicitStep( real64 const & time_n, // increment the solver statistics for reporting purposes m_solverStatistics.logOuterLoopIteration(); - GEOS_LOG_LEVEL_RANK_0( 1, " " ); GEOS_LOG_LEVEL_RANK_0( 1, "---------- Configuration did not converge. Testing new configuration. ----------" ); - GEOS_LOG_LEVEL_RANK_0( 1, " " ); } } else if( !attemptedSimplestConfiguration ) @@ -795,11 +837,13 @@ bool SolverBase::solveNonlinearSystem( real64 const & time_n, { GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( " Attempt: {:2}, ConfigurationIter: {:2}, NewtonIter: {:2}", dtAttempt, configurationLoopIter, newtonIter ) ); - // zero out matrix/rhs before assembly - m_localMatrix.zero(); - m_rhs.zero(); - { + Timer timer( m_timers["assemble"] ); + + // zero out matrix/rhs before assembly + m_localMatrix.zero(); + m_rhs.zero(); + arrayView1d< real64 > const localRhs = m_rhs.open(); // call assemble to fill the matrix and the rhs @@ -820,23 +864,27 @@ bool SolverBase::solveNonlinearSystem( real64 const & time_n, m_rhs.close(); - } - if( m_assemblyCallback ) - { - // Make a copy of LA objects and ship off to the callback - array1d< real64 > localRhsCopy( m_rhs.localSize() ); - localRhsCopy.setValues< parallelDevicePolicy<> >( m_rhs.values() ); - m_assemblyCallback( m_localMatrix, std::move( localRhsCopy ) ); + if( m_assemblyCallback ) + { + // Make a copy of LA objects and ship off to the callback + array1d< real64 > localRhsCopy( m_rhs.localSize() ); + localRhsCopy.setValues< parallelDevicePolicy<> >( m_rhs.values() ); + m_assemblyCallback( m_localMatrix, std::move( localRhsCopy ) ); + } } - // get residual norm - real64 residualNorm = calculateResidualNorm( time_n, stepDt, domain, m_dofManager, m_rhs.values() ); + real64 residualNorm = 0; + { + Timer timer( m_timers["convergence check"] ); + // get residual norm + residualNorm = calculateResidualNorm( time_n, stepDt, domain, m_dofManager, m_rhs.values() ); + GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( " ( R ) = ( {:4.2e} )", residualNorm ) ); + } - GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( " ( R ) = ( {:4.2e} ) ; ", residualNorm ) ); if( newtonIter > 0 ) { - GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( " Last LinSolve(iter,res) = ( {:3}, {:4.2e} ) ; ", + GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( " Last LinSolve(iter,res) = ( {:3}, {:4.2e} )", m_linearSolverResult.numIterations, m_linearSolverResult.residualReduction ) ); } @@ -917,21 +965,29 @@ bool SolverBase::solveNonlinearSystem( real64 const & time_n, krylovParams.relTolerance = eisenstatWalker( residualNorm, lastResidual, krylovParams.weakestTol ); } - // TODO: Trilinos currently requires this, re-evaluate after moving to Tpetra-based solvers - if( m_precond ) { - m_precond->clear(); - } + Timer timer( m_timers["linear solver total"] ); - // Compose parallel LA matrix/rhs out of local LA matrix/rhs - // - m_matrix.create( m_localMatrix.toViewConst(), m_dofManager.numLocalDofs(), MPI_COMM_GEOSX ); + // TODO: Trilinos currently requires this, re-evaluate after moving to Tpetra-based solvers + if( m_precond ) + { + m_precond->clear(); + } - // Output the linear system matrix/rhs for debugging purposes - debugOutputSystem( time_n, cycleNumber, newtonIter, m_matrix, m_rhs ); + { + Timer timer_setup( m_timers["linear solver create"] ); - // Solve the linear system - solveLinearSystem( m_dofManager, m_matrix, m_rhs, m_solution ); + // Compose parallel LA matrix/rhs out of local LA matrix/rhs + // + m_matrix.create( m_localMatrix.toViewConst(), m_dofManager.numLocalDofs(), MPI_COMM_GEOSX ); + } + + // Output the linear system matrix/rhs for debugging purposes + debugOutputSystem( time_n, cycleNumber, newtonIter, m_matrix, m_rhs ); + + // Solve the linear system + solveLinearSystem( m_dofManager, m_matrix, m_rhs, m_solution ); + } // Increment the solver statistics for reporting purposes m_solverStatistics.logNonlinearIteration( m_linearSolverResult.numIterations ); @@ -939,26 +995,34 @@ bool SolverBase::solveNonlinearSystem( real64 const & time_n, // Output the linear system solution for debugging purposes debugOutputSolution( time_n, cycleNumber, newtonIter, m_solution ); - // Compute the scaling factor for the Newton update - scaleFactor = scalingForSystemSolution( domain, m_dofManager, m_solution.values() ); - - if( getLogLevel() >= 1 ) { - GEOS_LOG_RANK_0( getName() + ": Global solution scaling factor = " << scaleFactor ); - } + Timer timer( m_timers["apply solution"] ); - if( !checkSystemSolution( domain, m_dofManager, m_solution.values(), scaleFactor ) ) - { - // TODO try chopping (similar to line search) - GEOS_LOG_RANK_0( " Solution check failed. Newton loop terminated." ); - break; + // Compute the scaling factor for the Newton update + scaleFactor = scalingForSystemSolution( domain, m_dofManager, m_solution.values() ); + + if( getLogLevel() >= 1 ) + { + GEOS_LOG_RANK_0( GEOS_FMT( " {}: Global solution scaling factor = {}", getName(), scaleFactor ) ); + } + + if( !checkSystemSolution( domain, m_dofManager, m_solution.values(), scaleFactor ) ) + { + // TODO try chopping (similar to line search) + GEOS_LOG_RANK_0( " Solution check failed. Newton loop terminated." ); + break; + } + + // apply the system solution to the fields/variables + applySystemSolution( m_dofManager, m_solution.values(), scaleFactor, stepDt, domain ); } - // apply the system solution to the fields/variables - applySystemSolution( m_dofManager, m_solution.values(), scaleFactor, stepDt, domain ); + { + Timer timer( m_timers["update state"] ); - // update non-primary variables (constitutive models) - updateState( domain ); + // update non-primary variables (constitutive models) + updateState( domain ); + } lastResidual = residualNorm; } @@ -1145,15 +1209,27 @@ void SolverBase::solveLinearSystem( DofManager const & dofManager, if( params.solverType == LinearSolverParameters::SolverType::direct || !m_precond ) { std::unique_ptr< LinearSolverBase< LAInterface > > solver = LAInterface::createSolver( params ); - solver->setup( matrix ); - solver->solve( rhs, solution ); + { + Timer timer_setup( m_timers["linear solver setup"] ); + solver->setup( matrix ); + } + { + Timer timer_setup( m_timers["linear solver solve"] ); + solver->solve( rhs, solution ); + } m_linearSolverResult = solver->result(); } else { - m_precond->setup( matrix ); + { + Timer timer_setup( m_timers["linear solver setup"] ); + m_precond->setup( matrix ); + } std::unique_ptr< KrylovSolver< ParallelVector > > solver = KrylovSolver< ParallelVector >::create( params, matrix, *m_precond ); - solver->solve( rhs, solution ); + { + Timer timer_setup( m_timers["linear solver solve"] ); + solver->solve( rhs, solution ); + } m_linearSolverResult = solver->result(); } @@ -1236,6 +1312,20 @@ void SolverBase::cleanup( real64 const GEOS_UNUSED_PARAM( time_n ), DomainPartition & GEOS_UNUSED_PARAM( domain ) ) { m_solverStatistics.outputStatistics(); + + if( getLogLevel() > 0 ) + { + for( auto & timer : m_timers ) + { + real64 const time = std::chrono::duration< double >( timer.second ).count(); + real64 const minTime = MpiWrapper::min( time ); + real64 const maxTime = MpiWrapper::max( time ); + if( maxTime > 0 ) + { + GEOS_LOG_RANK_0( GEOS_FMT( "{}: {} time = {} s (min), {} s (max)", getName(), timer.first, minTime, maxTime ) ); + } + } + } } Timestamp SolverBase::getMeshModificationTimestamp( DomainPartition & domain ) const diff --git a/src/coreComponents/physicsSolvers/SolverBase.hpp b/src/coreComponents/physicsSolvers/SolverBase.hpp index 02f8d812062..35d47e44204 100644 --- a/src/coreComponents/physicsSolvers/SolverBase.hpp +++ b/src/coreComponents/physicsSolvers/SolverBase.hpp @@ -795,6 +795,7 @@ class SolverBase : public ExecutableGroup std::function< void( CRSMatrix< real64, globalIndex >, array1d< real64 > ) > m_assemblyCallback; + std::map< std::string, std::chrono::system_clock::duration > m_timers; private: /// List of names of regions the solver will be applied to diff --git a/src/coreComponents/physicsSolvers/contact/LagrangianContactSolver.cpp b/src/coreComponents/physicsSolvers/contact/LagrangianContactSolver.cpp index 057cc2f97aa..74842e67b9f 100644 --- a/src/coreComponents/physicsSolvers/contact/LagrangianContactSolver.cpp +++ b/src/coreComponents/physicsSolvers/contact/LagrangianContactSolver.cpp @@ -258,8 +258,6 @@ void LagrangianContactSolver::implicitStepComplete( real64 const & time_n, true ); } ); - - GEOS_LOG_LEVEL_RANK_0( 1, " ***** ImplicitStepComplete *****" ); } void LagrangianContactSolver::postProcessInput() @@ -708,10 +706,14 @@ real64 LagrangianContactSolver::calculateResidualNorm( real64 const & GEOS_UNUSE // Add 0 just to match Matlab code results globalResidualNorm[2] /= (m_initialResidual[2]+1.0); } - GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( " ( Rdisplacement, Rtraction, Rtotal ) = ( {:15.6e}, {:15.6e}, {:15.6e} );", - globalResidualNorm[0], - globalResidualNorm[1], - globalResidualNorm[2] ) ); + if( getLogLevel() >= 1 && logger::internal::rank == 0 ) + { + std::cout<< GEOS_FMT( + " ( Rdisplacement, Rtraction, Rtotal ) = ( {:15.6e}, {:15.6e}, {:15.6e} )", + globalResidualNorm[0], + globalResidualNorm[1], + globalResidualNorm[2] ); + } return globalResidualNorm[2]; } diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsEmbeddedFractures.cpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsEmbeddedFractures.cpp index de43fd2eada..54aac882187 100644 --- a/src/coreComponents/physicsSolvers/contact/SolidMechanicsEmbeddedFractures.cpp +++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsEmbeddedFractures.cpp @@ -632,7 +632,7 @@ real64 SolidMechanicsEmbeddedFractures::calculateResidualNorm( real64 const & ti if( getLogLevel() >= 1 && logger::internal::rank==0 ) { - std::cout << GEOS_FMT( "( RFracture ) = ( {:4.2e} ) ; ", fractureResidualNorm ); + std::cout << GEOS_FMT( " ( RFracture ) = ( {:4.2e} )", fractureResidualNorm ); } return sqrt( solidResidualNorm * solidResidualNorm + fractureResidualNorm * fractureResidualNorm ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp index af65bc95f56..1baef8810cb 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp @@ -334,7 +334,7 @@ real64 CompositionalMultiphaseFVM::calculateResidualNorm( real64 const & GEOS_UN if( getLogLevel() >= 1 && logger::internal::rank == 0 ) { - std::cout << GEOS_FMT( " ( R{} ) = ( {:4.2e} ) ; ( Renergy ) = ( {:4.2e} ) ; ", + std::cout << GEOS_FMT( " ( R{} ) = ( {:4.2e} ) ( Renergy ) = ( {:4.2e} )", coupledSolverAttributePrefix(), globalResidualNorm[0], globalResidualNorm[1] ); } } @@ -352,7 +352,7 @@ real64 CompositionalMultiphaseFVM::calculateResidualNorm( real64 const & GEOS_UN if( getLogLevel() >= 1 && logger::internal::rank == 0 ) { - std::cout << GEOS_FMT( " ( R{} ) = ( {:4.2e} ) ; ", coupledSolverAttributePrefix(), residualNorm ); + std::cout << GEOS_FMT( " ( R{} ) = ( {:4.2e} )", coupledSolverAttributePrefix(), residualNorm ); } } return residualNorm; @@ -431,18 +431,18 @@ real64 CompositionalMultiphaseFVM::scalingForSystemSolution( DomainPartition & d if( m_isThermal ) { - GEOS_LOG_LEVEL_RANK_0( 1, getName() + ": Max deltaPres = " << maxDeltaPres << ", Max deltaCompDens = " << maxDeltaCompDens << ", Max deltaTemp = " << maxDeltaTemp << " (before scaling)" ); + GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( " {}: Max deltaPres = {}, Max deltaCompDens = {}, Max deltaTemp = {} (before scaling)", getName(), maxDeltaPres, maxDeltaCompDens, maxDeltaTemp ) ); } else { - GEOS_LOG_LEVEL_RANK_0( 1, getName() + ": Max deltaPres = " << maxDeltaPres << ", Max deltaCompDens = " << maxDeltaCompDens << " (before scaling)" ); + GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( " {}: Max deltaPres = {}, Max deltaCompDens = {} (before scaling)", getName(), maxDeltaPres, maxDeltaCompDens ) ); } if( m_scalingType == ScalingType::Local ) { - GEOS_LOG_LEVEL_RANK_0( 1, getName() + ": Min pressure scaling factor = " << minPresScalingFactor ); - GEOS_LOG_LEVEL_RANK_0( 1, getName() + ": Min comp dens scaling factor = " << minCompDensScalingFactor ); + GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( " {}: Min pressure scaling factor = {}", getName(), minPresScalingFactor ) ); + GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( " {}: Min comp dens scaling factor = {}", getName(), minCompDensScalingFactor ) ); if( m_isThermal ) - GEOS_LOG_LEVEL_RANK_0( 1, getName() + ": Min temperature scaling factor = " << minTempScalingFactor ); + GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( " {}: Min temperature scaling factor = {}", getName(), minTempScalingFactor ) ); } return LvArray::math::max( scalingFactor, m_minScalingFactor ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp index 74b2809a3df..80667badced 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp @@ -706,7 +706,7 @@ real64 CompositionalMultiphaseHybridFVM::calculateResidualNorm( real64 const & G if( getLogLevel() >= 1 && logger::internal::rank == 0 ) { - std::cout << GEOS_FMT( " ( R{} ) = ( {:4.2e} ) ; ", coupledSolverAttributePrefix(), residualNorm ); + std::cout << GEOS_FMT( " ( R{} ) = ( {:4.2e} )", coupledSolverAttributePrefix(), residualNorm ); } return residualNorm; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/ReactiveCompositionalMultiphaseOBL.cpp b/src/coreComponents/physicsSolvers/fluidFlow/ReactiveCompositionalMultiphaseOBL.cpp index 23aa73586f0..772a39586c4 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/ReactiveCompositionalMultiphaseOBL.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/ReactiveCompositionalMultiphaseOBL.cpp @@ -364,7 +364,10 @@ real64 ReactiveCompositionalMultiphaseOBL::calculateResidualNorm( real64 const & real64 const residual = m_useDARTSL2Norm ? MpiWrapper::max( localResidualNorm ) : std::sqrt( MpiWrapper::sum( localResidualNorm ) ); - GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( " ( Rflow ) = ( {:4.2e} ) ;", residual ) ); + if( getLogLevel() >= 1 && logger::internal::rank==0 ) + { + std::cout << GEOS_FMT( " ( Rflow ) = ( {:4.2e} )", residual ); + } return residual; } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp index eecbd88e191..3a4aa79bb10 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp @@ -228,7 +228,7 @@ real64 SinglePhaseFVM< BASE >::calculateResidualNorm( real64 const & GEOS_UNUSED if( getLogLevel() >= 1 && logger::internal::rank == 0 ) { - std::cout << GEOS_FMT( " ( R{} ) = ( {:4.2e} ) ; ( Renergy ) = ( {:4.2e} ) ; ", + std::cout << GEOS_FMT( " ( R{} ) = ( {:4.2e} ) ( Renergy ) = ( {:4.2e} )", FlowSolverBase::coupledSolverAttributePrefix(), globalResidualNorm[0], globalResidualNorm[1] ); } } @@ -246,7 +246,7 @@ real64 SinglePhaseFVM< BASE >::calculateResidualNorm( real64 const & GEOS_UNUSED if( getLogLevel() >= 1 && logger::internal::rank == 0 ) { - std::cout << GEOS_FMT( " ( R{} ) = ( {:4.2e} ) ; ", FlowSolverBase::coupledSolverAttributePrefix(), residualNorm ); + std::cout << GEOS_FMT( " ( R{} ) = ( {:4.2e} )", FlowSolverBase::coupledSolverAttributePrefix(), residualNorm ); } } return residualNorm; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.cpp index a67553034c3..f9096ef64d3 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.cpp @@ -561,7 +561,7 @@ real64 SinglePhaseHybridFVM::calculateResidualNorm( real64 const & GEOS_UNUSED_P if( getLogLevel() >= 1 && logger::internal::rank == 0 ) { - std::cout << GEOS_FMT( " ( R{} ) = ( {:4.2e} ) ; ", coupledSolverAttributePrefix(), residualNorm ); + std::cout << GEOS_FMT( " ( R{} ) = ( {:4.2e} )", coupledSolverAttributePrefix(), residualNorm ); } return residualNorm; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/proppantTransport/ProppantTransport.cpp b/src/coreComponents/physicsSolvers/fluidFlow/proppantTransport/ProppantTransport.cpp index 3b69eaec290..780d9bc5d17 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/proppantTransport/ProppantTransport.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/proppantTransport/ProppantTransport.cpp @@ -489,11 +489,10 @@ void ProppantTransport::implicitStepComplete( real64 const & GEOS_UNUSED_PARAM( void ProppantTransport::setupDofs( DomainPartition const & GEOS_UNUSED_PARAM( domain ), DofManager & dofManager ) const { - - for( auto const & meshTarget : getMeshTargets() ) + for( auto const & meshTarget: getMeshTargets()) { - printf( "(%s,%s):", meshTarget.first.first.c_str(), meshTarget.first.second.c_str() ); - std::cout<= 1 && logger::internal::rank == 0 ) { - std::cout << GEOS_FMT( " ( R{} ) = ( {:4.2e} ) ; ", ProppantTransport::coupledSolverAttributePrefix(), residualNorm ); + std::cout << GEOS_FMT( " ( R{} ) = ( {:4.2e} )", ProppantTransport::coupledSolverAttributePrefix(), residualNorm ); } return residualNorm; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp index 88f6e4dd276..045e9aee794 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp @@ -235,6 +235,25 @@ void CompositionalMultiphaseWell::registerDataOnMesh( Group & meshBodies ) reference().resizeDimension< 0 >( m_numComponents ); wellControls.registerWrapper< real64 >( viewKeyStruct::dCurrentTotalVolRate_dRateString() ). setRestartFlags( RestartFlags::NO_WRITE ); + + // write rates output header + // the rank that owns the reference well element is responsible + if( getLogLevel() > 0 && subRegion.isLocallyOwned() ) + { + string const wellControlsName = wellControls.getName(); + string const massUnit = m_useMass ? "kg" : "mol"; + integer const useSurfaceConditions = wellControls.useSurfaceConditions(); + string const conditionKey = useSurfaceConditions ? "surface" : "reservoir"; + string const unitKey = useSurfaceConditions ? "s" : "r"; + integer const numPhase = m_numPhases; + // format: time,bhp,total_rate,total_vol_rate,phase0_vol_rate,phase1_vol_rate,... + std::ofstream outputFile( m_ratesOutputDir + "/" + wellControlsName + ".csv" ); + outputFile << "time [s],bhp [Pa],total rate [" << massUnit << "/s],total " << conditionKey << " volumetric rate [" << unitKey << "m3/s]"; + for( integer ip = 0; ip < numPhase; ++ip ) + outputFile << ",phase" << ip << " " << conditionKey << " volumetric rate [" << unitKey << "m3/s]"; + outputFile << std::endl; + outputFile.close(); + } } ); } ); @@ -1208,7 +1227,7 @@ CompositionalMultiphaseWell::calculateResidualNorm( real64 const & time_n, if( getLogLevel() >= 1 && logger::internal::rank == 0 ) { - std::cout << GEOS_FMT( " ( R{} ) = ( {:4.2e} ) ; ", coupledSolverAttributePrefix(), residualNorm ); + std::cout << GEOS_FMT( " ( R{} ) = ( {:4.2e} )", coupledSolverAttributePrefix(), residualNorm ); } return residualNorm; } @@ -1794,6 +1813,16 @@ void CompositionalMultiphaseWell::implicitStepComplete( real64 const & time_n, { WellSolverBase::implicitStepComplete( time_n, dt, domain ); + if( getLogLevel() > 0 ) + { + printRates( time_n, dt, domain ); + } +} + +void CompositionalMultiphaseWell::printRates( real64 const & time_n, + real64 const & dt, + DomainPartition & domain ) +{ forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, MeshLevel & mesh, arrayView1d< string const > const & regionNames ) @@ -1813,29 +1842,36 @@ void CompositionalMultiphaseWell::implicitStepComplete( real64 const & time_n, } integer const numPhase = m_numPhases; - localIndex const iwelemRef = subRegion.getTopWellElementIndex(); - string const massUnit = m_useMass ? "kg" : "mol"; - - // subRegion data - - arrayView1d< real64 const > const & connRate = - subRegion.getField< fields::well::mixtureConnectionRate >(); // control data WellControls const & wellControls = getWellControls( subRegion ); string const wellControlsName = wellControls.getName(); - if( wellControls.getLogLevel() == 0 ) - { - return; - } + // format: time,total_rate,total_vol_rate,phase0_vol_rate,phase1_vol_rate,... + std::ofstream outputFile( m_ratesOutputDir + "/" + wellControlsName + ".csv", std::ios_base::app ); + + outputFile << time_n + dt; + if( !wellControls.isWellOpen( time_n + dt ) ) { GEOS_LOG( GEOS_FMT( "{}: well is shut", wellControlsName ) ); + // print all zeros in the rates file + outputFile << "0.0,0.0,0.0"; + for( integer ip = 0; ip < numPhase; ++ip ) + outputFile << ",0.0"; + outputFile << std::endl; return; } + localIndex const iwelemRef = subRegion.getTopWellElementIndex(); + string const massUnit = m_useMass ? "kg" : "mol"; + + // subRegion data + + arrayView1d< real64 const > const & connRate = + subRegion.getField< fields::well::mixtureConnectionRate >(); + integer const useSurfaceConditions = wellControls.useSurfaceConditions(); real64 const & currentBHP = @@ -1854,7 +1890,8 @@ void CompositionalMultiphaseWell::implicitStepComplete( real64 const & time_n, currentPhaseVolRate, &iwelemRef, &wellControlsName, - &massUnit] ( localIndex const ) + &massUnit, + &outputFile] ( localIndex const ) { string const conditionKey = useSurfaceConditions ? "surface" : "reservoir"; string const unitKey = useSurfaceConditions ? "s" : "r"; @@ -1862,15 +1899,19 @@ void CompositionalMultiphaseWell::implicitStepComplete( real64 const & time_n, real64 const currentTotalRate = connRate[iwelemRef]; GEOS_LOG( GEOS_FMT( "{}: BHP (at the specified reference elevation): {} Pa", wellControlsName, currentBHP ) ); + outputFile << "," << currentBHP; GEOS_LOG( GEOS_FMT( "{}: Total rate: {} {}/s; total {} volumetric rate: {} {}m3/s", wellControlsName, currentTotalRate, massUnit, conditionKey, currentTotalVolRate, unitKey ) ); + outputFile << "," << currentTotalRate << "," << currentTotalVolRate; for( integer ip = 0; ip < numPhase; ++ip ) { GEOS_LOG( GEOS_FMT( "{}: Phase {} {} volumetric rate: {} {}m3/s", wellControlsName, ip, conditionKey, currentPhaseVolRate[ip], unitKey ) ); + outputFile << "," << currentPhaseVolRate[ip]; } - + outputFile << std::endl; } ); + outputFile.close(); } ); } ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp index b6f156752ca..9154e28d0b7 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp @@ -312,8 +312,6 @@ class CompositionalMultiphaseWell : public WellSolverBase protected: - - virtual void postProcessInput() override; virtual void initializePostSubGroups() override; @@ -346,6 +344,10 @@ class CompositionalMultiphaseWell : public WellSolverBase real64 const & dt, WellElementSubRegion const & subRegion ); + void printRates( real64 const & time_n, + real64 const & dt, + DomainPartition & domain ) override; + private: /** diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp index fee24395952..e72ab845bb1 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp @@ -89,6 +89,19 @@ void SinglePhaseWell::registerDataOnMesh( Group & meshBodies ) wellControls.registerWrapper< real64 >( viewKeyStruct::dCurrentVolRate_dPresString() ); wellControls.registerWrapper< real64 >( viewKeyStruct::dCurrentVolRate_dRateString() ); + // write rates output header + if( getLogLevel() > 0 ) + { + string const wellControlsName = wellControls.getName(); + integer const useSurfaceConditions = wellControls.useSurfaceConditions(); + string const conditionKey = useSurfaceConditions ? "surface" : "reservoir"; + string const unitKey = useSurfaceConditions ? "s" : "r"; + // format: time,bhp,total_rate,total_vol_rate + std::ofstream outputFile( m_ratesOutputDir + "/" + wellControlsName + ".csv" ); + outputFile << "time [s],bhp [Pa],total rate [kg/s],total " << conditionKey << " volumetric rate ["<( viewKeyStruct::fluidNamesString() ); fluidName = getConstitutiveName< SingleFluidBase >( subRegion ); GEOS_ERROR_IF( fluidName.empty(), GEOS_FMT( "{}: Fluid model not found on subregion {}", @@ -299,7 +312,7 @@ void SinglePhaseWell::updateVolRateForConstraint( WellElementSubRegion & subRegi GEOS_LOG_RANK( GEOS_FMT( "{}: The total fluid density at surface conditions is {} kg/sm3. \n" "The total rate is {} kg/s, which corresponds to a total surface volumetric rate of {} sm3/s", wellControlsName, dens[iwelemRef][0], - currentVolRate, currentVolRate ) ); + connRate[iwelemRef], currentVolRate ) ); } } ); } ); @@ -827,7 +840,7 @@ SinglePhaseWell::calculateResidualNorm( real64 const & time_n, if( getLogLevel() >= 1 && logger::internal::rank == 0 ) { - std::cout << GEOS_FMT( " ( R{} ) = ( {:4.2e} ) ; ", coupledSolverAttributePrefix(), residualNorm ); + std::cout << GEOS_FMT( " ( R{} ) = ( {:4.2e} )", coupledSolverAttributePrefix(), residualNorm ); } return residualNorm; } @@ -995,6 +1008,16 @@ void SinglePhaseWell::implicitStepComplete( real64 const & time_n, { WellSolverBase::implicitStepComplete( time_n, dt, domain ); + if( getLogLevel() > 0 ) + { + printRates( time_n, dt, domain ); + } +} + +void SinglePhaseWell::printRates( real64 const & time_n, + real64 const & dt, + DomainPartition & domain ) +{ forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, MeshLevel & mesh, arrayView1d< string const > const & regionNames ) @@ -1025,13 +1048,16 @@ void SinglePhaseWell::implicitStepComplete( real64 const & time_n, WellControls const & wellControls = getWellControls( subRegion ); string const wellControlsName = wellControls.getName(); - if( wellControls.getLogLevel() == 0 ) - { - return; - } + // format: time,total_rate,total_vol_rate + std::ofstream outputFile( m_ratesOutputDir + "/" + wellControlsName + ".csv", std::ios_base::app ); + + outputFile << time_n + dt; + if( !wellControls.isWellOpen( time_n + dt ) ) { GEOS_LOG( GEOS_FMT( "{}: well is shut", wellControlsName ) ); + // print all zeros in the rates file + outputFile << "0.0,0.0,0.0" << std::endl; return; } @@ -1048,7 +1074,8 @@ void SinglePhaseWell::implicitStepComplete( real64 const & time_n, connRate, ¤tTotalVolRate, &iwelemRef, - &wellControlsName] ( localIndex const ) + &wellControlsName, + &outputFile] ( localIndex const ) { string const conditionKey = useSurfaceConditions ? "surface" : "reservoir"; string const unitKey = useSurfaceConditions ? "s" : "r"; @@ -1056,13 +1083,14 @@ void SinglePhaseWell::implicitStepComplete( real64 const & time_n, real64 const currentTotalRate = connRate[iwelemRef]; GEOS_LOG( GEOS_FMT( "{}: BHP (at the specified reference elevation): {} Pa", wellControlsName, currentBHP ) ); + outputFile << "," << currentBHP; GEOS_LOG( GEOS_FMT( "{}: Total rate: {} kg/s; total {} volumetric rate: {} {}m3/s", wellControlsName, currentTotalRate, conditionKey, currentTotalVolRate, unitKey ) ); + outputFile << "," << currentTotalRate << "," << currentTotalVolRate << std::endl; } ); } ); } ); } - REGISTER_CATALOG_ENTRY( SolverBase, SinglePhaseWell, string const &, Group * const ) }// namespace geos diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.hpp index acf604c34fe..ae71ca23f23 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.hpp @@ -254,6 +254,10 @@ class SinglePhaseWell : public WellSolverBase virtual void initializePostSubGroups() override; + void printRates( real64 const & time_n, + real64 const & dt, + DomainPartition & domain ) override; + private: /** diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.cpp index 6caff4769b1..cbf5fdf4e5d 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.cpp @@ -37,7 +37,8 @@ WellSolverBase::WellSolverBase( string const & name, Group * const parent ) : SolverBase( name, parent ), m_numDofPerWellElement( 0 ), - m_numDofPerResElement( 0 ) + m_numDofPerResElement( 0 ), + m_ratesOutputDir( name + "_rates" ) { this->getWrapper< string >( viewKeyStruct::discretizationString() ). setInputFlag( InputFlags::FALSE ); @@ -69,6 +70,17 @@ WellSolverBase::~WellSolverBase() = default; void WellSolverBase::postProcessInput() { SolverBase::postProcessInput(); + + // create dir for rates output + if( getLogLevel() > 0 ) + { + if( MpiWrapper::commRank() == 0 ) + { + makeDirsForPath( m_ratesOutputDir ); + } + // wait till the dir is created by rank 0 + MPI_Barrier( MPI_COMM_WORLD ); + } } void WellSolverBase::registerDataOnMesh( Group & meshBodies ) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.hpp index b5d55f0bf5b..be506667005 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.hpp @@ -299,6 +299,9 @@ class WellSolverBase : public SolverBase */ virtual void initializeWells( DomainPartition & domain ) = 0; + virtual void printRates( real64 const & time_n, + real64 const & dt, + DomainPartition & domain ) = 0; /// name of the flow solver string m_flowSolverName; @@ -309,6 +312,8 @@ class WellSolverBase : public SolverBase /// the number of Degrees of Freedom per reservoir element integer m_numDofPerResElement; + string const m_ratesOutputDir; + }; } diff --git a/src/coreComponents/physicsSolvers/multiphysics/CoupledSolver.hpp b/src/coreComponents/physicsSolvers/multiphysics/CoupledSolver.hpp index 124e9b81d4e..7c66bef98d9 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/CoupledSolver.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/CoupledSolver.hpp @@ -78,12 +78,15 @@ class CoupledSolver : public SolverBase { using SolverPtr = TYPEOFREF( solver ); using SolverType = TYPEOFPTR( SolverPtr {} ); - solver = this->getParent().template getGroupPointer< SolverType >( m_names[idx()] ); + auto const & solverName = m_names[idx()]; + auto const & solverType = LvArray::system::demangleType< SolverType >(); + solver = this->getParent().template getGroupPointer< SolverType >( solverName ); GEOS_THROW_IF( solver == nullptr, GEOS_FMT( "{}: Could not find solver '{}' of type {}", getDataContext(), - m_names[idx()], LvArray::system::demangleType< SolverType >() ), + solverName, solverType ), InputError ); + GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( "{}: found {} solver named {}", getName(), solver->catalogName(), solverName ) ); } ); } @@ -479,7 +482,7 @@ class CoupledSolver : public SolverBase if( params.m_subcyclingOption == 0 ) { - GEOS_LOG_LEVEL_RANK_0( 1, "***** Single Pass solver, no subcycling *****\n" ); + GEOS_LOG_LEVEL_RANK_0( 1, "***** Single Pass solver, no subcycling *****" ); } else { @@ -523,7 +526,7 @@ class CoupledSolver : public SolverBase // finally, we perform the convergence check on the multiphysics residual residualNorm = sqrt( residualNorm ); - GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( " ( R ) = ( {:4.2e} ) ; ", residualNorm ) ); + GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( " ( R ) = ( {:4.2e} )", residualNorm ) ); isConverged = ( residualNorm < params.m_newtonTol ); } @@ -545,7 +548,7 @@ class CoupledSolver : public SolverBase if( isConverged ) { - GEOS_LOG_LEVEL_RANK_0( 1, "***** The iterative coupling has converged in " << iter + 1 << " iteration(s)! *****\n" ); + GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( "***** The iterative coupling has converged in {} iteration(s) *****", iter + 1 ) ); } } return isConverged; diff --git a/src/coreComponents/physicsSolvers/multiphysics/FlowProppantTransportSolver.cpp b/src/coreComponents/physicsSolvers/multiphysics/FlowProppantTransportSolver.cpp index c67836ad49e..17ee134ba24 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/FlowProppantTransportSolver.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/FlowProppantTransportSolver.cpp @@ -101,7 +101,7 @@ real64 FlowProppantTransportSolver::sequentiallyCoupledSolverStep( real64 const resetStateToBeginningOfStep( domain ); } - GEOS_LOG_LEVEL_RANK_0( 1, "\tIteration: " << iter+1 << ", FlowSolver: " ); + GEOS_LOG_LEVEL_RANK_0( 1, " Iteration: " << iter+1 << ", FlowSolver: " ); dtReturnTemporary = flowSolver()->nonlinearImplicitStep( time_n, dtReturn, cycleNumber, domain ); @@ -116,11 +116,11 @@ real64 FlowProppantTransportSolver::sequentiallyCoupledSolverStep( real64 const if( fluidNonLinearParams.m_numNewtonIterations <= this->m_nonlinearSolverParameters.m_minIterNewton && iter > 0 ) { m_solverStatistics.logNonlinearIteration(); - GEOS_LOG_LEVEL_RANK_0( 1, "***** The iterative coupling has converged in " << iter << " iterations! *****\n" ); + GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( "***** The iterative coupling has converged in {} iterations *****", iter ) ); break; } - GEOS_LOG_LEVEL_RANK_0( 1, "\tIteration: " << iter+1 << ", Proppant Solver: " ); + GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( " Iteration: {}, Proppant Solver: ", iter+1 ) ); dtReturnTemporary = proppantTransportSolver()->nonlinearImplicitStep( time_n, dtReturn, cycleNumber, domain ); diff --git a/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp b/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp index 47f7d477285..822c9bfac8b 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp @@ -187,6 +187,8 @@ real64 HydrofractureSolver< POROMECHANICS_SOLVER >::fullyCoupledSolverStep( real int solveIter; for( solveIter=0; solveIter= 1 && logger::internal::rank==0 ) { - std::cout << GEOS_FMT( " ( R{} ) = ( {:4.2e} ) ; ", coupledSolverAttributePrefix(), residual ); + std::cout << GEOS_FMT( " ( R{} ) = ( {:4.2e} )", coupledSolverAttributePrefix(), residual ); } return residual; diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp index 0597ecd60bc..08f91d8fe89 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp @@ -1187,7 +1187,7 @@ SolidMechanicsLagrangianFEM:: if( getLogLevel() >= 1 && logger::internal::rank==0 ) { - std::cout << GEOS_FMT( " ( R{} ) = ( {:4.2e} ) ; ", coupledSolverAttributePrefix(), totalResidualNorm ); + std::cout << GEOS_FMT( " ( R{} ) = ( {:4.2e} )", coupledSolverAttributePrefix(), totalResidualNorm ); } return totalResidualNorm; diff --git a/src/coreComponents/physicsSolvers/surfaceGeneration/SurfaceGenerator.cpp b/src/coreComponents/physicsSolvers/surfaceGeneration/SurfaceGenerator.cpp index 6da89c891bb..4cb8cc65a84 100644 --- a/src/coreComponents/physicsSolvers/surfaceGeneration/SurfaceGenerator.cpp +++ b/src/coreComponents/physicsSolvers/surfaceGeneration/SurfaceGenerator.cpp @@ -1704,15 +1704,14 @@ void SurfaceGenerator::performFracture( const localIndex nodeID, // Split the node into two, using the original index, and a new one. localIndex newNodeIndex; - if( getLogLevel() ) + if( getLogLevel() > 0 ) { - GEOS_LOG_RANK( "" ); - std::cout<<"Splitting node "<::const_iterator i=separationPathFaces.begin(); i!=separationPathFaces.end(); ++i ) { - std::cout<<*i<<", "; + s << *i << " "; } - std::cout< 0 ) + { + GEOS_LOG_RANK( GEOS_FMT( "Done splitting node {} into nodes {} and {}", nodeID, nodeID, newNodeIndex ) ); + } // split edges map< localIndex, localIndex > splitEdges; @@ -1774,10 +1775,9 @@ void SurfaceGenerator::performFracture( const localIndex nodeID, edgeToFaceMap.clearSet( newEdgeIndex ); - if( getLogLevel() ) + if( getLogLevel() > 0 ) { - GEOS_LOG_RANK( "" ); - std::cout<<" Split edge "< 0 ) { - GEOS_LOG_RANK( "" ); - std::cout<<" Split face "<= 1 ) + { + SIFNode[nodeIndex] = *min_element( SIFNode_All[nodeIndex].begin(), SIFNode_All[nodeIndex].end()); + } for( localIndex const edgeIndex: m_tipEdges ) { @@ -3236,7 +3238,10 @@ void SurfaceGenerator::calculateNodeAndFaceSif( DomainPartition const & domain, { if( m_tipFaces.contains( faceIndex )) { - SIFonFace[faceIndex] = *max_element( SIFonFace_All[faceIndex].begin(), SIFonFace_All[faceIndex].end()); + if( SIFonFace_All[faceIndex].size() >= 1 ) + { + SIFonFace[faceIndex] = *max_element( SIFonFace_All[faceIndex].begin(), SIFonFace_All[faceIndex].end()); + } } } } diff --git a/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEM.cpp index dad568c8dde..fef2cb77c33 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEM.cpp @@ -323,7 +323,7 @@ void AcousticFirstOrderWaveEquationSEM::initializePostInitialConditionsPreSubGro CellElementSubRegion & elementSubRegion ) { arrayView2d< localIndex const, cells::NODE_MAP_USD > const elemsToNodes = elementSubRegion.nodeList(); - arrayView2d< localIndex const > const facesToElements = faceManager.elementList(); + arrayView2d< localIndex const > const elemsToFaces = elementSubRegion.faceList(); arrayView1d< real32 const > const velocity = elementSubRegion.getField< wavesolverfields::MediumVelocity >(); arrayView1d< real32 const > const density = elementSubRegion.getField< wavesolverfields::MediumDensity >(); @@ -344,9 +344,9 @@ void AcousticFirstOrderWaveEquationSEM::initializePostInitialConditionsPreSubGro acousticFirstOrderWaveEquationSEMKernels::DampingMatrixKernel< FE_TYPE > kernelD( finiteElement ); - kernelD.template launch< EXEC_POLICY, ATOMIC_POLICY >( faceManager.size(), + kernelD.template launch< EXEC_POLICY, ATOMIC_POLICY >( elementSubRegion.size(), X, - facesToElements, + elemsToFaces, facesToNodes, facesDomainBoundaryIndicator, freeSurfaceFaceIndicator, @@ -382,7 +382,7 @@ void AcousticFirstOrderWaveEquationSEM::applyFreeSurfaceBC( real64 const time, D fsManager.apply< FaceManager >( time, domain.getMeshBody( 0 ).getMeshLevel( m_discretizationName ), - string( "FreeSurface" ), + WaveSolverBase::viewKeyStruct::freeSurfaceString(), [&]( FieldSpecificationBase const & bc, string const &, SortedArrayView< localIndex const > const & targetSet, diff --git a/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEMKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEMKernel.hpp index 712e909fdb5..cc01a3ba96a 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEMKernel.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEMKernel.hpp @@ -268,8 +268,8 @@ struct DampingMatrixKernel * @tparam EXEC_POLICY the execution policy * @tparam ATOMIC_POLICY the atomic policy * @param[in] size the number of cells in the subRegion - * @param[in] X coordinates of the nodes - * @param[in] facesToElems map from faces to elements + * @param[in] nodeCoords coordinates of the nodes + * @param[in] elemsToFaces map from elements to faces * @param[in] facesToNodes map from face to nodes * @param[in] facesDomainBoundaryIndicator flag equal to 1 if the face is on the boundary, and to 0 otherwise * @param[in] freeSurfaceFaceIndicator flag equal to 1 if the face is on the free surface, and to 0 otherwise @@ -279,45 +279,41 @@ struct DampingMatrixKernel template< typename EXEC_POLICY, typename ATOMIC_POLICY > void launch( localIndex const size, - arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const X, - arrayView2d< localIndex const > const facesToElems, + arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords, + arrayView2d< localIndex const > const elemsToFaces, ArrayOfArraysView< localIndex const > const facesToNodes, arrayView1d< integer const > const facesDomainBoundaryIndicator, arrayView1d< localIndex const > const freeSurfaceFaceIndicator, arrayView1d< real32 const > const velocity, arrayView1d< real32 > const damping ) { - forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const f ) + forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const e ) { - // face on the domain boundary and not on free surface - if( facesDomainBoundaryIndicator[f] == 1 && freeSurfaceFaceIndicator[f] != 1 ) + for( localIndex i = 0; i < elemsToFaces.size( 1 ); ++i ) { - localIndex k = facesToElems( f, 0 ); - if( k < 0 ) - { - k = facesToElems( f, 1 ); - } - - real32 const alpha = 1.0 / velocity[k]; - - constexpr localIndex numNodesPerFace = FE_TYPE::numNodesPerFace; - - real64 xLocal[ numNodesPerFace ][ 3 ]; - for( localIndex a = 0; a < numNodesPerFace; ++a ) + localIndex const f = elemsToFaces( e, i ); + // face on the domain boundary and not on free surface + if( facesDomainBoundaryIndicator[f] == 1 && freeSurfaceFaceIndicator[f] != 1 ) { - for( localIndex i = 0; i < 3; ++i ) + constexpr localIndex numNodesPerFace = FE_TYPE::numNodesPerFace; + real64 xLocal[ numNodesPerFace ][ 3 ]; + for( localIndex a = 0; a < numNodesPerFace; ++a ) { - xLocal[a][i] = X( facesToNodes( f, a ), i ); + for( localIndex d = 0; d < 3; ++d ) + { + xLocal[a][d] = nodeCoords( facesToNodes( f, a ), d ); + } } - } - for( localIndex q = 0; q < numNodesPerFace; ++q ) - { - real32 const localIncrement = alpha * m_finiteElement.computeDampingTerm( q, xLocal ); - RAJA::atomicAdd< ATOMIC_POLICY >( &damping[facesToNodes[f][q]], localIncrement ); + real32 const alpha = 1.0 / velocity[e]; + for( localIndex q = 0; q < numNodesPerFace; ++q ) + { + real32 const localIncrement = alpha * m_finiteElement.computeDampingTerm( q, xLocal ); + RAJA::atomicAdd< ATOMIC_POLICY >( &damping[facesToNodes( f, q )], localIncrement ); + } } } - } ); // end loop over element + } ); } /// The finite element space/discretization object for the element type in the subRegion diff --git a/src/coreComponents/physicsSolvers/wavePropagation/AcousticVTIWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/AcousticVTIWaveEquationSEM.cpp new file mode 100644 index 00000000000..e17726bd496 --- /dev/null +++ b/src/coreComponents/physicsSolvers/wavePropagation/AcousticVTIWaveEquationSEM.cpp @@ -0,0 +1,726 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 TotalEnergies + * Copyright (c) 2019- GEOS Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + + +/** + * @file AcousticVTIWaveEquationSEM.cpp + */ + +#include "AcousticVTIWaveEquationSEM.hpp" +#include "AcousticVTIWaveEquationSEMKernel.hpp" + +#include "finiteElement/FiniteElementDiscretization.hpp" +#include "fieldSpecification/FieldSpecificationManager.hpp" +#include "fieldSpecification/PerfectlyMatchedLayer.hpp" +#include "mainInterface/ProblemManager.hpp" +#include "mesh/ElementType.hpp" +#include "mesh/mpiCommunications/CommunicationTools.hpp" +#include "WaveSolverUtils.hpp" +#include "WaveSolverBaseFields.hpp" + +namespace geos +{ + +using namespace dataRepository; + +AcousticVTIWaveEquationSEM::AcousticVTIWaveEquationSEM( const std::string & name, + Group * const parent ): + WaveSolverBase( name, + parent ) +{ + + registerWrapper( viewKeyStruct::pressureNp1AtReceiversString(), &m_pressureNp1AtReceivers ). + setInputFlag( InputFlags::FALSE ). + setSizedFromParent( 0 ). + setDescription( "Pressure value at each receiver for each timestep" ); +} + +void AcousticVTIWaveEquationSEM::initializePreSubGroups() +{ + WaveSolverBase::initializePreSubGroups(); + +} + +void AcousticVTIWaveEquationSEM::registerDataOnMesh( Group & meshBodies ) +{ + WaveSolverBase::registerDataOnMesh( meshBodies ); + + forDiscretizationOnMeshTargets( meshBodies, [&] ( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & ) + { + NodeManager & nodeManager = mesh.getNodeManager(); + + nodeManager.registerField< fields::wavesolverfields::Pressure_p_nm1, + fields::wavesolverfields::Pressure_p_n, + fields::wavesolverfields::Pressure_p_np1, + fields::wavesolverfields::Pressure_q_nm1, + fields::wavesolverfields::Pressure_q_n, + fields::wavesolverfields::Pressure_q_np1, + fields::wavesolverfields::ForcingRHS, + fields::wavesolverfields::MassVector, + fields::wavesolverfields::DampingVector_p, + fields::wavesolverfields::DampingVector_pq, + fields::wavesolverfields::DampingVector_q, + fields::wavesolverfields::DampingVector_qp, + fields::wavesolverfields::StiffnessVector_p, + fields::wavesolverfields::StiffnessVector_q, + fields::wavesolverfields::FreeSurfaceNodeIndicator, + fields::wavesolverfields::LateralSurfaceNodeIndicator, + fields::wavesolverfields::BottomSurfaceNodeIndicator >( this->getName() ); + + + FaceManager & faceManager = mesh.getFaceManager(); + faceManager.registerField< fields::wavesolverfields::FreeSurfaceFaceIndicator >( this->getName() ); + faceManager.registerField< fields::wavesolverfields::LateralSurfaceFaceIndicator >( this->getName() ); + faceManager.registerField< fields::wavesolverfields::BottomSurfaceFaceIndicator >( this->getName() ); + + ElementRegionManager & elemManager = mesh.getElemManager(); + + elemManager.forElementSubRegions< CellElementSubRegion >( [&]( CellElementSubRegion & subRegion ) + { + subRegion.registerField< fields::wavesolverfields::Delta >( this->getName() ); + subRegion.registerField< fields::wavesolverfields::Epsilon >( this->getName() ); + subRegion.registerField< fields::wavesolverfields::F >( this->getName() ); + subRegion.registerField< fields::wavesolverfields::MediumVelocity >( this->getName() ); + } ); + } ); +} + + +void AcousticVTIWaveEquationSEM::postProcessInput() +{ + + WaveSolverBase::postProcessInput(); + + localIndex const numReceiversGlobal = m_receiverCoordinates.size( 0 ); + + m_pressureNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal ); +} + +void AcousticVTIWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) +{ + NodeManager const & nodeManager = mesh.getNodeManager(); + FaceManager const & faceManager = mesh.getFaceManager(); + + arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const X32 = + nodeManager.getField< fields::referencePosition32 >().toViewConst(); + + arrayView2d< real64 const > const faceNormal = faceManager.faceNormal(); + arrayView2d< real64 const > const faceCenter = faceManager.faceCenter(); + + + arrayView2d< real64 const > const sourceCoordinates = m_sourceCoordinates.toViewConst(); + arrayView2d< localIndex > const sourceNodeIds = m_sourceNodeIds.toView(); + arrayView2d< real64 > const sourceConstants = m_sourceConstants.toView(); + arrayView1d< localIndex > const sourceIsAccessible = m_sourceIsAccessible.toView(); + sourceNodeIds.setValues< EXEC_POLICY >( -1 ); + sourceConstants.setValues< EXEC_POLICY >( -1 ); + sourceIsAccessible.zero(); + + arrayView2d< real64 const > const receiverCoordinates = m_receiverCoordinates.toViewConst(); + arrayView2d< localIndex > const receiverNodeIds = m_receiverNodeIds.toView(); + arrayView2d< real64 > const receiverConstants = m_receiverConstants.toView(); + arrayView1d< localIndex > const receiverIsLocal = m_receiverIsLocal.toView(); + receiverNodeIds.setValues< EXEC_POLICY >( -1 ); + receiverConstants.setValues< EXEC_POLICY >( -1 ); + receiverIsLocal.zero(); + + real32 const timeSourceFrequency = this->m_timeSourceFrequency; + localIndex const rickerOrder = this->m_rickerOrder; + arrayView2d< real32 > const sourceValue = m_sourceValue.toView(); + real64 dt = 0; + EventManager const & event = this->getGroupByPath< EventManager >( "/Problem/Events" ); + for( localIndex numSubEvent = 0; numSubEvent < event.numSubGroups(); ++numSubEvent ) + { + EventBase const * subEvent = static_cast< EventBase const * >( event.getSubGroups()[numSubEvent] ); + if( subEvent->getEventName() == "/Solvers/" + this->getName() ) + { + dt = subEvent->getReference< real64 >( EventBase::viewKeyStruct::forceDtString() ); + } + } + + mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, + CellElementSubRegion & elementSubRegion ) + { + GEOS_THROW_IF( elementSubRegion.getElementType() != ElementType::Hexahedron, + "Invalid type of element, the acoustic solver is designed for hexahedral meshes only (C3D8), using the SEM formulation", + InputError ); + + arrayView2d< localIndex const > const elemsToFaces = elementSubRegion.faceList(); + arrayView2d< localIndex const, cells::NODE_MAP_USD > const & elemsToNodes = elementSubRegion.nodeList(); + arrayView2d< real64 const > const elemCenter = elementSubRegion.getElementCenter(); + arrayView1d< integer const > const elemGhostRank = elementSubRegion.ghostRank(); + + finiteElement::FiniteElementBase const & + fe = elementSubRegion.getReference< finiteElement::FiniteElementBase >( getDiscretizationName() ); + finiteElement::FiniteElementDispatchHandler< SEM_FE_TYPES >::dispatch3D( fe, [&] ( auto const finiteElement ) + { + using FE_TYPE = TYPEOFREF( finiteElement ); + + constexpr localIndex numNodesPerElem = FE_TYPE::numNodes; + localIndex const numFacesPerElem = elementSubRegion.numFacesPerElement(); + + acousticVTIWaveEquationSEMKernels:: + PrecomputeSourceAndReceiverKernel:: + launch< EXEC_POLICY, FE_TYPE > + ( elementSubRegion.size(), + numNodesPerElem, + numFacesPerElem, + X32, + elemGhostRank, + elemsToNodes, + elemsToFaces, + elemCenter, + faceNormal, + faceCenter, + sourceCoordinates, + sourceIsAccessible, + sourceNodeIds, + sourceConstants, + receiverCoordinates, + receiverIsLocal, + receiverNodeIds, + receiverConstants, + sourceValue, + dt, + timeSourceFrequency, + rickerOrder ); + } ); + } ); +} + +void AcousticVTIWaveEquationSEM::addSourceToRightHandSide( integer const & cycleNumber, arrayView1d< real32 > const rhs ) +{ + arrayView2d< localIndex const > const sourceNodeIds = m_sourceNodeIds.toViewConst(); + arrayView2d< real64 const > const sourceConstants = m_sourceConstants.toViewConst(); + arrayView1d< localIndex const > const sourceIsAccessible = m_sourceIsAccessible.toViewConst(); + arrayView2d< real32 const > const sourceValue = m_sourceValue.toViewConst(); + + GEOS_THROW_IF( cycleNumber > sourceValue.size( 0 ), "Too many steps compared to array size", std::runtime_error ); + forAll< EXEC_POLICY >( sourceConstants.size( 0 ), [=] GEOS_HOST_DEVICE ( localIndex const isrc ) + { + if( sourceIsAccessible[isrc] == 1 ) + { + for( localIndex inode = 0; inode < sourceConstants.size( 1 ); ++inode ) + { + real32 const localIncrement = sourceConstants[isrc][inode] * sourceValue[cycleNumber][isrc]; + RAJA::atomicAdd< ATOMIC_POLICY >( &rhs[sourceNodeIds[isrc][inode]], localIncrement ); + } + } + } ); +} + +void AcousticVTIWaveEquationSEM::initializePostInitialConditionsPreSubGroups() +{ + + WaveSolverBase::initializePostInitialConditionsPreSubGroups(); + + DomainPartition & domain = this->getGroupByPath< DomainPartition >( "/Problem/domain" ); + + real64 const time = 0.0; + applyFreeSurfaceBC( time, domain ); + precomputeSurfaceFieldIndicator( domain ); + + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) + { + precomputeSourceAndReceiverTerm( mesh, regionNames ); + + NodeManager & nodeManager = mesh.getNodeManager(); + FaceManager & faceManager = mesh.getFaceManager(); + + /// get the array of indicators: 1 if the face is on the boundary; 0 otherwise + arrayView1d< integer > const & facesDomainBoundaryIndicator = faceManager.getDomainBoundaryIndicator(); + arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const X32 = + nodeManager.getField< fields::referencePosition32 >().toViewConst(); + + /// get face to node map + ArrayOfArraysView< localIndex const > const facesToNodes = faceManager.nodeList().toViewConst(); + + // mass matrix to be computed in this function + arrayView1d< real32 > const mass = nodeManager.getField< fields::wavesolverfields::MassVector >(); + mass.zero(); + /// damping matrices to be computed for each dof in the boundary of the mesh + arrayView1d< real32 > const damping_p = nodeManager.getField< fields::wavesolverfields::DampingVector_p >(); + arrayView1d< real32 > const damping_pq = nodeManager.getField< fields::wavesolverfields::DampingVector_pq >(); + arrayView1d< real32 > const damping_q = nodeManager.getField< fields::wavesolverfields::DampingVector_q >(); + arrayView1d< real32 > const damping_qp = nodeManager.getField< fields::wavesolverfields::DampingVector_qp >(); + damping_p.zero(); + damping_pq.zero(); + damping_q.zero(); + damping_qp.zero(); + + /// get array of indicators: 1 if face is on the free surface; 0 otherwise + arrayView1d< localIndex const > const freeSurfaceFaceIndicator = faceManager.getField< fields::wavesolverfields::FreeSurfaceFaceIndicator >(); + arrayView1d< localIndex const > const lateralSurfaceFaceIndicator = faceManager.getField< fields::wavesolverfields::LateralSurfaceFaceIndicator >(); + arrayView1d< localIndex const > const bottomSurfaceFaceIndicator = faceManager.getField< fields::wavesolverfields::BottomSurfaceFaceIndicator >(); + + mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, + CellElementSubRegion & elementSubRegion ) + { + + arrayView2d< localIndex const, cells::NODE_MAP_USD > const elemsToNodes = elementSubRegion.nodeList(); + arrayView2d< localIndex const > const facesToElements = faceManager.elementList(); + arrayView1d< real32 const > const velocity = elementSubRegion.getField< fields::wavesolverfields::MediumVelocity >(); + arrayView1d< real32 const > const epsilon = elementSubRegion.getField< fields::wavesolverfields::Epsilon >(); + arrayView1d< real32 const > const delta = elementSubRegion.getField< fields::wavesolverfields::Delta >(); + arrayView1d< real32 const > const vti_f = elementSubRegion.getField< fields::wavesolverfields::F >(); + + finiteElement::FiniteElementBase const & + fe = elementSubRegion.getReference< finiteElement::FiniteElementBase >( getDiscretizationName() ); + finiteElement::FiniteElementDispatchHandler< SEM_FE_TYPES >::dispatch3D( fe, [&] ( auto const finiteElement ) + { + using FE_TYPE = TYPEOFREF( finiteElement ); + + acousticVTIWaveEquationSEMKernels::MassMatrixKernel< FE_TYPE > kernelM( finiteElement ); + + kernelM.template launch< EXEC_POLICY, ATOMIC_POLICY >( elementSubRegion.size(), + X32, + elemsToNodes, + velocity, + mass ); + + acousticVTIWaveEquationSEMKernels::DampingMatrixKernel< FE_TYPE > kernelD( finiteElement ); + + kernelD.template launch< EXEC_POLICY, ATOMIC_POLICY >( faceManager.size(), + X32, + facesToElements, + facesToNodes, + facesDomainBoundaryIndicator, + freeSurfaceFaceIndicator, + lateralSurfaceFaceIndicator, + bottomSurfaceFaceIndicator, + velocity, + epsilon, + delta, + vti_f, + damping_p, + damping_q, + damping_pq, + damping_qp ); + } ); + } ); + } ); + +} + +void AcousticVTIWaveEquationSEM::precomputeSurfaceFieldIndicator( DomainPartition & domain ) +{ + real64 const time = 0.0; + FieldSpecificationManager & fsManager = FieldSpecificationManager::getInstance(); + FunctionManager const & functionManager = FunctionManager::getInstance(); + + FaceManager & faceManager = domain.getMeshBody( 0 ).getMeshLevel( m_discretizationName ).getFaceManager(); + NodeManager & nodeManager = domain.getMeshBody( 0 ).getMeshLevel( m_discretizationName ).getNodeManager(); + + ArrayOfArraysView< localIndex const > const faceToNodeMap = faceManager.nodeList().toViewConst(); + + /// array of indicators: 1 if a face is on on lateral surface; 0 otherwise + arrayView1d< localIndex > const lateralSurfaceFaceIndicator = faceManager.getField< fields::wavesolverfields::LateralSurfaceFaceIndicator >(); + /// array of indicators: 1 if a node is on on lateral surface; 0 otherwise + arrayView1d< localIndex > const lateralSurfaceNodeIndicator = nodeManager.getField< fields::wavesolverfields::LateralSurfaceNodeIndicator >(); + + /// array of indicators: 1 if a face is on on bottom surface; 0 otherwise + arrayView1d< localIndex > const bottomSurfaceFaceIndicator = faceManager.getField< fields::wavesolverfields::BottomSurfaceFaceIndicator >(); + /// array of indicators: 1 if a node is on on bottom surface; 0 otherwise + arrayView1d< localIndex > const bottomSurfaceNodeIndicator = nodeManager.getField< fields::wavesolverfields::BottomSurfaceNodeIndicator >(); + + // Lateral surfaces + fsManager.apply< FaceManager >( time, + domain.getMeshBody( 0 ).getMeshLevel( m_discretizationName ), + viewKeyStruct::lateralSurfaceString(), + [&]( FieldSpecificationBase const & bc, + string const &, + SortedArrayView< localIndex const > const & targetSet, + FaceManager &, + string const & ) + { + string const & functionName = bc.getFunctionName(); + + if( functionName.empty() || functionManager.getGroup< FunctionBase >( functionName ).isFunctionOfTime() == 2 ) + { + for( localIndex i = 0; i < targetSet.size(); ++i ) + { + localIndex const kf = targetSet[ i ]; + lateralSurfaceFaceIndicator[kf] = 1; + + localIndex const numNodes = faceToNodeMap.sizeOfArray( kf ); + for( localIndex a=0; a < numNodes; ++a ) + { + localIndex const dof = faceToNodeMap( kf, a ); + lateralSurfaceNodeIndicator[dof] = 1; + } + } + } + else + { + GEOS_ERROR( "This option is not supported yet" ); + } + } ); + + // For the Bottom surface + fsManager.apply< FaceManager >( time, + domain.getMeshBody( 0 ).getMeshLevel( m_discretizationName ), + viewKeyStruct::bottomSurfaceString(), + [&]( FieldSpecificationBase const & bc, + string const &, + SortedArrayView< localIndex const > const & targetSet, + FaceManager &, + string const & ) + { + string const & functionName = bc.getFunctionName(); + + if( functionName.empty() || functionManager.getGroup< FunctionBase >( functionName ).isFunctionOfTime() == 2 ) + { + for( localIndex i = 0; i < targetSet.size(); ++i ) + { + localIndex const kf = targetSet[ i ]; + bottomSurfaceFaceIndicator[kf] = 1; + + localIndex const numNodes = faceToNodeMap.sizeOfArray( kf ); + for( localIndex a=0; a < numNodes; ++a ) + { + localIndex const dof = faceToNodeMap( kf, a ); + bottomSurfaceNodeIndicator[dof] = 1; + } + } + } + else + { + GEOS_ERROR( "This option is not supported yet" ); + } + } ); +} + +void AcousticVTIWaveEquationSEM::applyFreeSurfaceBC( real64 time, DomainPartition & domain ) +{ + FieldSpecificationManager & fsManager = FieldSpecificationManager::getInstance(); + FunctionManager const & functionManager = FunctionManager::getInstance(); + + FaceManager & faceManager = domain.getMeshBody( 0 ).getMeshLevel( m_discretizationName ).getFaceManager(); + NodeManager & nodeManager = domain.getMeshBody( 0 ).getMeshLevel( m_discretizationName ).getNodeManager(); + + arrayView1d< real32 > const p_nm1 = nodeManager.getField< fields::wavesolverfields::Pressure_p_nm1 >(); + arrayView1d< real32 > const p_n = nodeManager.getField< fields::wavesolverfields::Pressure_p_n >(); + arrayView1d< real32 > const p_np1 = nodeManager.getField< fields::wavesolverfields::Pressure_p_np1 >(); + + arrayView1d< real32 > const q_nm1 = nodeManager.getField< fields::wavesolverfields::Pressure_q_nm1 >(); + arrayView1d< real32 > const q_n = nodeManager.getField< fields::wavesolverfields::Pressure_q_n >(); + arrayView1d< real32 > const q_np1 = nodeManager.getField< fields::wavesolverfields::Pressure_q_np1 >(); + + ArrayOfArraysView< localIndex const > const faceToNodeMap = faceManager.nodeList().toViewConst(); + + /// array of indicators: 1 if a face is on on free surface; 0 otherwise + arrayView1d< localIndex > const freeSurfaceFaceIndicator = faceManager.getField< fields::wavesolverfields::FreeSurfaceFaceIndicator >(); + + /// array of indicators: 1 if a node is on on free surface; 0 otherwise + arrayView1d< localIndex > const freeSurfaceNodeIndicator = nodeManager.getField< fields::wavesolverfields::FreeSurfaceNodeIndicator >(); + + fsManager.apply< FaceManager >( time, + domain.getMeshBody( 0 ).getMeshLevel( m_discretizationName ), + WaveSolverBase::viewKeyStruct::freeSurfaceString(), + [&]( FieldSpecificationBase const & bc, + string const &, + SortedArrayView< localIndex const > const & targetSet, + FaceManager &, + string const & ) + { + string const & functionName = bc.getFunctionName(); + + if( functionName.empty() || functionManager.getGroup< FunctionBase >( functionName ).isFunctionOfTime() == 2 ) + { + real64 const value = bc.getScale(); + + for( localIndex i = 0; i < targetSet.size(); ++i ) + { + localIndex const kf = targetSet[ i ]; + freeSurfaceFaceIndicator[kf] = 1; + + localIndex const numNodes = faceToNodeMap.sizeOfArray( kf ); + for( localIndex a=0; a < numNodes; ++a ) + { + localIndex const dof = faceToNodeMap( kf, a ); + freeSurfaceNodeIndicator[dof] = 1; + + p_np1[dof] = value; + p_n[dof] = value; + p_nm1[dof] = value; + + q_np1[dof] = value; + q_n[dof] = value; + q_nm1[dof] = value; + } + } + } + else + { + GEOS_ERROR( "This option is not supported yet" ); + } + } ); +} + +real64 AcousticVTIWaveEquationSEM::explicitStepForward( real64 const & time_n, + real64 const & dt, + integer cycleNumber, + DomainPartition & domain, + bool computeGradient ) +{ + real64 dtOut = explicitStepInternal( time_n, dt, cycleNumber, domain ); + + forDiscretizationOnMeshTargets( domain.getMeshBodies(), + [&] ( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & GEOS_UNUSED_PARAM ( regionNames ) ) + { + NodeManager & nodeManager = mesh.getNodeManager(); + + arrayView1d< real32 > const p_nm1 = nodeManager.getField< fields::wavesolverfields::Pressure_p_nm1 >(); + arrayView1d< real32 > const p_n = nodeManager.getField< fields::wavesolverfields::Pressure_p_n >(); + arrayView1d< real32 > const p_np1 = nodeManager.getField< fields::wavesolverfields::Pressure_p_np1 >(); + + arrayView1d< real32 > const q_nm1 = nodeManager.getField< fields::wavesolverfields::Pressure_q_nm1 >(); + arrayView1d< real32 > const q_n = nodeManager.getField< fields::wavesolverfields::Pressure_q_n >(); + arrayView1d< real32 > const q_np1 = nodeManager.getField< fields::wavesolverfields::Pressure_q_np1 >(); + + if( computeGradient ) + { + GEOS_ERROR( "This option is not supported yet" ); + } + + forAll< EXEC_POLICY >( nodeManager.size(), [=] GEOS_HOST_DEVICE ( localIndex const a ) + { + p_nm1[a] = p_n[a]; + p_n[a] = p_np1[a]; + + q_nm1[a] = q_n[a]; + q_n[a] = q_np1[a]; + + } ); + } ); + + return dtOut; +} + +void AcousticVTIWaveEquationSEM::initializePML() +{ + GEOS_ERROR( "This option is not supported yet" ); + return; +} + +void AcousticVTIWaveEquationSEM::applyPML( real64 const GEOS_UNUSED_PARAM( time ), + DomainPartition & GEOS_UNUSED_PARAM( domain )) +{ + GEOS_ERROR( "This option is not supported yet" ); + return; +} + +real64 AcousticVTIWaveEquationSEM::explicitStepBackward( real64 const & GEOS_UNUSED_PARAM( time_n ), + real64 const & GEOS_UNUSED_PARAM( dt ), + integer GEOS_UNUSED_PARAM( cycleNumber ), + DomainPartition & GEOS_UNUSED_PARAM( domain ), + bool GEOS_UNUSED_PARAM( computeGradient ) ) +{ + GEOS_ERROR( "This option is not supported yet" ); + return -1; +} + +real64 AcousticVTIWaveEquationSEM::explicitStepInternal( real64 const & time_n, + real64 const & dt, + integer cycleNumber, + DomainPartition & domain ) +{ + GEOS_MARK_FUNCTION; + + GEOS_LOG_RANK_0_IF( dt < epsilonLoc, "Warning! Value for dt: " << dt << "s is smaller than local threshold: " << epsilonLoc ); + + forDiscretizationOnMeshTargets( domain.getMeshBodies(), + [&] ( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) + { + NodeManager & nodeManager = mesh.getNodeManager(); + + arrayView1d< real32 const > const mass = nodeManager.getField< fields::wavesolverfields::MassVector >(); + arrayView1d< real32 const > const damping_p = nodeManager.getField< fields::wavesolverfields::DampingVector_p >(); + arrayView1d< real32 const > const damping_q = nodeManager.getField< fields::wavesolverfields::DampingVector_q >(); + arrayView1d< real32 const > const damping_pq = nodeManager.getField< fields::wavesolverfields::DampingVector_pq >(); + arrayView1d< real32 const > const damping_qp = nodeManager.getField< fields::wavesolverfields::DampingVector_qp >(); + + arrayView1d< real32 > const p_nm1 = nodeManager.getField< fields::wavesolverfields::Pressure_p_nm1 >(); + arrayView1d< real32 > const p_n = nodeManager.getField< fields::wavesolverfields::Pressure_p_n >(); + arrayView1d< real32 > const p_np1 = nodeManager.getField< fields::wavesolverfields::Pressure_p_np1 >(); + + arrayView1d< real32 > const q_nm1 = nodeManager.getField< fields::wavesolverfields::Pressure_q_nm1 >(); + arrayView1d< real32 > const q_n = nodeManager.getField< fields::wavesolverfields::Pressure_q_n >(); + arrayView1d< real32 > const q_np1 = nodeManager.getField< fields::wavesolverfields::Pressure_q_np1 >(); + + arrayView1d< localIndex const > const freeSurfaceNodeIndicator = nodeManager.getField< fields::wavesolverfields::FreeSurfaceNodeIndicator >(); + arrayView1d< localIndex const > const lateralSurfaceNodeIndicator = nodeManager.getField< fields::wavesolverfields::LateralSurfaceNodeIndicator >(); + arrayView1d< localIndex const > const bottomSurfaceNodeIndicator = nodeManager.getField< fields::wavesolverfields::BottomSurfaceNodeIndicator >(); + arrayView1d< real32 > const stiffnessVector_p = nodeManager.getField< fields::wavesolverfields::StiffnessVector_p >(); + arrayView1d< real32 > const stiffnessVector_q = nodeManager.getField< fields::wavesolverfields::StiffnessVector_q >(); + arrayView1d< real32 > const rhs = nodeManager.getField< fields::wavesolverfields::ForcingRHS >(); + + auto kernelFactory = acousticVTIWaveEquationSEMKernels::ExplicitAcousticVTISEMFactory( dt ); + + finiteElement:: + regionBasedKernelApplication< EXEC_POLICY, + constitutive::NullModel, + CellElementSubRegion >( mesh, + regionNames, + getDiscretizationName(), + "", + kernelFactory ); + + addSourceToRightHandSide( cycleNumber, rhs ); + + /// calculate your time integrators + real64 const dt2 = dt*dt; + + GEOS_MARK_SCOPE ( updateP ); + forAll< EXEC_POLICY >( nodeManager.size(), [=] GEOS_HOST_DEVICE ( localIndex const a ) + { + if( freeSurfaceNodeIndicator[a] != 1 ) + { + p_np1[a] = 2.0*mass[a]*p_n[a]/dt2; + p_np1[a] -= mass[a]*p_nm1[a]/dt2; + p_np1[a] += stiffnessVector_p[a]; + p_np1[a] += rhs[a]; + + q_np1[a] = 2.0*mass[a]*q_n[a]/dt2; + q_np1[a] -= mass[a]*q_nm1[a]/dt2; + q_np1[a] += stiffnessVector_q[a]; + q_np1[a] += rhs[a]; + + if( lateralSurfaceNodeIndicator[a] != 1 && bottomSurfaceNodeIndicator[a] != 1 ) + { + // Interior node, no boundary terms + p_np1[a] /= mass[a]/dt2; + q_np1[a] /= mass[a]/dt2; + } + else + { + // Boundary node + p_np1[a] += damping_p[a]*p_nm1[a]/dt/2; + p_np1[a] += damping_pq[a]*q_nm1[a]/dt/2; + + q_np1[a] += damping_q[a]*q_nm1[a]/dt/2; + q_np1[a] += damping_qp[a]*p_nm1[a]/dt/2; + // Hand-made Inversion of 2x2 matrix + real32 coef_pp = mass[a]/dt2; + coef_pp += damping_p[a]/dt/2; + real32 coef_pq = damping_pq[a]/dt/2; + + real32 coef_qq = mass[a]/dt2; + coef_qq += damping_q[a]/2/dt; + real32 coef_qp = damping_qp[a]/dt/2; + + real32 det_pq = 1/(coef_pp * coef_qq - coef_pq*coef_qp); + + real32 aux_p_np1 = p_np1[a]; + p_np1[a] = det_pq*(coef_qq*p_np1[a] - coef_pq*q_np1[a]); + q_np1[a] = det_pq*(coef_pp*q_np1[a] - coef_qp*aux_p_np1); + } + } + } ); + + /// synchronize pressure fields + FieldIdentifiers fieldsToBeSync; + fieldsToBeSync.addFields( FieldLocation::Node, { fields::wavesolverfields::Pressure_p_np1::key() } ); + fieldsToBeSync.addFields( FieldLocation::Node, { fields::wavesolverfields::Pressure_q_np1::key() } ); + + CommunicationTools & syncFields = CommunicationTools::getInstance(); + syncFields.synchronizeFields( fieldsToBeSync, + mesh, + domain.getNeighbors(), + true ); + + // compute the seismic traces since last step. + arrayView2d< real32 > const pReceivers = m_pressureNp1AtReceivers.toView(); + + computeAllSeismoTraces( time_n, dt, p_np1, p_n, pReceivers ); + + /// prepare next step + forAll< EXEC_POLICY >( nodeManager.size(), [=] GEOS_HOST_DEVICE ( localIndex const a ) + { + stiffnessVector_p[a] = 0.0; + stiffnessVector_q[a] = 0.0; + rhs[a] = 0.0; + } ); + + } ); + + return dt; +} + +void AcousticVTIWaveEquationSEM::cleanup( real64 const time_n, + integer const cycleNumber, + integer const eventCounter, + real64 const eventProgress, + DomainPartition & domain ) +{ + // call the base class cleanup (for reporting purposes) + SolverBase::cleanup( time_n, cycleNumber, eventCounter, eventProgress, domain ); + + // compute the remaining seismic traces, if needed + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & ) + { + NodeManager & nodeManager = mesh.getNodeManager(); + arrayView1d< real32 const > const p_n = nodeManager.getField< fields::wavesolverfields::Pressure_p_n >(); + arrayView1d< real32 const > const p_np1 = nodeManager.getField< fields::wavesolverfields::Pressure_p_np1 >(); + arrayView2d< real32 > const pReceivers = m_pressureNp1AtReceivers.toView(); + computeAllSeismoTraces( time_n, 0, p_np1, p_n, pReceivers ); + } ); +} + +void AcousticVTIWaveEquationSEM::computeAllSeismoTraces( real64 const time_n, + real64 const dt, + arrayView1d< real32 const > const var_np1, + arrayView1d< real32 const > const var_n, + arrayView2d< real32 > varAtReceivers ) +{ + + /* + * In forward case we compute seismo if time_n + dt is the first time + * step after the timeSeismo to write. + * + * time_n timeSeismo time_n + dt + * ---|--------------|-------------| + * + * In backward (time_n goes decreasing) case we compute seismo if + * time_n is the last time step before the timeSeismo to write. + * + * time_n - dt timeSeismo time_n + * ---|--------------|-------------| + */ + for( real64 timeSeismo; + (m_forward)?((timeSeismo = m_dtSeismoTrace*m_indexSeismoTrace) <= (time_n + dt + epsilonLoc) && m_indexSeismoTrace < m_nsamplesSeismoTrace): + ((timeSeismo = m_dtSeismoTrace*(m_nsamplesSeismoTrace-m_indexSeismoTrace-1)) >= (time_n - dt - epsilonLoc) && m_indexSeismoTrace < m_nsamplesSeismoTrace); + m_indexSeismoTrace++ ) + { + WaveSolverUtils::computeSeismoTrace( time_n, (m_forward)?dt:-dt, timeSeismo, (m_forward)?m_indexSeismoTrace:(m_nsamplesSeismoTrace-m_indexSeismoTrace-1), m_receiverNodeIds, m_receiverConstants, + m_receiverIsLocal, + m_nsamplesSeismoTrace, m_outputSeismoTrace, var_np1, var_n, varAtReceivers ); + } +} + + +REGISTER_CATALOG_ENTRY( SolverBase, AcousticVTIWaveEquationSEM, string const &, dataRepository::Group * const ) + +} /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/wavePropagation/AcousticVTIWaveEquationSEM.hpp b/src/coreComponents/physicsSolvers/wavePropagation/AcousticVTIWaveEquationSEM.hpp new file mode 100644 index 00000000000..1fe0ede62af --- /dev/null +++ b/src/coreComponents/physicsSolvers/wavePropagation/AcousticVTIWaveEquationSEM.hpp @@ -0,0 +1,163 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 TotalEnergies + * Copyright (c) 2019- GEOS Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + + +/** + * @file AcousticVTIWaveEquationSEM.hpp + */ + +#ifndef GEOS_PHYSICSSOLVERS_WAVEPROPAGATION_ACOUSTICVTIWAVEEQUATIONSEM_HPP_ +#define GEOS_PHYSICSSOLVERS_WAVEPROPAGATION_ACOUSTICVTIWAVEEQUATIONSEM_HPP_ + +#include "mesh/MeshFields.hpp" +#include "physicsSolvers/SolverBase.hpp" +#include "WaveSolverBase.hpp" +#include "WaveSolverBaseFields.hpp" + +namespace geos +{ + +class AcousticVTIWaveEquationSEM : public WaveSolverBase +{ +public: + + using EXEC_POLICY = parallelDevicePolicy< 32 >; + using ATOMIC_POLICY = AtomicPolicy< EXEC_POLICY >; + + AcousticVTIWaveEquationSEM( const std::string & name, + Group * const parent ); + + static string catalogName() { return "AcousticVTISEM"; } + + virtual void initializePreSubGroups() override; + + virtual void registerDataOnMesh( Group & meshBodies ) override final; + + /** + * @defgroup Solver Interface Functions + * + * These functions provide the primary interface that is required for derived classes + */ + /**@{*/ + virtual real64 explicitStepForward( real64 const & time_n, + real64 const & dt, + integer const cycleNumber, + DomainPartition & domain, + bool const computeGradient ) override; + + + virtual real64 explicitStepBackward( real64 const & GEOS_UNUSED_PARAM( time_n ), + real64 const & GEOS_UNUSED_PARAM( dt ), + integer const GEOS_UNUSED_PARAM( cycleNumber ), + DomainPartition & GEOS_UNUSED_PARAM( domain ), + bool const GEOS_UNUSED_PARAM( computeGradient ) ) override; + + /**@}*/ + + /** + * @brief Multiply the precomputed term by the Ricker and add to the right-hand side + * @param cycleNumber the cycle number/step number of evaluation of the source + * @param rhs the right hand side vector to be computed + */ + virtual void addSourceToRightHandSide( integer const & cycleNumber, arrayView1d< real32 > const rhs ); + + /** + * TODO: move implementation into WaveSolverBase + * @brief Computes the traces on all receivers (see @computeSeismoTraces) up to time_n+dt + * @param time_n the time corresponding to the field values pressure_n + * @param dt the simulation timestep + * @param var_np1 the field values at time_n + dt + * @param var_n the field values at time_n + * @param varAtReceivers the array holding the trace values, where the output is written + */ + virtual void computeAllSeismoTraces( real64 const time_n, + real64 const dt, + arrayView1d< real32 const > const var_np1, + arrayView1d< real32 const > const var_n, + arrayView2d< real32 > varAtReceivers ); + + /** + * @brief Overridden from ExecutableGroup. Used to write last seismogram if needed. + */ + virtual void cleanup( real64 const time_n, integer const cycleNumber, integer const eventCounter, real64 const eventProgress, DomainPartition & domain ) override; + + struct viewKeyStruct : WaveSolverBase::viewKeyStruct + { + static constexpr char const * pressureNp1AtReceiversString() { return "pressureNp1AtReceivers"; } + static constexpr char const * lateralSurfaceString() { return "LateralSurface"; } + static constexpr char const * bottomSurfaceString() { return "BottomSurface"; } + } waveEquationViewKeys; + + + /** internal function to the class to compute explicitStep either for backward or forward. + * (requires not to be private because it is called from GEOS_HOST_DEVICE method) + * @param time_n time at the beginning of the step + * @param dt the perscribed timestep + * @param cycleNumber the current cycle number + * @param domain the domain object + * @return return the timestep that was achieved during the step. + */ + real64 explicitStepInternal( real64 const & time_n, + real64 const & dt, + integer const cycleNumber, + DomainPartition & domain ); + + /** + * @brief (Empty but must be defined) Initialize Perfectly Matched Layer (PML) information + */ + virtual void initializePML() override; + + /** + * @brief (Empty but must be defined) Apply Perfectly Matched Layer (PML) to the regions defined in the geometry box from the xml + * @param time the time to apply the BC + * @param domain the partition domain + */ + virtual void applyPML( real64 const GEOS_UNUSED_PARAM( time ), DomainPartition & GEOS_UNUSED_PARAM( domain ) ) override; + +protected: + + virtual void postProcessInput() override final; + + virtual void initializePostInitialConditionsPreSubGroups() override final; + +private: + + /** + * @brief Locate sources and receivers position in the mesh elements, evaluate the basis functions at each point and save them to the + * corresponding elements nodes. + * @param mesh mesh of the computational domain + */ + virtual void precomputeSourceAndReceiverTerm( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; + + /** + * @brief Compute the lateral and bottom surface Field indicators of the boxed domain + * @param domain the partition domain + */ + virtual void precomputeSurfaceFieldIndicator( DomainPartition & domain ); + + /** + * @brief Apply free surface condition to the face define in the geometry box from the xml + * @param time the time to apply the BC + * @param domain the partition domain + */ + virtual void applyFreeSurfaceBC( real64 const time, DomainPartition & domain ) override; + + /// Pressure_p_np1 at the receiver location for each time step for each receiver + array2d< real32 > m_pressureNp1AtReceivers; + +}; + +} /* namespace geos */ + +#endif /* GEOS_PHYSICSSOLVERS_WAVEPROPAGATION_AcousticVTIWaveEquationSEM_HPP_ */ diff --git a/src/coreComponents/physicsSolvers/wavePropagation/AcousticVTIWaveEquationSEMKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/AcousticVTIWaveEquationSEMKernel.hpp new file mode 100644 index 00000000000..2972bfb3420 --- /dev/null +++ b/src/coreComponents/physicsSolvers/wavePropagation/AcousticVTIWaveEquationSEMKernel.hpp @@ -0,0 +1,567 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 TotalEnergies + * Copyright (c) 2019- GEOS Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file AcousticVTIWaveEquationSEMKernel.hpp + */ + +#ifndef GEOS_PHYSICSSOLVERS_WAVEPROPAGATION_ACOUSTICVTIWAVEEQUATIONSEMKERNEL_HPP_ +#define GEOS_PHYSICSSOLVERS_WAVEPROPAGATION_ACOUSTICVTIWAVEEQUATIONSEMKERNEL_HPP_ + +#include "finiteElement/elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp" +#include "finiteElement/kernelInterface/KernelBase.hpp" +#include "WaveSolverBaseFields.hpp" +#include "WaveSolverUtils.hpp" + + +namespace geos +{ + +/// Namespace to contain the acoustic wave kernels. +namespace acousticVTIWaveEquationSEMKernels +{ + +struct PrecomputeSourceAndReceiverKernel +{ + + + + /** + * @brief Launches the precomputation of the source and receiver terms + * @tparam EXEC_POLICY execution policy + * @tparam FE_TYPE finite element type + * @param[in] size the number of cells in the subRegion + * @param[in] numNodesPerElem number of nodes per element + * @param[in] numFacesPerElem number of faces per element + * @param[in] X coordinates of the nodes + * @param[in] elemGhostRank the ghost ranks + * @param[in] elemsToNodes map from element to nodes + * @param[in] elemsToFaces map from element to faces + * @param[in] facesToNodes map from faces to nodes + * @param[in] elemCenter coordinates of the element centers + * @param[in] sourceCoordinates coordinates of the source terms + * @param[out] sourceIsAccessible flag indicating whether the source is accessible or not + * @param[out] sourceNodeIds indices of the nodes of the element where the source is located + * @param[out] sourceNodeConstants constant part of the source terms + * @param[in] receiverCoordinates coordinates of the receiver terms + * @param[out] receiverIsLocal flag indicating whether the receiver is local or not + * @param[out] receiverNodeIds indices of the nodes of the element where the receiver is located + * @param[out] receiverNodeConstants constant part of the receiver term + * @param[out] sourceValue the value of the source + * @param[in] dt the time step size + * @param[in] timeSourceFrequency the time frequency of the source + * @param[in] rickerOrder the order of the ricker + */ + template< typename EXEC_POLICY, typename FE_TYPE > + static void + launch( localIndex const size, + localIndex const numNodesPerElem, + localIndex const numFacesPerElem, + arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const X, + arrayView1d< integer const > const elemGhostRank, + arrayView2d< localIndex const, cells::NODE_MAP_USD > const & elemsToNodes, + arrayView2d< localIndex const > const elemsToFaces, + arrayView2d< real64 const > const & elemCenter, + arrayView2d< real64 const > const faceNormal, + arrayView2d< real64 const > const faceCenter, + arrayView2d< real64 const > const sourceCoordinates, + arrayView1d< localIndex > const sourceIsAccessible, + arrayView2d< localIndex > const sourceNodeIds, + arrayView2d< real64 > const sourceConstants, + arrayView2d< real64 const > const receiverCoordinates, + arrayView1d< localIndex > const receiverIsLocal, + arrayView2d< localIndex > const receiverNodeIds, + arrayView2d< real64 > const receiverConstants, + arrayView2d< real32 > const sourceValue, + real64 const dt, + real32 const timeSourceFrequency, + localIndex const rickerOrder ) + { + + forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + real64 const center[3] = { elemCenter[k][0], + elemCenter[k][1], + elemCenter[k][2] }; + + // Step 1: locate the sources, and precompute the source term + + /// loop over all the source that haven't been found yet + for( localIndex isrc = 0; isrc < sourceCoordinates.size( 0 ); ++isrc ) + { + if( sourceIsAccessible[isrc] == 0 ) + { + real64 const coords[3] = { sourceCoordinates[isrc][0], + sourceCoordinates[isrc][1], + sourceCoordinates[isrc][2] }; + + bool const sourceFound = + WaveSolverUtils::locateSourceElement( numFacesPerElem, + center, + faceNormal, + faceCenter, + elemsToFaces[k], + coords ); + if( sourceFound ) + { + real64 coordsOnRefElem[3]{}; + + + WaveSolverUtils::computeCoordinatesOnReferenceElement< FE_TYPE >( coords, + elemsToNodes[k], + X, + coordsOnRefElem ); + + sourceIsAccessible[isrc] = 1; + real64 Ntest[FE_TYPE::numNodes]; + FE_TYPE::calcN( coordsOnRefElem, Ntest ); + + for( localIndex a = 0; a < numNodesPerElem; ++a ) + { + sourceNodeIds[isrc][a] = elemsToNodes[k][a]; + sourceConstants[isrc][a] = Ntest[a]; + } + + for( localIndex cycle = 0; cycle < sourceValue.size( 0 ); ++cycle ) + { + real64 const time = cycle*dt; + sourceValue[cycle][isrc] = WaveSolverUtils::evaluateRicker( time, timeSourceFrequency, rickerOrder ); + } + } + } + } // end loop over all sources + + + // Step 2: locate the receivers, and precompute the receiver term + + /// loop over all the receivers that haven't been found yet + for( localIndex ircv = 0; ircv < receiverCoordinates.size( 0 ); ++ircv ) + { + if( receiverIsLocal[ircv] == 0 ) + { + real64 const coords[3] = { receiverCoordinates[ircv][0], + receiverCoordinates[ircv][1], + receiverCoordinates[ircv][2] }; + + real64 coordsOnRefElem[3]{}; + bool const receiverFound = + WaveSolverUtils::locateSourceElement( numFacesPerElem, + center, + faceNormal, + faceCenter, + elemsToFaces[k], + coords ); + + if( receiverFound && elemGhostRank[k] < 0 ) + { + WaveSolverUtils::computeCoordinatesOnReferenceElement< FE_TYPE >( coords, + elemsToNodes[k], + X, + coordsOnRefElem ); + + receiverIsLocal[ircv] = 1; + + real64 Ntest[FE_TYPE::numNodes]; + FE_TYPE::calcN( coordsOnRefElem, Ntest ); + + for( localIndex a = 0; a < numNodesPerElem; ++a ) + { + receiverNodeIds[ircv][a] = elemsToNodes[k][a]; + receiverConstants[ircv][a] = Ntest[a]; + } + } + } + } // end loop over receivers + + } ); + + } +}; + +template< typename FE_TYPE > +struct MassMatrixKernel +{ + + MassMatrixKernel( FE_TYPE const & finiteElement ) + : m_finiteElement( finiteElement ) + {} + + /** + * @brief Launches the precomputation of the mass matrices + * @tparam EXEC_POLICY the execution policy + * @tparam ATOMIC_POLICY the atomic policy + * @param[in] size the number of cells in the subRegion + * @param[in] X coordinates of the nodes + * @param[in] elemsToNodes map from element to nodes + * @param[in] velocity cell-wise velocity + * @param[out] mass diagonal of the mass matrix + */ + template< typename EXEC_POLICY, typename ATOMIC_POLICY > + void + launch( localIndex const size, + arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const X, + arrayView2d< localIndex const, cells::NODE_MAP_USD > const elemsToNodes, + arrayView1d< real32 const > const velocity, + arrayView1d< real32 > const mass ) + + { + forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + + constexpr localIndex numNodesPerElem = FE_TYPE::numNodes; + constexpr localIndex numQuadraturePointsPerElem = FE_TYPE::numQuadraturePoints; + + real32 const invC2 = 1.0 / ( velocity[k] * velocity[k] ); + real64 xLocal[ numNodesPerElem ][ 3 ]; + for( localIndex a = 0; a < numNodesPerElem; ++a ) + { + for( localIndex i = 0; i < 3; ++i ) + { + xLocal[a][i] = X( elemsToNodes( k, a ), i ); + } + } + + for( localIndex q = 0; q < numQuadraturePointsPerElem; ++q ) + { + real32 const localIncrement = invC2 * m_finiteElement.computeMassTerm( q, xLocal ); + RAJA::atomicAdd< ATOMIC_POLICY >( &mass[elemsToNodes[k][q]], localIncrement ); + } + } ); // end loop over element + } + + /// The finite element space/discretization object for the element type in the subRegion + FE_TYPE const & m_finiteElement; + +}; + +template< typename FE_TYPE > +struct DampingMatrixKernel +{ + + DampingMatrixKernel( FE_TYPE const & finiteElement ) + : m_finiteElement( finiteElement ) + {} + + /** + * @brief Launches the precomputation of the damping matrices + * @tparam EXEC_POLICY the execution policy + * @tparam ATOMIC_POLICY the atomic policy + * @param[in] size the number of cells in the subRegion + * @param[in] X coordinates of the nodes + * @param[in] facesToElems map from faces to elements + * @param[in] facesToNodes map from face to nodes + * @param[in] facesDomainBoundaryIndicator flag equal to 1 if the face is on the boundary, and to 0 otherwise + * @param[in] freeSurfaceFaceIndicator flag equal to 1 if the face is on the free surface, and to 0 otherwise + * @param[in] lateralSurfaceFaceIndicator flag equal to 1 if the face is on a lateral surface, and to 0 otherwise + * @param[in] bottomSurfaceFaceIndicator flag equal to 1 if the face is on the bottom surface, and to 0 otherwise + * @param[in] velocity cell-wise velocity + * @param[in] epsilon cell-wise Thomsen epsilon parameter + * @param[in] delta cell-wise Thomsen delta parameter + * @param[in] vti_f cell-wise f parameter in Fletcher's equation + * @param[out] damping_p diagonal of the damping matrix for quantities in p in p-equation + * @param[out] damping_q diagonal of the damping matrix for quantities in q in q-equation + * @param[out] damping_pq diagonal of the damping matrix for quantities in q in p-equation + * @param[out] damping_qp diagonal of the damping matrix for quantities in p in q-equation + */ + template< typename EXEC_POLICY, typename ATOMIC_POLICY > + void + launch( localIndex const size, + arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const X, + arrayView2d< localIndex const > const facesToElems, + ArrayOfArraysView< localIndex const > const facesToNodes, + arrayView1d< integer const > const facesDomainBoundaryIndicator, + arrayView1d< localIndex const > const freeSurfaceFaceIndicator, + arrayView1d< localIndex const > const lateralSurfaceFaceIndicator, + arrayView1d< localIndex const > const bottomSurfaceFaceIndicator, + arrayView1d< real32 const > const velocity, + arrayView1d< real32 const > const epsilon, + arrayView1d< real32 const > const delta, + arrayView1d< real32 const > const vti_f, + arrayView1d< real32 > const damping_p, + arrayView1d< real32 > const damping_q, + arrayView1d< real32 > const damping_pq, + arrayView1d< real32 > const damping_qp ) + { + forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const f ) + { + // face on the domain boundary and not on free surface + if( facesDomainBoundaryIndicator[f] == 1 && freeSurfaceFaceIndicator[f] != 1 ) + { + localIndex k = facesToElems( f, 0 ); + if( k < 0 ) + { + k = facesToElems( f, 1 ); + } + + // ABC coefficients + real32 alpha = 1.0 / velocity[k]; + // VTI coefficients + real32 vti_p_xy= 0, vti_p_z = 0, vti_pq_z = 0; + real32 vti_q_xy= 0, vti_q_z = 0, vti_qp_xy= 0; + if( lateralSurfaceFaceIndicator[f] == 1 ) + { + // ABC coefficients updated to fit horizontal velocity + alpha /= sqrt( 1+2*epsilon[k] ); + + vti_p_xy = (1+2*epsilon[k]); + vti_q_xy = -(vti_f[k]-1); + vti_qp_xy = (vti_f[k]+2*delta[k]); + } + if( bottomSurfaceFaceIndicator[f] == 1 ) + { + // ABC coefficients updated to fit horizontal velocity + alpha /= sqrt( 1+2*delta[k] ); + vti_p_z = -(vti_f[k]-1); + vti_pq_z = vti_f[k]; + vti_q_z = 1; + } + + constexpr localIndex numNodesPerFace = FE_TYPE::numNodesPerFace; + + real64 xLocal[ numNodesPerFace ][ 3 ]; + for( localIndex a = 0; a < numNodesPerFace; ++a ) + { + for( localIndex i = 0; i < 3; ++i ) + { + xLocal[a][i] = X( facesToNodes( k, a ), i ); + } + } + + for( localIndex q = 0; q < numNodesPerFace; ++q ) + { + real32 const localIncrement_p = alpha*(vti_p_xy + vti_p_z) * m_finiteElement.computeDampingTerm( q, xLocal ); + RAJA::atomicAdd< ATOMIC_POLICY >( &damping_p[facesToNodes[f][q]], localIncrement_p ); + + real32 const localIncrement_pq = alpha*vti_pq_z * m_finiteElement.computeDampingTerm( q, xLocal ); + RAJA::atomicAdd< ATOMIC_POLICY >( &damping_pq[facesToNodes[f][q]], localIncrement_pq ); + + real32 const localIncrement_q = alpha*(vti_q_xy + vti_q_z) * m_finiteElement.computeDampingTerm( q, xLocal ); + RAJA::atomicAdd< ATOMIC_POLICY >( &damping_q[facesToNodes[f][q]], localIncrement_q ); + + real32 const localIncrement_qp = alpha*vti_qp_xy * m_finiteElement.computeDampingTerm( q, xLocal ); + RAJA::atomicAdd< ATOMIC_POLICY >( &damping_qp[facesToNodes[f][q]], localIncrement_qp ); + } + } + } ); // end loop over element + } + + /// The finite element space/discretization object for the element type in the subRegion + FE_TYPE const & m_finiteElement; + +}; + + + +/** + * @brief Implements kernels for solving the acoustic wave equations + * explicit central FD method and SEM + * @copydoc geos::finiteElement::KernelBase + * @tparam SUBREGION_TYPE The type of subregion that the kernel will act on. + * + * ### AcousticVTIWaveEquationSEMKernel Description + * Implements the KernelBase interface functions required for solving + * the acoustic wave equations using the + * "finite element kernel application" functions such as + * geos::finiteElement::RegionBasedKernelApplication. + * + * The number of degrees of freedom per support point for both + * the test and trial spaces are specified as `1`. + */ + + +template< typename SUBREGION_TYPE, + typename CONSTITUTIVE_TYPE, + typename FE_TYPE > +class ExplicitAcousticVTISEM : public finiteElement::KernelBase< SUBREGION_TYPE, + CONSTITUTIVE_TYPE, + FE_TYPE, + 1, + 1 > +{ +public: + + /// Alias for the base class; + using Base = finiteElement::KernelBase< SUBREGION_TYPE, + CONSTITUTIVE_TYPE, + FE_TYPE, + 1, + 1 >; + + /// Maximum number of nodes per element, which is equal to the maxNumTestSupportPointPerElem and + /// maxNumTrialSupportPointPerElem by definition. When the FE_TYPE is not a Virtual Element, this + /// will be the actual number of nodes per element. + static constexpr int numNodesPerElem = Base::maxNumTestSupportPointsPerElem; + + using Base::numDofPerTestSupportPoint; + using Base::numDofPerTrialSupportPoint; + using Base::m_elemsToNodes; + using Base::m_elemGhostRank; + using Base::m_constitutiveUpdate; + using Base::m_finiteElementSpace; + +//***************************************************************************** + /** + * @brief Constructor + * @copydoc geos::finiteElement::KernelBase::KernelBase + * @param nodeManager Reference to the NodeManager object. + * @param edgeManager Reference to the EdgeManager object. + * @param faceManager Reference to the FaceManager object. + * @param targetRegionIndex Index of the region the subregion belongs to. + * @param dt The time interval for the step. + * elements to be processed during this kernel launch. + */ + ExplicitAcousticVTISEM( NodeManager & nodeManager, + EdgeManager const & edgeManager, + FaceManager const & faceManager, + localIndex const targetRegionIndex, + SUBREGION_TYPE const & elementSubRegion, + FE_TYPE const & finiteElementSpace, + CONSTITUTIVE_TYPE & inputConstitutiveType, + real64 const dt ): + Base( elementSubRegion, + finiteElementSpace, + inputConstitutiveType ), + m_X( nodeManager.getField< fields::referencePosition32 >() ), + m_p_n( nodeManager.getField< fields::wavesolverfields::Pressure_p_n >() ), + m_q_n( nodeManager.getField< fields::wavesolverfields::Pressure_q_n >() ), + m_stiffnessVector_p( nodeManager.getField< fields::wavesolverfields::StiffnessVector_p >() ), + m_stiffnessVector_q( nodeManager.getField< fields::wavesolverfields::StiffnessVector_q >() ), + m_epsilon( elementSubRegion.template getField< fields::wavesolverfields::Epsilon >() ), + m_delta( elementSubRegion.template getField< fields::wavesolverfields::Delta >() ), + m_vti_f( elementSubRegion.template getField< fields::wavesolverfields::F >() ), + m_dt( dt ) + { + GEOS_UNUSED_VAR( edgeManager ); + GEOS_UNUSED_VAR( faceManager ); + GEOS_UNUSED_VAR( targetRegionIndex ); + } + + //***************************************************************************** + /** + * @copydoc geos::finiteElement::KernelBase::StackVariables + * + * ### ExplicitAcousticVTISEM Description + * Adds a stack arrays for the nodal force, primary displacement variable, etc. + */ + struct StackVariables : Base::StackVariables + { +public: + GEOS_HOST_DEVICE + StackVariables(): + xLocal() + {} + + /// C-array stack storage for element local the nodal positions. + real64 xLocal[ numNodesPerElem ][ 3 ]; + }; + //*************************************************************************** + + + /** + * @copydoc geos::finiteElement::KernelBase::setup + * + * Copies the primary variable, and position into the local stack array. + */ + GEOS_HOST_DEVICE + GEOS_FORCE_INLINE + void setup( localIndex const k, + StackVariables & stack ) const + { + /// numDofPerTrialSupportPoint = 1 + for( localIndex a=0; a< numNodesPerElem; ++a ) + { + localIndex const nodeIndex = m_elemsToNodes( k, a ); + for( int i=0; i< 3; ++i ) + { + stack.xLocal[ a ][ i ] = m_X[ nodeIndex ][ i ]; + } + } + } + + /** + * @copydoc geos::finiteElement::KernelBase::quadraturePointKernel + * + * ### ExplicitAcousticVTISEM Description + * Calculates stiffness vector + * + */ + GEOS_HOST_DEVICE + GEOS_FORCE_INLINE + void quadraturePointKernel( localIndex const k, + localIndex const q, + StackVariables & stack ) const + { + // Pseudo Stiffness xy + m_finiteElementSpace.template computeStiffnessxyTerm( q, stack.xLocal, [&] ( int i, int j, real64 val ) + { + real32 const localIncrement_p = val*(-1-2*m_epsilon[k])*m_p_n[m_elemsToNodes[k][j]]; + RAJA::atomicAdd< parallelDeviceAtomic >( &m_stiffnessVector_p[m_elemsToNodes[k][i]], localIncrement_p ); + real32 const localIncrement_q = val*((-2*m_delta[k]-m_vti_f[k])*m_p_n[m_elemsToNodes[k][j]] +(m_vti_f[k]-1)*m_q_n[m_elemsToNodes[k][j]]); + RAJA::atomicAdd< parallelDeviceAtomic >( &m_stiffnessVector_q[m_elemsToNodes[k][i]], localIncrement_q ); + } ); + + // Pseudo-Stiffness z + + m_finiteElementSpace.template computeStiffnesszTerm( q, stack.xLocal, [&] ( int i, int j, real64 val ) + { + real32 const localIncrement_p = val*((m_vti_f[k]-1)*m_p_n[m_elemsToNodes[k][j]] - m_vti_f[k]*m_q_n[m_elemsToNodes[k][j]]); + RAJA::atomicAdd< parallelDeviceAtomic >( &m_stiffnessVector_p[m_elemsToNodes[k][i]], localIncrement_p ); + + real32 const localIncrement_q = -val*m_q_n[m_elemsToNodes[k][j]]; + RAJA::atomicAdd< parallelDeviceAtomic >( &m_stiffnessVector_q[m_elemsToNodes[k][i]], localIncrement_q ); + } ); + } + +protected: + /// The array containing the nodal position array. + arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const m_X; + + /// The array containing the nodal pressure array. + arrayView1d< real32 const > const m_p_n; + + /// The array containing the nodal auxiliary variable array. + arrayView1d< real32 const > const m_q_n; + + /// The array containing the product of the stiffness matrix and the nodal pressure for the equation in p. + arrayView1d< real32 > const m_stiffnessVector_p; + + /// The array containing the product of the stiffness matrix and the nodal pressure for the equation in q. + arrayView1d< real32 > const m_stiffnessVector_q; + + /// The array containing the epsilon Thomsen parameter. + arrayView1d< real32 const > const m_epsilon; + + /// The array containing the delta Thomsen parameter. + arrayView1d< real32 const > const m_delta; + + /// The array containing the f parameter. + arrayView1d< real32 const > const m_vti_f; + + /// The time increment for this time integration step. + real64 const m_dt; + + +}; + + + +/// The factory used to construct a ExplicitAcousticWaveEquation kernel. +using ExplicitAcousticVTISEMFactory = finiteElement::KernelFactory< ExplicitAcousticVTISEM, + real64 >; + +} // namespace acousticVTIWaveEquationSEMKernels + +} // namespace geos + +#endif //GEOS_PHYSICSSOLVERS_WAVEPROPAGATION_ACOUSTICVTIWAVEEQUATIONSEMKERNEL_HPP_ diff --git a/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEM.cpp index f8fef564ce3..414d35729d9 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEM.cpp @@ -260,10 +260,11 @@ void AcousticWaveEquationSEM::initializePostInitialConditionsPreSubGroups() NodeManager & nodeManager = mesh.getNodeManager(); FaceManager & faceManager = mesh.getFaceManager(); + ElementRegionManager & elemManager = mesh.getElemManager(); /// get the array of indicators: 1 if the face is on the boundary; 0 otherwise arrayView1d< integer const > const & facesDomainBoundaryIndicator = faceManager.getDomainBoundaryIndicator(); - arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const X32 = nodeManager.getField< fields::referencePosition32 >().toViewConst(); + arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords = nodeManager.getField< fields::referencePosition32 >().toViewConst(); /// get face to node map ArrayOfArraysView< localIndex const > const facesToNodes = faceManager.nodeList().toViewConst(); @@ -283,28 +284,29 @@ void AcousticWaveEquationSEM::initializePostInitialConditionsPreSubGroups() /// get array of indicators: 1 if face is on the free surface; 0 otherwise arrayView1d< localIndex const > const freeSurfaceFaceIndicator = faceManager.getField< fields::FreeSurfaceFaceIndicator >(); - mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, - CellElementSubRegion & elementSubRegion ) + elemManager.forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, + CellElementSubRegion & elementSubRegion ) { + finiteElement::FiniteElementBase const & + fe = elementSubRegion.getReference< finiteElement::FiniteElementBase >( getDiscretizationName() ); + arrayView2d< localIndex const, cells::NODE_MAP_USD > const elemsToNodes = elementSubRegion.nodeList(); - arrayView2d< localIndex const > const facesToElements = faceManager.elementList(); + arrayView2d< localIndex const > const elemsToFaces = elementSubRegion.faceList(); + arrayView1d< real32 const > const velocity = elementSubRegion.getField< fields::MediumVelocity >(); arrayView1d< real32 const > const density = elementSubRegion.getField< fields::MediumDensity >(); /// Partial gradient if gradient as to be computed arrayView1d< real32 > grad = elementSubRegion.getField< fields::PartialGradient >(); - { - grad.zero(); - } - finiteElement::FiniteElementBase const & - fe = elementSubRegion.getReference< finiteElement::FiniteElementBase >( getDiscretizationName() ); + grad.zero(); + finiteElement::FiniteElementDispatchHandler< SEM_FE_TYPES >::dispatch3D( fe, [&] ( auto const finiteElement ) { using FE_TYPE = TYPEOFREF( finiteElement ); acousticWaveEquationSEMKernels::MassMatrixKernel< FE_TYPE > kernelM( finiteElement ); kernelM.template launch< EXEC_POLICY, ATOMIC_POLICY >( elementSubRegion.size(), - X32, + nodeCoords, elemsToNodes, velocity, density, @@ -312,10 +314,9 @@ void AcousticWaveEquationSEM::initializePostInitialConditionsPreSubGroups() { GEOS_MARK_SCOPE( DampingMatrixKernel ); acousticWaveEquationSEMKernels::DampingMatrixKernel< FE_TYPE > kernelD( finiteElement ); - - kernelD.template launch< EXEC_POLICY, ATOMIC_POLICY >( faceManager.size(), - X32, - facesToElements, + kernelD.template launch< EXEC_POLICY, ATOMIC_POLICY >( elementSubRegion.size(), + nodeCoords, + elemsToFaces, facesToNodes, facesDomainBoundaryIndicator, freeSurfaceFaceIndicator, @@ -747,18 +748,6 @@ void AcousticWaveEquationSEM::applyPML( real64 const time, DomainPartition & dom } -/** - * Checks if a directory exists. - * - * @param dirName Directory name to check existence of. - * @return true is dirName exists and is a directory. - */ -bool dirExists( const std::string & dirName ) -{ - struct stat buffer; - return stat( dirName.c_str(), &buffer ) == 0; -} - real64 AcousticWaveEquationSEM::explicitStepForward( real64 const & time_n, real64 const & dt, integer cycleNumber, @@ -813,8 +802,10 @@ real64 AcousticWaveEquationSEM::explicitStepForward( real64 const & time_n, std::string fileName = GEOS_FMT( "lifo/rank_{:05}/pressuredt2_{: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 && !dirExists( dirName )) + if( string::npos != (size_t)lastDirSeparator && !directoryExists( dirName )) + { makeDirsForPath( dirName ); + } //std::string fileName = GEOS_FMT( "pressuredt2_{:06}_{:08}_{:04}.dat", m_shotIndex, cycleNumber, rank ); //const int fileDesc = open( fileName.c_str(), O_CREAT | O_WRONLY | O_DIRECT | O_TRUNC, S_IRUSR | S_IWUSR | S_IRGRP | S_IWGRP | diff --git a/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEMKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEMKernel.hpp index c0a9044f7b6..b9815ce1b94 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEMKernel.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEMKernel.hpp @@ -211,25 +211,25 @@ struct MassMatrixKernel arrayView1d< real32 > const mass ) { - forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k ) + forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const e ) { constexpr localIndex numNodesPerElem = FE_TYPE::numNodes; constexpr localIndex numQuadraturePointsPerElem = FE_TYPE::numQuadraturePoints; - real32 const invC2 = 1.0 / ( density[k] * velocity[k] * velocity[k] ); + real32 const invC2 = 1.0 / ( density[e] * velocity[e] * velocity[e] ); real64 xLocal[ numNodesPerElem ][ 3 ]; for( localIndex a = 0; a < numNodesPerElem; ++a ) { for( localIndex i = 0; i < 3; ++i ) { - xLocal[a][i] = X( elemsToNodes( k, a ), i ); + xLocal[a][i] = X( elemsToNodes( e, a ), i ); } } for( localIndex q = 0; q < numQuadraturePointsPerElem; ++q ) { real32 const localIncrement = invC2 * m_finiteElement.computeMassTerm( q, xLocal ); - RAJA::atomicAdd< ATOMIC_POLICY >( &mass[elemsToNodes[k][q]], localIncrement ); + RAJA::atomicAdd< ATOMIC_POLICY >( &mass[elemsToNodes( e, q )], localIncrement ); } } ); // end loop over element } @@ -252,8 +252,8 @@ struct DampingMatrixKernel * @tparam EXEC_POLICY the execution policy * @tparam ATOMIC_POLICY the atomic policy * @param[in] size the number of cells in the subRegion - * @param[in] X coordinates of the nodes - * @param[in] facesToElems map from faces to elements + * @param[in] nodeCoords coordinates of the nodes + * @param[in] elemsToFaces map from elements to faces * @param[in] facesToNodes map from face to nodes * @param[in] facesDomainBoundaryIndicator flag equal to 1 if the face is on the boundary, and to 0 otherwise * @param[in] freeSurfaceFaceIndicator flag equal to 1 if the face is on the free surface, and to 0 otherwise @@ -264,8 +264,8 @@ struct DampingMatrixKernel template< typename EXEC_POLICY, typename ATOMIC_POLICY > void launch( localIndex const size, - arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const X, - arrayView2d< localIndex const > const facesToElems, + arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords, + arrayView2d< localIndex const > const elemsToFaces, ArrayOfArraysView< localIndex const > const facesToNodes, arrayView1d< integer const > const facesDomainBoundaryIndicator, arrayView1d< localIndex const > const freeSurfaceFaceIndicator, @@ -273,37 +273,33 @@ struct DampingMatrixKernel arrayView1d< real32 const > const density, arrayView1d< real32 > const damping ) { - forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const f ) + forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const e ) { - // face on the domain boundary and not on free surface - if( facesDomainBoundaryIndicator[f] == 1 && freeSurfaceFaceIndicator[f] != 1 ) + for( localIndex i = 0; i < elemsToFaces.size( 1 ); ++i ) { - localIndex k = facesToElems( f, 0 ); - if( k < 0 ) - { - k = facesToElems( f, 1 ); - } - - real32 const alpha = 1.0 / (density[k] * velocity[k]); - - constexpr localIndex numNodesPerFace = FE_TYPE::numNodesPerFace; - - real64 xLocal[ numNodesPerFace ][ 3 ]; - for( localIndex a = 0; a < numNodesPerFace; ++a ) + localIndex const f = elemsToFaces( e, i ); + // face on the domain boundary and not on free surface + if( facesDomainBoundaryIndicator[f] == 1 && freeSurfaceFaceIndicator[f] != 1 ) { - for( localIndex i = 0; i < 3; ++i ) + constexpr localIndex numNodesPerFace = FE_TYPE::numNodesPerFace; + real64 xLocal[ numNodesPerFace ][ 3 ]; + for( localIndex a = 0; a < numNodesPerFace; ++a ) { - xLocal[a][i] = X( facesToNodes( f, a ), i ); + for( localIndex d = 0; d < 3; ++d ) + { + xLocal[a][d] = nodeCoords( facesToNodes( f, a ), d ); + } } - } + real32 const alpha = 1.0 / (density[e] * velocity[e]); - for( localIndex q = 0; q < numNodesPerFace; ++q ) - { - real32 const localIncrement = alpha * m_finiteElement.computeDampingTerm( q, xLocal ); - RAJA::atomicAdd< ATOMIC_POLICY >( &damping[facesToNodes[f][q]], localIncrement ); + for( localIndex q = 0; q < numNodesPerFace; ++q ) + { + real32 const localIncrement = alpha * m_finiteElement.computeDampingTerm( q, xLocal ); + RAJA::atomicAdd< ATOMIC_POLICY >( &damping[facesToNodes( f, q )], localIncrement ); + } } } - } ); // end loop over element + } ); } /// The finite element space/discretization object for the element type in the subRegion diff --git a/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEM.cpp index e423f910c82..9a724f22f93 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEM.cpp @@ -383,7 +383,7 @@ void ElasticFirstOrderWaveEquationSEM::initializePostInitialConditionsPreSubGrou { arrayView2d< localIndex const, cells::NODE_MAP_USD > const elemsToNodes = elementSubRegion.nodeList(); - arrayView2d< localIndex const > const facesToElements = faceManager.elementList(); + arrayView2d< localIndex const > const elemsToFaces = elementSubRegion.faceList(); arrayView1d< real32 > const density = elementSubRegion.getField< wavesolverfields::MediumDensity >(); arrayView1d< real32 > const velocityVp = elementSubRegion.getField< wavesolverfields::MediumVelocityVp >(); arrayView1d< real32 > const velocityVs = elementSubRegion.getField< wavesolverfields::MediumVelocityVs >(); @@ -407,9 +407,9 @@ void ElasticFirstOrderWaveEquationSEM::initializePostInitialConditionsPreSubGrou elasticFirstOrderWaveEquationSEMKernels::DampingMatrixKernel< FE_TYPE > kernelD( finiteElement ); - kernelD.template launch< EXEC_POLICY, ATOMIC_POLICY >( faceManager.size(), + kernelD.template launch< EXEC_POLICY, ATOMIC_POLICY >( elementSubRegion.size(), X, - facesToElements, + elemsToFaces, facesToNodes, facesDomainBoundaryIndicator, freeSurfaceFaceIndicator, @@ -451,7 +451,7 @@ void ElasticFirstOrderWaveEquationSEM::applyFreeSurfaceBC( real64 const time, Do fsManager.apply( time, domain.getMeshBody( 0 ).getMeshLevel( m_discretizationName ), - string( "FreeSurface" ), + WaveSolverBase::viewKeyStruct::freeSurfaceString(), [&]( FieldSpecificationBase const & bc, string const &, SortedArrayView< localIndex const > const & targetSet, diff --git a/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEMKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEMKernel.hpp index c4564bbf791..973ca5c8403 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEMKernel.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEMKernel.hpp @@ -262,8 +262,8 @@ struct DampingMatrixKernel * @tparam EXEC_POLICY the execution policy * @tparam ATOMIC_POLICY the atomic policy * @param[in] size the number of cells in the subRegion - * @param[in] X coordinates of the nodes - * @param[in] facesToElems map from faces to elements + * @param[in] nodeCoords coordinates of the nodes + * @param[in] elemsToFaces map from elements to faces * @param[in] facesToNodes map from face to nodes * @param[in] facesDomainBoundaryIndicator flag equal to 1 if the face is on the boundary, and to 0 otherwise * @param[in] freeSurfaceFaceIndicator flag equal to 1 if the face is on the free surface, and to 0 otherwise @@ -273,8 +273,8 @@ struct DampingMatrixKernel template< typename EXEC_POLICY, typename ATOMIC_POLICY > void launch( localIndex const size, - arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const X, - arrayView2d< localIndex const > const facesToElems, + arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords, + arrayView2d< localIndex const > const elemsToFaces, ArrayOfArraysView< localIndex const > const facesToNodes, arrayView1d< integer const > const facesDomainBoundaryIndicator, arrayView1d< localIndex const > const freeSurfaceFaceIndicator, @@ -286,49 +286,39 @@ struct DampingMatrixKernel arrayView1d< real32 > const dampingy, arrayView1d< real32 > const dampingz ) { - forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const f ) + forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const e ) { - // face on the domain boundary and not on free surface - if( facesDomainBoundaryIndicator[f] == 1 && freeSurfaceFaceIndicator[f] != 1 ) + for( localIndex i = 0; i < elemsToFaces.size( 1 ); ++i ) { - localIndex k = facesToElems( f, 0 ); - if( k < 0 ) - { - k = facesToElems( f, 1 ); - } - - constexpr localIndex numNodesPerFace = FE_TYPE::numNodesPerFace; - - real64 xLocal[ numNodesPerFace ][ 3 ]; - for( localIndex a = 0; a < numNodesPerFace; ++a ) + localIndex const f = elemsToFaces( e, i ); + // face on the domain boundary and not on free surface + if( facesDomainBoundaryIndicator[f] == 1 && freeSurfaceFaceIndicator[f] != 1 ) { - for( localIndex i = 0; i < 3; ++i ) + constexpr localIndex numNodesPerFace = FE_TYPE::numNodesPerFace; + real64 xLocal[ numNodesPerFace ][ 3 ]; + for( localIndex a = 0; a < numNodesPerFace; ++a ) { - xLocal[a][i] = X( facesToNodes( f, a ), i ); + for( localIndex d = 0; d < 3; ++d ) + { + xLocal[a][d] = nodeCoords( facesToNodes( f, a ), d ); + } } - } - for( localIndex q = 0; q < numNodesPerFace; ++q ) - { - real32 const localIncrementx = density[k] * (velocityVp[k]*LvArray::math::abs( faceNormal[f][0] ) + velocityVs[k]*sqrt( faceNormal[f][1]*faceNormal[f][1] + - faceNormal[f][2]*faceNormal[f][2] ) )* m_finiteElement.computeDampingTerm( - q, - xLocal ); - real32 const localIncrementy = density[k] * (velocityVp[k]*LvArray::math::abs( faceNormal[f][1] ) + velocityVs[k]*sqrt( faceNormal[f][0]*faceNormal[f][0] + - faceNormal[f][2]*faceNormal[f][2] ) )* m_finiteElement.computeDampingTerm( - q, - xLocal ); - real32 const localIncrementz = density[k] * (velocityVp[k]*LvArray::math::abs( faceNormal[f][2] ) + velocityVs[k]*sqrt( faceNormal[f][0]*faceNormal[f][0] + - faceNormal[f][1]*faceNormal[f][1] ) )* m_finiteElement.computeDampingTerm( - q, - xLocal ); - - RAJA::atomicAdd< ATOMIC_POLICY >( &dampingx[facesToNodes[f][q]], localIncrementx ); - RAJA::atomicAdd< ATOMIC_POLICY >( &dampingy[facesToNodes[f][q]], localIncrementy ); - RAJA::atomicAdd< ATOMIC_POLICY >( &dampingz[facesToNodes[f][q]], localIncrementz ); + real32 const nx = faceNormal( f, 0 ), ny = faceNormal( f, 1 ), nz = faceNormal( f, 2 ); + for( localIndex q = 0; q < numNodesPerFace; ++q ) + { + real32 const aux = density[e] * m_finiteElement.computeDampingTerm( q, xLocal ); + real32 const localIncrementx = density[e] * (velocityVp[e] * abs( nx ) + velocityVs[e] * sqrt( pow( ny, 2 ) + pow( nz, 2 ) ) ) * aux; + real32 const localIncrementy = density[e] * (velocityVp[e] * abs( ny ) + velocityVs[e] * sqrt( pow( nx, 2 ) + pow( nz, 2 ) ) ) * aux; + real32 const localIncrementz = density[e] * (velocityVp[e] * abs( nz ) + velocityVs[e] * sqrt( pow( nx, 2 ) + pow( ny, 2 ) ) ) * aux; + + RAJA::atomicAdd< ATOMIC_POLICY >( &dampingx[facesToNodes( f, q )], localIncrementx ); + RAJA::atomicAdd< ATOMIC_POLICY >( &dampingy[facesToNodes( f, q )], localIncrementy ); + RAJA::atomicAdd< ATOMIC_POLICY >( &dampingz[facesToNodes( f, q )], localIncrementz ); + } } } - } ); // end loop over element + } ); } /// The finite element space/discretization object for the element type in the subRegion diff --git a/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEM.cpp index ca71b4c7989..01ccb470d86 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEM.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEM.cpp @@ -525,14 +525,9 @@ void ElasticWaveEquationSEM::initializePostInitialConditionsPreSubGroups() NodeManager & nodeManager = mesh.getNodeManager(); FaceManager & faceManager = mesh.getFaceManager(); + ElementRegionManager & elemManager = mesh.getElemManager(); - /// get the array of indicators: 1 if the face is on the boundary; 0 otherwise - arrayView1d< integer const > const & facesDomainBoundaryIndicator = faceManager.getDomainBoundaryIndicator(); - arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const X = nodeManager.getField< fields::referencePosition32 >().toViewConst(); - arrayView2d< real64 const > const faceNormal = faceManager.faceNormal(); - - /// get face to node map - ArrayOfArraysView< localIndex const > const facesToNodes = faceManager.nodeList().toViewConst(); + arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords = nodeManager.getField< fields::referencePosition32 >().toViewConst(); // mass matrix to be computed in this function arrayView1d< real32 > const mass = nodeManager.getField< fields::MassVector >(); @@ -546,40 +541,39 @@ void ElasticWaveEquationSEM::initializePostInitialConditionsPreSubGroups() dampingz.zero(); /// get array of indicators: 1 if face is on the free surface; 0 otherwise - arrayView1d< localIndex const > const freeSurfaceFaceIndicator = faceManager.getField< fields::FreeSurfaceFaceIndicator >(); + arrayView1d< localIndex const > const freeSurfaceFaceIndicator = faceManager.getField< fields::FreeSurfaceFaceIndicator >(); + arrayView1d< integer const > const & facesDomainBoundaryIndicator = faceManager.getDomainBoundaryIndicator(); + ArrayOfArraysView< localIndex const > const facesToNodes = faceManager.nodeList().toViewConst(); + arrayView2d< real64 const > const faceNormal = faceManager.faceNormal(); - mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, - CellElementSubRegion & elementSubRegion ) + elemManager.forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, + CellElementSubRegion & elementSubRegion ) { + finiteElement::FiniteElementBase const & + fe = elementSubRegion.getReference< finiteElement::FiniteElementBase >( getDiscretizationName() ); arrayView2d< localIndex const, cells::NODE_MAP_USD > const elemsToNodes = elementSubRegion.nodeList(); - arrayView2d< localIndex const > const facesToElements = faceManager.elementList(); - arrayView1d< real32 > const density = elementSubRegion.getField< fields::MediumDensity >(); - arrayView1d< real32 > const velocityVp = elementSubRegion.getField< fields::MediumVelocityVp >(); - arrayView1d< real32 > const velocityVs = elementSubRegion.getField< fields::MediumVelocityVs >(); + arrayView2d< localIndex const > const elemsToFaces = elementSubRegion.faceList(); - finiteElement::FiniteElementBase const & - fe = elementSubRegion.getReference< finiteElement::FiniteElementBase >( getDiscretizationName() ); + arrayView1d< real32 const > const density = elementSubRegion.getField< fields::MediumDensity >(); + arrayView1d< real32 const > const velocityVp = elementSubRegion.getField< fields::MediumVelocityVp >(); + arrayView1d< real32 const > const velocityVs = elementSubRegion.getField< fields::MediumVelocityVs >(); - finiteElement::FiniteElementDispatchHandler< SEM_FE_TYPES >::dispatch3D( fe, - [&] - ( auto const finiteElement ) + finiteElement::FiniteElementDispatchHandler< SEM_FE_TYPES >::dispatch3D( fe, [&] ( auto const finiteElement ) { using FE_TYPE = TYPEOFREF( finiteElement ); elasticWaveEquationSEMKernels::MassMatrixKernel< FE_TYPE > kernelM( finiteElement ); - kernelM.template launch< EXEC_POLICY, ATOMIC_POLICY >( elementSubRegion.size(), - X, + nodeCoords, elemsToNodes, density, mass ); elasticWaveEquationSEMKernels::DampingMatrixKernel< FE_TYPE > kernelD( finiteElement ); - - kernelD.template launch< EXEC_POLICY, ATOMIC_POLICY >( faceManager.size(), - X, - facesToElements, + kernelD.template launch< EXEC_POLICY, ATOMIC_POLICY >( elementSubRegion.size(), + nodeCoords, + elemsToFaces, facesToNodes, facesDomainBoundaryIndicator, freeSurfaceFaceIndicator, @@ -608,9 +602,9 @@ void ElasticWaveEquationSEM::applyFreeSurfaceBC( real64 const time, DomainPartit arrayView1d< real32 > const ux_np1 = nodeManager.getField< fields::Displacementx_np1 >(); arrayView1d< real32 > const uy_np1 = nodeManager.getField< fields::Displacementy_np1 >(); arrayView1d< real32 > const uz_np1 = nodeManager.getField< fields::Displacementz_np1 >(); - arrayView1d< real32 > const ux_n = nodeManager.getField< fields::Displacementx_n >(); - arrayView1d< real32 > const uy_n = nodeManager.getField< fields::Displacementy_n >(); - arrayView1d< real32 > const uz_n = nodeManager.getField< fields::Displacementz_n >(); + arrayView1d< real32 > const ux_n = nodeManager.getField< fields::Displacementx_n >(); + arrayView1d< real32 > const uy_n = nodeManager.getField< fields::Displacementy_n >(); + arrayView1d< real32 > const uz_n = nodeManager.getField< fields::Displacementz_n >(); arrayView1d< real32 > const ux_nm1 = nodeManager.getField< fields::Displacementx_nm1 >(); arrayView1d< real32 > const uy_nm1 = nodeManager.getField< fields::Displacementy_nm1 >(); arrayView1d< real32 > const uz_nm1 = nodeManager.getField< fields::Displacementz_nm1 >(); @@ -626,7 +620,7 @@ void ElasticWaveEquationSEM::applyFreeSurfaceBC( real64 const time, DomainPartit fsManager.apply( time, domain.getMeshBody( 0 ).getMeshLevel( m_discretizationName ), - string( "FreeSurface" ), + WaveSolverBase::viewKeyStruct::freeSurfaceString(), [&]( FieldSpecificationBase const & bc, string const &, SortedArrayView< localIndex const > const & targetSet, diff --git a/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEMKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEMKernel.hpp index a776dec8a6c..5f9ad660db3 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEMKernel.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEMKernel.hpp @@ -245,7 +245,7 @@ struct MassMatrixKernel arrayView1d< real32 > const mass ) { - forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k ) + forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const e ) { constexpr localIndex numNodesPerElem = FE_TYPE::numNodes; @@ -256,14 +256,14 @@ struct MassMatrixKernel { for( localIndex i = 0; i < 3; ++i ) { - xLocal[a][i] = X( elemsToNodes( k, a ), i ); + xLocal[a][i] = X( elemsToNodes( e, a ), i ); } } for( localIndex q = 0; q < numQuadraturePointsPerElem; ++q ) { - real32 const localIncrement = density[k] * m_finiteElement.computeMassTerm( q, xLocal ); - RAJA::atomicAdd< ATOMIC_POLICY >( &mass[elemsToNodes[k][q]], localIncrement ); + real32 const localIncrement = density[e] * m_finiteElement.computeMassTerm( q, xLocal ); + RAJA::atomicAdd< ATOMIC_POLICY >( &mass[elemsToNodes( e, q )], localIncrement ); } } ); // end loop over element } @@ -286,8 +286,8 @@ struct DampingMatrixKernel * @tparam EXEC_POLICY the execution policy * @tparam ATOMIC_POLICY the atomic policy * @param[in] size the number of cells in the subRegion - * @param[in] X coordinates of the nodes - * @param[in] facesToElems map from faces to elements + * @param[in] nodeCoords coordinates of the nodes + * @param[in] elemsToFaces map from elements to faces * @param[in] facesToNodes map from face to nodes * @param[in] facesDomainBoundaryIndicator flag equal to 1 if the face is on the boundary, and to 0 otherwise * @param[in] freeSurfaceFaceIndicator flag equal to 1 if the face is on the free surface, and to 0 otherwise @@ -298,8 +298,8 @@ struct DampingMatrixKernel //std::enable_if_t< geos::is_sem_formulation< std::remove_cv_t< FE_TYPE_ > >::value, void > void launch( localIndex const size, - arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const X, - arrayView2d< localIndex const > const facesToElems, + arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords, + arrayView2d< localIndex const > const elemsToFaces, ArrayOfArraysView< localIndex const > const facesToNodes, arrayView1d< integer const > const facesDomainBoundaryIndicator, arrayView1d< localIndex const > const freeSurfaceFaceIndicator, @@ -311,49 +311,39 @@ struct DampingMatrixKernel arrayView1d< real32 > const dampingy, arrayView1d< real32 > const dampingz ) { - forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const f ) + forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const e ) { - // face on the domain boundary and not on free surface - if( facesDomainBoundaryIndicator[f] == 1 && freeSurfaceFaceIndicator[f] != 1 ) + for( localIndex i = 0; i < elemsToFaces.size( 1 ); ++i ) { - localIndex k = facesToElems( f, 0 ); - if( k < 0 ) - { - k = facesToElems( f, 1 ); - } - - constexpr localIndex numNodesPerFace = FE_TYPE::numNodesPerFace; - - real64 xLocal[ numNodesPerFace ][ 3 ]; - for( localIndex a = 0; a < numNodesPerFace; ++a ) + localIndex const f = elemsToFaces( e, i ); + // face on the domain boundary and not on free surface + if( facesDomainBoundaryIndicator[f] == 1 && freeSurfaceFaceIndicator[f] != 1 ) { - for( localIndex i = 0; i < 3; ++i ) + constexpr localIndex numNodesPerFace = FE_TYPE::numNodesPerFace; + real64 xLocal[ numNodesPerFace ][ 3 ]; + for( localIndex a = 0; a < numNodesPerFace; ++a ) { - xLocal[a][i] = X( facesToNodes( f, a ), i ); + for( localIndex d = 0; d < 3; ++d ) + { + xLocal[a][d] = nodeCoords( facesToNodes( f, a ), d ); + } } - } - for( localIndex q = 0; q < numNodesPerFace; ++q ) - { - real32 const localIncrementx = density[k] * (velocityVp[k]*LvArray::math::abs( faceNormal[f][0] ) + velocityVs[k]*sqrt( faceNormal[f][1]*faceNormal[f][1] + - faceNormal[f][2]*faceNormal[f][2] ) )* m_finiteElement.computeDampingTerm( - q, - xLocal ); - real32 const localIncrementy = density[k] * (velocityVp[k]*LvArray::math::abs( faceNormal[f][1] ) + velocityVs[k]*sqrt( faceNormal[f][0]*faceNormal[f][0] + - faceNormal[f][2]*faceNormal[f][2] ) )* m_finiteElement.computeDampingTerm( - q, - xLocal ); - real32 const localIncrementz = density[k] * (velocityVp[k]*LvArray::math::abs( faceNormal[f][2] ) + velocityVs[k]*sqrt( faceNormal[f][0]*faceNormal[f][0] + - faceNormal[f][1]*faceNormal[f][1] ) )* m_finiteElement.computeDampingTerm( - q, - xLocal ); - - RAJA::atomicAdd< ATOMIC_POLICY >( &dampingx[facesToNodes[f][q]], localIncrementx ); - RAJA::atomicAdd< ATOMIC_POLICY >( &dampingy[facesToNodes[f][q]], localIncrementy ); - RAJA::atomicAdd< ATOMIC_POLICY >( &dampingz[facesToNodes[f][q]], localIncrementz ); + real32 const nx = faceNormal( f, 0 ), ny = faceNormal( f, 1 ), nz = faceNormal( f, 2 ); + for( localIndex q = 0; q < numNodesPerFace; ++q ) + { + real32 const aux = density[e] * m_finiteElement.computeDampingTerm( q, xLocal ); + real32 const localIncrementx = aux * ( velocityVp[e] * abs( nx ) + velocityVs[e] * sqrt( pow( ny, 2 ) + pow( nz, 2 ) ) ); + real32 const localIncrementy = aux * ( velocityVp[e] * abs( ny ) + velocityVs[e] * sqrt( pow( nx, 2 ) + pow( nz, 2 ) ) ); + real32 const localIncrementz = aux * ( velocityVp[e] * abs( nz ) + velocityVs[e] * sqrt( pow( nx, 2 ) + pow( ny, 2 ) ) ); + + RAJA::atomicAdd< ATOMIC_POLICY >( &dampingx[facesToNodes( f, q )], localIncrementx ); + RAJA::atomicAdd< ATOMIC_POLICY >( &dampingy[facesToNodes( f, q )], localIncrementy ); + RAJA::atomicAdd< ATOMIC_POLICY >( &dampingz[facesToNodes( f, q )], localIncrementz ); + } } } - } ); // end loop over element + } ); } /// The finite element space/discretization object for the element type in the subRegion diff --git a/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverBase.cpp b/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverBase.cpp index 793d3f27bb4..9f3840c27c5 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverBase.cpp @@ -394,4 +394,11 @@ localIndex WaveSolverBase::getNumNodesPerElem() } +bool WaveSolverBase::directoryExists( std::string const & directoryName ) +{ + struct stat buffer; + return stat( directoryName.c_str(), &buffer ) == 0; +} + + } /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverBase.hpp b/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverBase.hpp index 14beb4b005f..f1ea47acb54 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverBase.hpp @@ -112,6 +112,7 @@ class WaveSolverBase : public SolverBase static constexpr char const * usePMLString() { return "usePML"; } static constexpr char const * parametersPMLString() { return "parametersPML"; } + static constexpr char const * freeSurfaceString() { return "FreeSurface"; } }; /** @@ -128,6 +129,13 @@ class WaveSolverBase : public SolverBase virtual void postProcessInput() override; + /** + * @brief Utility function to check if a directory exists + * @param directoryName the name of the directory + * @return true if the directory exists, false otherwise + */ + bool directoryExists( std::string const & directoryName ); + /** * @brief Apply free surface condition to the face defined in the geometry box of the xml * @param time the time to apply the BC diff --git a/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverBaseFields.hpp b/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverBaseFields.hpp index 0b04145b044..c55fdb8cc2d 100644 --- a/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverBaseFields.hpp +++ b/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverBaseFields.hpp @@ -14,11 +14,11 @@ /** - * @file AcousticFirstOrderWaveEquationSEM.hpp + * @file WaveSolverBaseFields.hpp */ -#ifndef GEOS_PHYSICSSOLVERS_WAVEPROPAGATION__HPP_WAVESOLVERBASEFIELDS -#define GEOS_PHYSICSSOLVERS_WAVEPROPAGATION__HPP_WAVESOLVERBASEFIELDS +#ifndef GEOS_PHYSICSSOLVERS_WAVEPROPAGATION_WAVESOLVERBASEFIELDS_HPP_ +#define GEOS_PHYSICSSOLVERS_WAVEPROPAGATION_WAVESOLVERBASEFIELDS_HPP_ #include "common/DataLayouts.hpp" #include "mesh/MeshFields.hpp" @@ -32,6 +32,46 @@ namespace fields namespace wavesolverfields { +DECLARE_FIELD( Delta, + "delta", + array1d< real32 >, + 0, + NOPLOT, + WRITE_AND_READ, + "Delta thomsen anisotropy parameter" ); + +DECLARE_FIELD( Epsilon, + "epsilon", + array1d< real32 >, + 0, + NOPLOT, + WRITE_AND_READ, + "Epsilon thomsen anisotropy parameter" ); + +DECLARE_FIELD( F, + "f", + array1d< real32 >, + 0, + NOPLOT, + WRITE_AND_READ, + "f quantity in VTI/TTI Fletcher's equations" ); + +DECLARE_FIELD( StiffnessVector_p, + "stiffnessVector_p", + array1d< real32 >, + 0, + NOPLOT, + WRITE_AND_READ, + "Stiffness vector contains R_h*Pressure_n." ); + +DECLARE_FIELD( StiffnessVector_q, + "stiffnessVector_q", + array1d< real32 >, + 0, + NOPLOT, + WRITE_AND_READ, + "Stiffness vector contains R_h*Pressure_n." ); + DECLARE_FIELD( Pressure_np1, "pressure_np1", array1d< real32 >, @@ -40,6 +80,54 @@ DECLARE_FIELD( Pressure_np1, WRITE_AND_READ, "Scalar pressure at time n+1." ); +DECLARE_FIELD( Pressure_p_nm1, + "pressure_p_nm1", + array1d< real32 >, + 0, + NOPLOT, + WRITE_AND_READ, + "Scalar pressure at time n-1." ); + +DECLARE_FIELD( Pressure_p_n, + "pressure_p_n", + array1d< real32 >, + 0, + NOPLOT, + WRITE_AND_READ, + "Scalar pressure at time n." ); + +DECLARE_FIELD( Pressure_p_np1, + "pressure_p_np1", + array1d< real32 >, + 0, + LEVEL_0, + WRITE_AND_READ, + "Scalar pressure at time n+1." ); + +DECLARE_FIELD( Pressure_q_nm1, + "pressure_q_nm1", + array1d< real32 >, + 0, + NOPLOT, + WRITE_AND_READ, + "Scalar auxiliary pressure q at time n-1." ); + +DECLARE_FIELD( Pressure_q_n, + "pressure_q_n", + array1d< real32 >, + 0, + NOPLOT, + WRITE_AND_READ, + "Scalar auxiliary pressure q at time n." ); + +DECLARE_FIELD( Pressure_q_np1, + "pressure_q_np1", + array1d< real32 >, + 0, + LEVEL_0, + WRITE_AND_READ, + "Scalar auxiliary pressure q at time n+1." ); + DECLARE_FIELD( Velocity_x, "velocity_x", array2d< real32 >, @@ -88,6 +176,38 @@ DECLARE_FIELD( DampingVector, WRITE_AND_READ, "Diagonal of the Damping Matrix." ); +DECLARE_FIELD( DampingVector_p, + "dampingVector_p", + array1d< real32 >, + 0, + NOPLOT, + WRITE_AND_READ, + "Diagonal of the Damping Matrix for p terms in p equation." ); + +DECLARE_FIELD( DampingVector_pq, + "dampingVector_pq", + array1d< real32 >, + 0, + NOPLOT, + WRITE_AND_READ, + "Diagonal of the Damping Matrix for q terms in p equation." ); + +DECLARE_FIELD( DampingVector_q, + "dampingVector_q", + array1d< real32 >, + 0, + NOPLOT, + WRITE_AND_READ, + "Diagonal of the Damping Matrix for q terms in q equation." ); + +DECLARE_FIELD( DampingVector_qp, + "dampingVector_qp", + array1d< real32 >, + 0, + NOPLOT, + WRITE_AND_READ, + "Diagonal of the Damping Matrix for p terms in q equation." ); + DECLARE_FIELD( MediumVelocity, "mediumVelocity", array1d< real32 >, @@ -120,6 +240,38 @@ DECLARE_FIELD( FreeSurfaceNodeIndicator, WRITE_AND_READ, "Free surface indicator, 1 if a node is on free surface 0 otherwise." ); +DECLARE_FIELD( LateralSurfaceFaceIndicator, + "lateralSurfaceFaceIndicator", + array1d< localIndex >, + 0, + NOPLOT, + WRITE_AND_READ, + "Free surface indicator, 1 if a face is on a lateral surface 0 otherwise." ); + +DECLARE_FIELD( LateralSurfaceNodeIndicator, + "lateralSurfaceNodeIndicator", + array1d< localIndex >, + 0, + NOPLOT, + WRITE_AND_READ, + "Lateral surface indicator, 1 if a face is on a lateral surface 0 otherwise." ); + +DECLARE_FIELD( BottomSurfaceFaceIndicator, + "bottomSurfaceFaceIndicator", + array1d< localIndex >, + 0, + NOPLOT, + WRITE_AND_READ, + "Bottom surface indicator, 1 if a face is on the bottom surface 0 otherwise." ); + +DECLARE_FIELD( BottomSurfaceNodeIndicator, + "bottomSurfaceNodeIndicator", + array1d< localIndex >, + 0, + NOPLOT, + WRITE_AND_READ, + "Bottom surface indicator, 1 if a face is on the bottom surface 0 otherwise." ); + DECLARE_FIELD( Displacementx_np1, "displacementx_np1", array1d< real32 >, diff --git a/src/coreComponents/schema/docs/AcousticVTISEM.rst b/src/coreComponents/schema/docs/AcousticVTISEM.rst new file mode 100644 index 00000000000..000b0ace852 --- /dev/null +++ b/src/coreComponents/schema/docs/AcousticVTISEM.rst @@ -0,0 +1,30 @@ + + +========================= ============== ========== ======================================================================================================================================================================================================================================================================================================================== +Name Type Default Description +========================= ============== ========== ======================================================================================================================================================================================================================================================================================================================== +cflFactor real64 0.5 Factor to apply to the `CFL condition `_ when calculating the maximum allowable time step. Values should be in the interval (0,1] +discretization string required Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified. +dtSeismoTrace real64 0 Time step for output pressure at receivers +enableLifo integer 0 Set to 1 to enable LIFO storage feature +forward integer 1 Set to 1 to compute forward propagation +initialDt real64 1e+99 Initial time-step value required by the solver to the event manager. +lifoOnDevice integer -80 Set the capacity of the lifo device storage (if negative, opposite of percentage of remaining memory) +lifoOnHost integer -80 Set the capacity of the lifo host storage (if negative, opposite of percentage of remaining memory) +lifoSize integer 2147483647 Set the capacity of the lifo storage (should be the total number of buffers to store in the LIFO) +linearDASGeometry real64_array2d {{0}} Geometry parameters for a linear DAS fiber (dip, azimuth, gauge length) +logLevel integer 0 Log level +name string required A name is required for any non-unique nodes +outputSeismoTrace integer 0 Flag that indicates if we write the seismo trace in a file .txt, 0 no output, 1 otherwise +receiverCoordinates real64_array2d required Coordinates (x,y,z) of the receivers +rickerOrder integer 2 Flag that indicates the order of the Ricker to be used o, 1 or 2. Order 2 by default +saveFields integer 0 Set to 1 to save fields during forward and restore them during backward +shotIndex integer 0 Set the current shot for temporary files +sourceCoordinates real64_array2d required Coordinates (x,y,z) of the sources +targetRegions string_array required Allowable regions that the solver may be applied to. Note that this does not indicate that the solver will be applied to these regions, only that allocation will occur such that the solver may be applied to these regions. The decision about what regions this solver will beapplied to rests in the EventManager. +timeSourceFrequency real32 required Central frequency for the time source +LinearSolverParameters node unique :ref:`XML_LinearSolverParameters` +NonlinearSolverParameters node unique :ref:`XML_NonlinearSolverParameters` +========================= ============== ========== ======================================================================================================================================================================================================================================================================================================================== + + diff --git a/src/coreComponents/schema/docs/AcousticVTISEM_other.rst b/src/coreComponents/schema/docs/AcousticVTISEM_other.rst new file mode 100644 index 00000000000..a5f09a365c9 --- /dev/null +++ b/src/coreComponents/schema/docs/AcousticVTISEM_other.rst @@ -0,0 +1,23 @@ + + +========================= ============================================================================================================================================================== ======================================================================= +Name Type Description +========================= ============================================================================================================================================================== ======================================================================= +indexSeismoTrace integer Count for output pressure at receivers +maxStableDt real64 Value of the Maximum Stable Timestep for this solver. +meshTargets geos_mapBase< std_pair< string, string >, LvArray_Array< string, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. +pressureNp1AtReceivers real32_array2d Pressure value at each receiver for each timestep +receiverIsLocal integer_array Flag that indicates whether the receiver is local to this MPI rank +receiverNodeIds integer_array2d Indices of the nodes (in the right order) for each receiver point +sourceConstants real64_array2d Constant part of the receiver for the nodes listed in m_receiverNodeIds +sourceIsAccessible integer_array Flag that indicates whether the source is local to this MPI rank +sourceNodeIds integer_array2d Indices of the nodes (in the right order) for each source point +sourceValue real32_array2d Source Value of the sources +useDAS integer Flag to indicate if DAS type of data will be modeled +usePML integer Flag to apply PML +LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` +NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` +SolverStatistics node :ref:`DATASTRUCTURE_SolverStatistics` +========================= ============================================================================================================================================================== ======================================================================= + + diff --git a/src/coreComponents/schema/docs/Solvers.rst b/src/coreComponents/schema/docs/Solvers.rst index 0983f14a28f..d153cb46725 100644 --- a/src/coreComponents/schema/docs/Solvers.rst +++ b/src/coreComponents/schema/docs/Solvers.rst @@ -6,6 +6,7 @@ Name Type Default Description gravityVector R1Tensor {0,0,-9.81} Gravity vector used in the physics solvers AcousticFirstOrderSEM node :ref:`XML_AcousticFirstOrderSEM` AcousticSEM node :ref:`XML_AcousticSEM` +AcousticVTISEM node :ref:`XML_AcousticVTISEM` CompositionalMultiphaseFVM node :ref:`XML_CompositionalMultiphaseFVM` CompositionalMultiphaseHybridFVM node :ref:`XML_CompositionalMultiphaseHybridFVM` CompositionalMultiphaseReservoir node :ref:`XML_CompositionalMultiphaseReservoir` diff --git a/src/coreComponents/schema/docs/Solvers_other.rst b/src/coreComponents/schema/docs/Solvers_other.rst index c30e95606f1..aff8135dedc 100644 --- a/src/coreComponents/schema/docs/Solvers_other.rst +++ b/src/coreComponents/schema/docs/Solvers_other.rst @@ -5,6 +5,7 @@ Name Type Description =========================================== ==== ================================================================ AcousticFirstOrderSEM node :ref:`DATASTRUCTURE_AcousticFirstOrderSEM` AcousticSEM node :ref:`DATASTRUCTURE_AcousticSEM` +AcousticVTISEM node :ref:`DATASTRUCTURE_AcousticVTISEM` CompositionalMultiphaseFVM node :ref:`DATASTRUCTURE_CompositionalMultiphaseFVM` CompositionalMultiphaseHybridFVM node :ref:`DATASTRUCTURE_CompositionalMultiphaseHybridFVM` CompositionalMultiphaseReservoir node :ref:`DATASTRUCTURE_CompositionalMultiphaseReservoir` diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd index 6e073e4a00d..e823dde2bc2 100644 --- a/src/coreComponents/schema/schema.xsd +++ b/src/coreComponents/schema/schema.xsd @@ -298,6 +298,10 @@ + + + + @@ -1914,6 +1918,7 @@ the relative residual norm satisfies: + @@ -2050,6 +2055,52 @@ the relative residual norm satisfies: + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/coreComponents/schema/schema.xsd.other b/src/coreComponents/schema/schema.xsd.other index d6e5a858135..7f1d776a570 100644 --- a/src/coreComponents/schema/schema.xsd.other +++ b/src/coreComponents/schema/schema.xsd.other @@ -474,6 +474,7 @@ + @@ -600,6 +601,37 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/coreComponents/unitTests/constitutiveTests/testCO2BrinePVTModels.cpp b/src/coreComponents/unitTests/constitutiveTests/testCO2BrinePVTModels.cpp index 34c57b1d4f1..24cee2654a9 100644 --- a/src/coreComponents/unitTests/constitutiveTests/testCO2BrinePVTModels.cpp +++ b/src/coreComponents/unitTests/constitutiveTests/testCO2BrinePVTModels.cpp @@ -744,7 +744,6 @@ TEST_F( EzrokhiBrineDensityTest, brineCO2DensityMolarValuesAndDeriv ) } } - class SpanWagnerCO2DensityTest : public ::testing::Test { public: @@ -844,7 +843,6 @@ TEST_F( SpanWagnerCO2DensityTest, spanWagnerCO2DensityMolarValuesAndDeriv ) } } - class CO2SolubilityTest : public ::testing::Test { public: @@ -945,12 +943,10 @@ TEST_F( BrineEnthalpyTest, BrineEnthalpyMassValuesAndDeriv ) BrineEnthalpy::KernelWrapper pvtFunctionWrapper = pvtFunction->createKernelWrapper(); - real64 const savedValues[] = { 279338.177769, 281327.917595, 282984.982790, 273523.601187, 275547.730627, - 277232.967486, 259305.594861, 261448.359022, 263229.980673, 207243.035969, - 208857.558651, 210202.000531, 197603.080058, 199274.617099, 200665.764632, - 174031.122202, 175899.343122, 177450.286494, 135147.894170, 136387.199707, - 137419.018273, 121682.558928, 123001.503570, 124098.561778, 88756.649542, - 90350.327221, 91670.592316 }; + real64 const savedValues[] = + { 433114, 435103, 436760, 427299, 429323, 431008, 413081, 415224, 417005, 462186, 463801, 465145, 452546, + 454218, 455609, 428974, 430843, 432394, 491259, 492499, 493530, 477794, 479113, 480210, 444868, 446462, 447782 }; + integer counter = 0; for( integer iComp = 0; iComp < 3; ++iComp ) { @@ -982,14 +978,12 @@ TEST_F( BrineEnthalpyTest, BrineEnthalpyMolarValuesAndDeriv ) real64 const eps = sqrt( std::numeric_limits< real64 >::epsilon()); real64 const relTol = 5e-5; - - BrineEnthalpy::KernelWrapper pvtFunctionWrapper = pvtFunction->createKernelWrapper(); - real64 const savedValues[] = { 15234891.499346, 15338606.577025, 15424985.891636, 15102742.031582, 15207238.691387, 15294258.271084, - 14779605.524173, 14886798.427634, 14976008.570783, 11042832.057960, 11121210.933610, 11186485.534383, - 10823742.150877, 10903416.807419, 10969752.900309, 10288015.835964, 10372160.580672, 10442128.397178, - 6850772.616574, 6903815.290194, 6947985.177129, 6544742.270173, 6599594.923452, 6645247.529535, - 5796426.147754, 5857522.733710, 5908248.223573}; + real64 const savedValues[] = + { 1.87298e+07, 1.88335e+07, 1.89199e+07, 1.85976e+07, 1.87021e+07, 1.87892e+07, 1.82745e+07, 1.83817e+07, 1.84709e+07, + 1.6837e+07, 1.69154e+07, 1.69807e+07, 1.66179e+07, 1.66976e+07, 1.67639e+07, 1.60822e+07, 1.61663e+07, 1.62363e+07, + 1.49442e+07, 1.49973e+07, 1.50414e+07, 1.46382e+07, 1.4693e+07, 1.47387e+07, 1.38899e+07, 1.3951e+07, 1.40017e+07 }; + integer counter = 0; for( integer iComp = 0; iComp < 3; ++iComp ) { @@ -1008,7 +1002,6 @@ TEST_F( BrineEnthalpyTest, BrineEnthalpyMolarValuesAndDeriv ) } } - class CO2EnthalpyTest : public ::testing::Test { public: @@ -1044,13 +1037,10 @@ TEST_F( CO2EnthalpyTest, CO2EnthalpyMassValuesAndDeriv ) real64 const eps = sqrt( std::numeric_limits< real64 >::epsilon()); real64 const relTol = 5e-5; - - CO2Enthalpy::KernelWrapper pvtFunctionWrapper = pvtFunction->createKernelWrapper(); - real64 const savedValues[] = { 28447.084306, 29131.068469, 29700.204529, 9320.187656, 10117.295548, 10779.101555, -37449.569995, - -36262.216311, -35283.355068, 28447.084306, 29131.068469, 29700.204529, 9320.187656, 10117.295548, - 10779.101555, -37449.569995, -36262.216311, -35283.355068, 28447.084306, 29131.068469, 29700.204529, - 9320.187656, 10117.295548, 10779.101555, -37449.569995, -36262.216311, -35283.355068}; + real64 const savedValues[] = + { 534287, 534971, 535540, 515160, 515957, 516619, 468390, 469578, 470557, 534287, 534971, 535540, 515160, 515957, + 516619, 468390, 469578, 470557, 534287, 534971, 535540, 515160, 515957, 516619, 468390, 469578, 470557 }; integer counter = 0; for( integer iComp = 0; iComp < 3; ++iComp ) @@ -1083,14 +1073,14 @@ TEST_F( CO2EnthalpyTest, CO2EnthalpyMolarValuesAndDeriv ) real64 const eps = sqrt( std::numeric_limits< real64 >::epsilon()); real64 const relTol = 5e-5; - - CO2Enthalpy::KernelWrapper pvtFunctionWrapper = pvtFunction->createKernelWrapper(); - real64 const savedValues[] = { 646524.643323, 662069.737939, 675004.648394, 211822.446731, 229938.535180, 244979.580788, - -851126.590796, -824141.279795, -801894.433361, 646524.643323, 662069.737939, 675004.648394, - 211822.446731, 229938.535180, 244979.580788, -851126.590796, -824141.279795, -801894.433361, - 646524.643323, 662069.737939, 675004.648394, 211822.446731, 229938.535180, 244979.580788, - -851126.590796, -824141.279795, -801894.433361 }; + real64 const savedValues[] = + { 12142888.279686311, 12158433.374301625, 12171368.284757353, 11708186.083185773, 11726302.171616826, 11741343.217216015, + 10645237.045567034, 10672222.35656886, 10694469.203002082, 12142888.279686311, 12158433.374301625, 12171368.284757353, + 11708186.083185773, 11726302.171616826, 11741343.217216015, 10645237.045567034, 10672222.35656886, 10694469.203002082, + 12142888.279686311, 12158433.374301625, 12171368.284757353, 11708186.083185773, 11726302.171616826, 11741343.217216015, + 10645237.045567034, 10672222.35656886, 10694469.203002082 }; + integer counter = 0; for( integer iComp = 0; iComp < 3; ++iComp ) { @@ -1109,7 +1099,6 @@ TEST_F( CO2EnthalpyTest, CO2EnthalpyMolarValuesAndDeriv ) } } - int main( int argc, char * * argv ) { ::testing::InitGoogleTest( &argc, argv ); diff --git a/src/docs/sphinx/CompleteXMLSchema.rst b/src/docs/sphinx/CompleteXMLSchema.rst index eb52d4dbff6..828584302fe 100644 --- a/src/docs/sphinx/CompleteXMLSchema.rst +++ b/src/docs/sphinx/CompleteXMLSchema.rst @@ -1,9 +1,11 @@ -====================== +###################### Datastructure Index -====================== +###################### + +************************** Input Schema Definitions -======================== +************************** :download:`XML Schema <../../coreComponents/schema/docs/../schema.xsd>` @@ -22,6 +24,13 @@ Element: AcousticSEM .. include:: ../../coreComponents/schema/docs/AcousticSEM.rst +.. _XML_AcousticVTISEM: + +Element: AcousticVTISEM +======================= +.. include:: ../../coreComponents/schema/docs/AcousticVTISEM.rst + + .. _XML_Aquifer: Element: Aquifer @@ -1251,6 +1260,13 @@ Datastructure: AcousticSEM .. include:: ../../coreComponents/schema/docs/AcousticSEM_other.rst +.. _DATASTRUCTURE_AcousticVTISEM: + +Datastructure: AcousticVTISEM +============================= +.. include:: ../../coreComponents/schema/docs/AcousticVTISEM_other.rst + + .. _DATASTRUCTURE_Aquifer: Datastructure: Aquifer diff --git a/src/docs/sphinx/advancedExamples/validationStudies/thermoPoromechanics/thermalConsolidation/Example.rst b/src/docs/sphinx/advancedExamples/validationStudies/thermoPoromechanics/thermalConsolidation/Example.rst index 7ef20823419..74e9a10750d 100644 --- a/src/docs/sphinx/advancedExamples/validationStudies/thermoPoromechanics/thermalConsolidation/Example.rst +++ b/src/docs/sphinx/advancedExamples/validationStudies/thermoPoromechanics/thermalConsolidation/Example.rst @@ -20,7 +20,7 @@ This example uses no external input files and everything required is contained w .. code-block:: console - inputFiles/thermoPoromechanics/ThermoPoroElastic_consolidation_benchmark.xml + inputFiles/thermoPoromechanics/ThermoPoroElastic_consolidation_benchmark_fim.xml --------------------------------------------------- Description of the case diff --git a/src/docs/sphinx/basicExamples/co2Injection/Example.rst b/src/docs/sphinx/basicExamples/co2Injection/Example.rst index a0f5a4c81a4..cec5f2b464c 100644 --- a/src/docs/sphinx/basicExamples/co2Injection/Example.rst +++ b/src/docs/sphinx/basicExamples/co2Injection/Example.rst @@ -260,24 +260,24 @@ The simulation can be launched with 4 cores using MPI-parallelism: .. code-block:: console - mpirun -np 4 geosx -i SimpleCo2InjTutorial.xml -x 1 -y 1 -z 4 + mpirun -np 4 geosx -i simpleCo2InjTutorial.xml -x 1 -y 1 -z 4 -A restart from a checkpoint file `SimpleCo2InjTutorial_restart_000000024.root` is always available thanks to the following command line : +A restart from a checkpoint file `simpleCo2InjTutorial_restart_000000024.root` is always available thanks to the following command line : .. code-block:: console - mpirun -np 4 geosx -i SimpleCo2InjTutorial.xml -r SimpleCo2InjTutorial_restart_000000024 -x 1 -y 1 -z 4 + mpirun -np 4 geosx -i simpleCo2InjTutorial.xml -r simpleCo2InjTutorial_restart_000000024 -x 1 -y 1 -z 4 The output then shows the loading of HDF5 restart files by each core. .. code-block:: console - Loading restart file SimpleCo2InjTutorial_restart_000000024 - Rank 0: rankFilePattern = SimpleCo2InjTutorial_restart_000000024/rank_%07d.hdf5 - Rank 0: Reading in restart file at SimpleCo2InjTutorial_restart_000000024/rank_0000000.hdf5 - Rank 1: Reading in restart file at SimpleCo2InjTutorial_restart_000000024/rank_0000001.hdf5 - Rank 3: Reading in restart file at SimpleCo2InjTutorial_restart_000000024/rank_0000003.hdf5 - Rank 2: Reading in restart file at SimpleCo2InjTutorial_restart_000000024/rank_0000002.hdf5 + Loading restart file simpleCo2InjTutorial_restart_000000024 + Rank 0: rankFilePattern = simpleCo2InjTutorial_restart_000000024/rank_%07d.hdf5 + Rank 0: Reading in restart file at simpleCo2InjTutorial_restart_000000024/rank_0000000.hdf5 + Rank 1: Reading in restart file at simpleCo2InjTutorial_restart_000000024/rank_0000001.hdf5 + Rank 3: Reading in restart file at simpleCo2InjTutorial_restart_000000024/rank_0000003.hdf5 + Rank 2: Reading in restart file at simpleCo2InjTutorial_restart_000000024/rank_0000002.hdf5 and the simulation restarts from this point in time. diff --git a/src/docs/sphinx/basicExamples/multiphaseFlowWithWells/Example.rst b/src/docs/sphinx/basicExamples/multiphaseFlowWithWells/Example.rst index f49b2cdabb3..f0fdd28eece 100644 --- a/src/docs/sphinx/basicExamples/multiphaseFlowWithWells/Example.rst +++ b/src/docs/sphinx/basicExamples/multiphaseFlowWithWells/Example.rst @@ -27,14 +27,14 @@ This example is based on the XML file located at .. code-block:: console - ../../../../../inputFiles/compositionalMultiphaseWell/benchmarks/Egg/deadOilEgg_base_iterative.xml + ../../../../../inputFiles/compositionalMultiphaseWell/benchmarks/Egg/deadOilEgg_benchmark.xml The mesh file corresponding to the Egg model is stored in the GEOSXDATA repository. Therefore, you must first download the GEOSXDATA repository in the same folder as the GEOS repository to run this test case. .. note:: - GEOSXDATA is a separate repository in which we store large mesh files in order to keep the main GEOS repository lightweight. + `GEOSXDATA `_ is a separate repository in which we store large mesh files in order to keep the main GEOS repository lightweight. The XML file considered here follows the typical structure of the GEOS input files: diff --git a/src/docs/sphinx/requirements.txt b/src/docs/sphinx/requirements.txt index 97f1025dbf7..0fb4f2f3c11 100644 --- a/src/docs/sphinx/requirements.txt +++ b/src/docs/sphinx/requirements.txt @@ -3,12 +3,16 @@ matplotlib numpy h5py mpmath -docutils<0.18 +docutils>=0.18 pandas # using plantuml for diagrams +Sphinx>=7.0.0 +# pydata-sphinx-theme +sphinx_rtd_theme sphinxcontrib-plantuml sphinx-argparse +sphinx-design # Running CLI programs and capture outputs -sphinxcontrib-programoutput==0.17 +sphinxcontrib-programoutput>=0.17 # Installing the mesh_doctor requirements to be able to load all the modules and run the help. -r ../../coreComponents/python/modules/geosx_mesh_doctor/requirements.txt \ No newline at end of file diff --git a/src/index.rst b/src/index.rst index d3093503a63..6fe0dd5862e 100644 --- a/src/index.rst +++ b/src/index.rst @@ -1,7 +1,6 @@ +######################## GEOS Documentation -================================= - -Welcome to our documentation pages! +######################## GEOS is a code framework focused on enabling streamlined development of physics simulations on high performance computing platforms. Our documentation @@ -26,6 +25,80 @@ High quality documentation is a critical component of a successful code. If you have suggestions for improving the guides below, please post an issue on our `issue tracker `_. + + +.. grid:: 2 + :gutter: 4 + + .. grid-item-card:: + + Quick Start Guide + ^^^^^^^^^^^^^^^^^^^ + + New to GEOS? We will walk you through downloading the source, compiling the code, and testing the installation. + + +++ + + .. button-ref:: QuickStart + :expand: + :color: info + :click-parent: + + To the Quick Start + + .. grid-item-card:: + + Tutorials + ^^^^^^^^^^ + + Working tutorials that show how to run some common problems. After going through these examples, you should have a good understanding of how to set up and solve your own models. + + +++ + + .. button-ref:: Tutorials + :expand: + :color: info + :click-parent: + + To the Tutorials + + .. grid-item-card:: + + Basic Examples + ^^^^^^^^^^^^^^^ + + Example problems that are organized around physical processes (fluid flow, mechanics, etc.). + + +++ + + .. button-ref:: BasicExamples + :expand: + :color: info + :click-parent: + + To the Basic Examples + + .. grid-item-card:: + + Advanced Examples + ^^^^^^^^^^^^^^^^^^^ + + Example problems that demonstrate additional physical models, constitutive models, advanced features, etc. + + +++ + + .. button-ref:: AdvancedExamples + :expand: + :color: info + :click-parent: + + To the Advanced Examples + + +******************** +Table of Contents +******************** + .. toctree:: :maxdepth: 2 @@ -55,8 +128,9 @@ you have suggestions for improving the guides below, please post an issue on our docs/sphinx/Acknowledgements +********************* Indices and tables -================== +********************* * :ref:`genindex` * :ref:`modindex`