From c347dc0bb0cd658c057986f62750b7b2a2d79121 Mon Sep 17 00:00:00 2001 From: MelReyCG Date: Wed, 26 Nov 2025 09:32:49 +0100 Subject: [PATCH 01/13] =?UTF-8?q?=F0=9F=93=9D=20Code=20Rules=20page=20firs?= =?UTF-8?q?t=20draft=20version?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../developerGuide/Contributing/CodeRules.rst | 1244 +++++++++++++++++ 1 file changed, 1244 insertions(+) create mode 100644 src/docs/sphinx/developerGuide/Contributing/CodeRules.rst diff --git a/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst b/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst new file mode 100644 index 00000000000..84aed982b7e --- /dev/null +++ b/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst @@ -0,0 +1,1244 @@ +############################################################################### +Code Rules +############################################################################### + +.. note:: + This document uses collapsible code blocks. Click on "Show/Hide Code" to expand examples. + +Introduction +============ +This document provides comprehensive coding rules and best practices for developing with GEOS. These rules complement the :doc:`CodeStyle` guidelines and focus on typing, error handling, parallelism, performance, and general architecture principles. + +All contributors should familiarize themselves with these rules to maintain code quality, consistency, and reliability across the GEOS codebase. + +Typing +====== + +Type Correspondence Table +------------------------- + +GEOS defines its own type system for **portability**, **clarity**, and **performance**. + +**Always use GEOS types** instead of standard C++ types. + +Basic Types +^^^^^^^^^^^ + +.. list-table:: Standard C++ vs GEOS Types + :header-rows: 1 + :widths: 30 30 40 + + * - Standard C++ Type + - GEOS Type + - Description + * - ``int`` + - ``integer`` + - Signed integer + * - ``std::size_t`` + - ``localIndex`` + - Local array indexing (MPI partition) + * - N/A + - ``globalIndex`` + - Global indexing (across MPI) + * - ``float`` + - ``real32`` + - 32-bit floating point + * - ``double`` + - ``real64`` + - 64-bit floating point + * - ``std::string`` + - ``string`` + - String type + * - ``std::string_view`` + - ``string_view`` + - Non-owning string view + +.. dropdown:: Example: Using GEOS basic types + :icon: code + + .. code-block:: c++ + + // BAD - Do not use standard types + int count = 0; + double value = 1.5; + + // GOOD - Use GEOS types + integer count = 0; + real64 value = 1.5; + localIndex elementIndex = 42; + +Array Types +^^^^^^^^^^^ + +GEOS provides multi-dimensional CHAI array types for managed host-device data. + +Refer to the rule of :ref:`CHAI Memory Management`. + +.. list-table:: GEOS Array Types + :header-rows: 1 + :widths: 35 65 + + * - GEOS Type + - Description + * - ``array1d`` + - 1D owning array + * - ``arrayView1d`` + - 1D non-owning view (modifiable) + * - ``arraySlice1d`` + - 1D slice (from multi-dim array) + * - ``array2d`` + - 2D owning array + * - ``arrayView2d`` + - 2D non-owning view + * - ``arraySlice2d`` + - 2D slice + * - ``array3d``, ``array4d``, ``array5d`` + - 3D/4D/5D owning arrays + * - ``arrayView3d``, etc. + - 3D/4D/5D views + * - ``SortedArray`` + - Sorted array with search capabilities + * - ``ArrayOfArrays`` + - Jagged 2D array (variable row sizes) + * - ``ArrayOfSets`` + - Array of sets (unique elements per row) + * - ``stackArray1d`` + - Stack-allocated array (max size N) + +.. dropdown:: Example: Using GEOS arrays + :icon: code + + .. code-block:: c++ + + // BAD - Do not use std::vector + std::vector values; + + // GOOD - Use GEOS arrays + array1d values; + arrayView1d valuesView = values.toView(); + arrayView1d constView = values.toViewConst(); + +Tensor Types +^^^^^^^^^^^^ + +Use GEOS tensor types for geometric and mechanical calculations: + +.. list-table:: GEOS Tensor Types + :header-rows: 1 + :widths: 30 70 + + * - Type + - Description + * - ``R1Tensor`` + - 3D vector (real64) + * - ``R1Tensor32`` + - 3D vector (real32) + * - ``R2SymTensor`` + - Symmetric 6-component tensor (Voigt notation) + +``std`` Standard Library +------------------------ + +It is possible to use ``std`` library components for data on host memory, for doing so, the rule is +to **only use GEOS container wrappers** instead of direct standard library containers. + + * - Standard C++ Type + - GEOS Type + - Description + * - ``std::vector`` + - ``stdVector`` + - ``std`` dynamically allocated array + * - ``std::array`` + - ``stdArray`` + - ``std`` fixed constexpr sized array + * - ``std::map`` + - ``map`` + - ``std`` sorted dictionary + * - ``std::unordered_map`` + - ``unordered_map`` + - ``std`` unsorted dictionary (hash-map) + +.. dropdown:: Example: Container usage + :icon: code + + .. code-block:: c++ + + // BAD - Direct std containers + std::vector data; + std::array fixedData; + std::map orderedMap; + + // GOOD - GEOS wrappers + array1d data; + stackArray1d fixedData; + map orderedMap; + +The following standard library components are allowed: +- ``std::pair``, ``std::tuple`` for temporary structures, but beware of memory allocation, sometimes struct are to prefer. +- ``std::function`` for callbacks +- ``std::optional`` for optional return values +- ``std::move``, ``std::forward`` for move semantics +- Algorithms from ```` where appropriate +- ``std::unique_ptr``, ``std::shared_ptr`` for memory management + +Contribution +============ + +Documentation +------------- + +Always Provide Comprehensive Documentation +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +All new features must include appropriate documentation at multiple levels: + +Wrapper Documentation (User-Oriented) +"""""""""""""""""""""""""""""""""""""" + +All data repository wrappers must be documented with ``setDescription()``. +As much as possible, valid values rules should be provided, and default values should be provided with ``setApplyDefaultValue()``. + +.. dropdown:: Example: Wrapper documentation + :icon: code + + .. code-block:: c++ + + registerWrapper( "timeStep", &m_timeStep ). + setInputFlag( InputFlags::REQUIRED ). + setApplyDefaultValue( 1.0e-6 ). + setDescription( "Initial time step value. This parameter controls " + "the starting time increment for the simulation. " + "Valid range: (0, 1e-3]. See User Guide Chapter 5." ); + +RST Documentation (User Guide) +""""""""""""""""""""""""""""""" + +User-facing features must be documented in the Sphinx RST documentation (``src/docs/sphinx/``): + +- Explain **what** the feature does +- Provide **usage examples** +- Document **input parameters** +- Include **expected output** or behaviour +- Add **references** to related features + +Doxygen & Naming (Developer-Oriented) +"""""""""""""""""""""""""""""""""""""" + +Keep Doxygen comments and naming as clear as possible. Keep in mind that it is target to developers from other domains. + +.. dropdown:: Example: Doxygen documentation + :icon: code + + .. code-block:: c++ + // GOOD - Documentation that gives useful information about what the function *does*, without + // going into theorical knowledge nor implementation details. + /** + * @brief Computes the residual for the nonlinear solver. + * @param[in] domain The domain partition containing mesh and field data + * @param[in] time Current simulation time + * @param[in] dt Time step size + * @param[out] residual The computed residual vector + * @return Residual L2 norm + * + * This method assembles the discrete residual for the physics equations + * by looping over all elements in the target regions. The residual is + * stored in the provided vector and its L2 norm is returned. + */ + real64 computeResidual( DomainPartition & domain, + real64 const time, + real64 const dt, + arrayView1d const & residual ); + + // BAD - Documentation does not tell anything about the point of this function. + /** + * @brief Compute CFL numbers + * @param[in] time current time + * @param[in] dt the time step size + * @param[in] domain the domain partition + */ + void computeCFLNumbers( real64 const time, + real64 const dt, + DomainPartition & domain ) const; + + +Testing +------- + +Unit Tests +^^^^^^^^^^ + +**Always provide unit tests** to ensure code behaves according to expectations. +Unit testing is not about testing every single line of code (an unrealistic goal), but about *verifying assumptions and preventing regressions*. Don't assume your code works as intended—prove it with tests. +- Focus on **critical logic** (core functionality and algorithms), **edge cases** (extreme, erroneous, empty values). +- **Use GoogleTest framework**. +- In GEOS, the goal of unit test is never to take position on the validity of physical results (this is the point of integrated-tests). +- Place tests in ``src/coreComponents//tests/``. +- Test GPU code paths if applicable (use ``#ifdef GEOS_USE_DEVICE``). + +.. dropdown:: Example: Unit test structure + :icon: code + + .. code-block:: c++ + + #include + + // Test core functionality + TEST( MyComponentTest, ComputesCorrectResult ) + { + MyComponent component; + EXPECT_DOUBLE_EQ( component.compute( 2.0 ), 4.0 ); + } + + // Test edge case + TEST( MyComponentTest, HandlesZeroInput ) + { + MyComponent component; + EXPECT_DOUBLE_EQ( component.compute( 0.0 ), 0.0 ); + } + + // Test error detection + TEST( MyComponentTest, RejectsInvalidInput ) + { + MyComponent component; + EXPECT_THROW( component.compute( -1.0 ), std::invalid_argument ); + } + + // entry point of the test binary + int main( int argc, char * * argv ) + { + ::testing::InitGoogleTest( &argc, argv ); + g_commandLineOptions = *geos::basicSetup( argc, argv ); + int const result = RUN_ALL_TESTS(); + geos::basicCleanup(); + return result; + } + + +Code Coverage +^^^^^^^^^^^^^ + +**Code coverage should never decrease.** New code contributions must maintain or improve overall code coverage. +Use code-cov reports to identify untested code paths. + +Integrated Tests & examples +^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +**Every XML-usable & serializable component must have examples (at least, have an integrated-test case):** +- Each different usage pattern requires its own example, +- Examples must be ran as integrated tests to guaranty they keep valid, +- Place examples in ``examples/`` directory, and integrated tests in ``inputFiles/``, + +For more info, please refer to :ref:`IntegratedTests.rst`. + +Logging +======= + +Rules +----- + +Use GEOS Logging Macros +^^^^^^^^^^^^^^^^^^^^^^^ + +**Never use** ``std::cout``, ``std::cerr``, or ``printf`` for logging. Always use GEOS logging macros defined in ``Logger.hpp``: + +.. dropdown:: Example: Correct logging usage + :icon: code + + .. code-block:: c++ + + // BAD - Do not use + std::cout << "Starting simulation" << std::endl; + printf("Value = %f\n", value); + + // GOOD - Use GEOS macros + GEOS_LOG("Starting simulation"); + GEOS_LOG_RANK_0("Master rank initialization complete"); +**Available logging macros:** + +- ``GEOS_LOG(...)`` - Log on all ranks +- ``GEOS_LOG_VAR(...)`` - Log variable name and value +- ``GEOS_LOG_IF(condition, msg)`` - Conditional logging +- ``GEOS_LOG_RANK_0(msg)`` - Log only on rank 0 (MPI master) +- ``GEOS_LOG_RANK_0_IF(condition, msg)`` - Conditional logging on rank 0 +- ``GEOS_LOG_RANK(msg)`` - Log to rank-specific output stream +- ``GEOS_LOG_RANK_VAR(var)`` - Log variable on rank stream + +Use Log Levels +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +**Logging should be meaningful and controlled**: +- When possible, **use ``GEOS_LOG_LEVEL_RANK*`` macros with appropriate ``logInfo``** to allow output filtering, +- **Avoid unnecessary logging**, prefer rank 0 logging for global information to avoid redundant output, +- **Consider logs performance impact** (output flush frequency, formatting). + +See :doc:`LogLevels` for using / adding log level (e.g., ``logInfo::Convergence``, ``logInfo::TimeStep``). + +.. dropdown:: Example: Using log levels + :icon: code + + .. code-block:: c++ + + // physics solver cpp - effective logging + GEOS_LOG_LEVEL_RANK_0( logInfo::Convergence, + GEOS_FMT( "Newton iteration {}: residual = {:4.2e}", + iter, residual ) ); + + // src/coreComponents/physicsSolvers/LogLevelsInfo.hpp - convergence logging config + struct Convergence + { + static constexpr int getMinLogLevel() { return 1; } + static constexpr std::string_view getDescription() { return "Convergence information"; } + }; + +.. dropdown:: Example: Avoiding log spam + :icon: code + + .. code-block:: c++ + + // BAD - Will spam output for every element on every ranks + forAll< parallelDevicePolicy<> >( numElems, [=] GEOS_HOST_DEVICE ( localIndex const ei ) + { + GEOS_LOG("Processing element " << ei); // DO NOT DO THIS + }); + + // GOOD - Log summary information + GEOS_LOG_RANK_0(GEOS_FMT("Processed {} elements", numElems)); + +Logging in RAJA Kernels +^^^^^^^^^^^^^^^^^^^^^^^^ + +**Logging inside RAJA kernels only for DEBUG purposes** and must never appear on the develop branch. + +**Rationale:** Logging on GPU can: +- Reserve cache lines and degrade performance, +- Cause race conditions in output, +- Generate massive amounts of output, +- Produce unpredictable behaviour (e.g. runtime crashes on AMD). + +.. dropdown:: Example: Debug-only kernel logging + :icon: code + + .. code-block:: c++ + + forAll< parallelDevicePolicy<> >( numElems, [=] GEOS_HOST_DEVICE ( localIndex const ei ) + { + GEOS_LOG("processing..."); // Compiles, but remove before merge to develop! + GEOS_LOG(GEOS_FMT("processing elem {}/{}", ei, numElems)); // Will crash on CUDA builds (unsupported GEOS_FMT) + printf("DEBUG: ei=%d\n", ei); // Compiles, but remove before merge to develop! + }); + +Structured Data Output +---------------------- + +For tabular or statistical output, use structured logging facilities instead of plain text logging: + +- Use ``TableTextFormatter`` with a ``logLevel`` XML setting for formatted text table output, +- Use ``TableCSVFormatter`` with a ``writeCSV()`` XML setting for structured data output. + +.. dropdown:: Example: Log table output + :icon: code + + .. code-block:: c++ + // log output preparation, layout precomputation can be kept for periodic outputs + TableLayout statsLogLayout( "", { "Flux(es)", "Region", "Element Count", massColumn, rateColumn } ); + m_logLayout = statsLogLayout; + + // filling table data (log format is for current timestep) + m_logData.addRow( fluxName, + elementSetName, + GEOS_FMT( "{}", wrappedStats.stats().m_elementCount ), + GEOS_FMT( "{}", wrappedStats.stats().m_producedMass ), + GEOS_FMT( "{}", wrappedStats.stats().m_productionRate ) ); + + // statistics CSV output (only if enabled and from rank 0) + if( logLevelActive && MpiWrapper::commRank() == 0 ) + { + string const title = GEOS_FMT( "{}, flux statistics for: {}", + getName(), fluxesStr ); + string const multilineTitle = wrapTextToMaxLength( title, 80 ); // enabling automatic line wrapping + m_logLayout.setTitle( multilineTitle ); + TableTextFormatter const tableStatFormatter( m_logLayout ); + GEOS_LOG_RANK( tableStatFormatter.toString( statsData ) ); + } + +.. dropdown:: Example: CSV output + :icon: code + + .. code-block:: c++ + + // CSV output preparation, layout precomputation can be kept for periodic outputs + TableLayout const statsCSVLayout( "", {"Time [s]", "Flux(es)", "Region", "Element Count", massColumn, rateColumn} ); + m_csvLayout = statsCSVLayout; + + // filling table data (CSV format is all timestep series) + m_csvData.addRow( wrappedStats.getStatsPeriodStart(), + fluxName, + elementSetName, + wrappedStats.stats().m_elementCount, + wrappedStats.stats().m_producedMass, + wrappedStats.stats().m_productionRate ); + + // statistics CSV output (only if enabled and from rank 0) + if( m_writeCSV > 0 && MpiWrapper::commRank() == 0 ) + { + std::ofstream outputFile( m_csvFilename ); + TableCSVFormatter const tableStatFormatter( m_csvLayout ); + outputFile << tableStatFormatter.toString( csvData ); + outputFile.close(); + csvData.clear(); + } + +Error Management +================ + +- **Developpers must provide runtime validations** of their code behaviour, +- for runtime validation fails that can be encountered by users, **meaningful messages and context must be provided**. + +Errors (GEOS_ERROR*) +-------------------- + +Why Use Errors? +^^^^^^^^^^^^^^^ + +Use ``GEOS_ERROR*`` macros when encountering a **blocking error** where: + +- The simulation result would be certifiably invalid, +- The application state is unrecoverable, +- Continuing execution would be unsafe or meaningless. + +.. dropdown:: Example: Error handling + :icon: code + + .. code-block:: c++ + + // Example: not implemented feature + GEOS_ERROR( "Rounding interpolation with derivatives not implemented" ); + + // Example: invalid mesh topology + GEOS_ERROR_IF( numFaces == 0, + GEOS_FMT("Element {} has zero faces - invalid mesh", elemId) ); + + // Example: comparison-based error + GEOS_ERROR_IF_LT( dt, 0.0, "Time step must be positive" ); + +Exceptions (GEOS_THROW*) +------------------------ + +Why Use Exceptions? +^^^^^^^^^^^^^^^^^^^ + +Use ``GEOS_THROW*`` macros for the same reasons as ``GEOS_ERROR*`` (**unrecoverable state**), and when: + +- You want to **add context** at higher call stack levels, +- You need to **propagate detailed error information** upward +- If you use custom exception class, you need to document them here. + +.. list-table:: Available exception types + :header-rows: 1 + :widths: 35 65 + + * - Error class + - Used for + * - ``std::runtime_error`` + - General runtime errors + * - ``geos::InputError`` + - Errors in user input (XML, data files) + * - ``geos::SimulationError`` + - Errors during simulation execution + * - ``geos::BadTypeError`` + - Type conversion errors + * - ``geos::NotAnError`` + - Non critical, volontary stopping state + +.. dropdown:: Example: Throwing exceptions + :icon: code + + .. code-block:: c++ + + // Example 1: throw with context + GEOS_THROW_IF( !file.is_open(), + GEOS_FMT("Cannot open file {}", filename), + InputError ); + + // Example 2: rethrow with context info (because of unclear exception) + try + { + tolerance = stod( inputParams[8] ); + } + catch( const std::invalid_argument & e ) + { + GEOS_THROW( GEOS_FMT( "{}: invalid model parameter value: {}", functionName, e.what() ), + InputError, getDataContext() ); + } + +Never Catch and Continue +^^^^^^^^^^^^^^^^^^^^^^^^^ + +**Do not catch exceptions and continue execution**. Exceptions in GEOS indicate serious problems, they +do not serve as an alternative code path to consult a behaviour state. + +**Rationale:** + +- **Performance:** Exception handling is expensive, especially with stack unwinding, +- **Predictability:** Catching and continuing makes behaviour unpredictable and hard to debug, +- **Safety:** The error state may have corrupted simulation data, broke the code path at an unexpected place, and/or left the system in an invalid state. + +.. dropdown:: Example: Exception handling anti-pattern + :icon: code + + .. code-block:: c++ + + // BAD - Do not catch to continue + // Here, the "solver" should have a "isFailed" flag or a returned resulting state + try + { + solver.solve(); + } + catch( std::exception const & e ) // <- Note: this exception could have been throw for unexpected reason + { + GEOS_LOG("Solver failed, trying alternative"); + alternativeSolver.solve(); + } + + // GOOD - Manage your code paths manually + if(! solver.solve() ) + { + GEOS_LOG("Solver failed, trying alternative"); + alternativeSolver.solve(); + } + +Warnings (GEOS_WARNING*) +------------------------ + +Why Use Warnings? +^^^^^^^^^^^^^^^^^ + +Use ``GEOS_WARNING*`` macros when: + +- An issue is detected but **simulation can continue** +- Simulation may be **affected but not completly invalid** +- User should be notified of a **potential problem** + +.. dropdown:: Example: Warning usage + :icon: code + + .. code-block:: c++ + + // Example: suboptimal but valid configuration + GEOS_WARNING_IF( dt > recommendedDt, + GEOS_FMT("Time step {} exceeds recommended value {}. This may affect solution accuracy.", + dt, recommendedDt), + getDataContext() ); + +Errors & Warning in RAJA Kernels +-------------------------------- + +Same rule as logging: **Error and warning output inside RAJA kernels only for DEBUG purposes** and must never appear on the develop branch. + +- GPU kernel errors cause immediate kernel termination, adding potentially costly branches, +- Error handling on device has significant performance & cache impact, +- GPU error handling may be unsupported depending on the platform / device, +- Can cause deadlocks in parallel execution. + +.. dropdown:: Example: Kernel error handling (debug only) + :icon: code + + .. code-block:: c++ + + // On GPU, errors will cause kernel abort + forAll< parallelDevicePolicy<> >( numElems, [=] GEOS_HOST_DEVICE ( localIndex const ei ) + { + #if !defined( NDEBUG ) + // This will trap on GPU - use sparingly for debugging only + GEOS_ERROR_IF( badCondition, "Bad condition detected" ); + #endif + }); + +Contextualization with DataContext +----------------------------------- + +GEOS provides a abstract contextualization mechanism using ``DataContext`` to provide detailed error +information on the source of the data, which can be: + +- the file & line of the simulation XML file, +- if no source exists, ``Group`` / ``Wrapper`` path in the data repository. + +**Add ``DataContext`` when** Errors occur related to a **data repository component**: + +- a ``Wrapper`` when validation fails for a given user input setting, +- a ``Group`` when the error is implied by the instance state, +- multiple data-contexts when the error is due to multiple objects / settings (first = more important). + +Use ``getDataContext()`` or ``getWrapperDataContext()`` as last parameters of the logging macros to give context. + +.. dropdown:: Example: Using DataContext + :icon: code + + .. code-block:: c++ + + // Example 1, using the getWrapperDataContext() because the error is + // related to a precise input setting. + GEOS_ERROR_IF( !meshBodies.hasGroup( meshBodyName ), + GEOS_FMT( "MeshBody ({}) is specified in targetRegions, but does not exist.", + meshBodyName ), + getWrapperDataContext( viewKeyStruct::targetRegionsString() ) ); + + // Exemple 2, using getDataContext() because the error is related with the entire instance + // state (a physics solver for instance) + GEOS_ERROR_IF( convergenceFailed, + "Failed to converge!", + getDataContext() ); + + // Exemple 3, using multiple DataContext because the error is related with multiple objects + GEOS_THROW_IF( m_isThermal && !isFluidModelThermal, + GEOS_FMT( "CompositionalMultiphaseBase {}: the thermal option is enabled in the solver, but the fluid model {} is incompatible with the thermal option", + getDataContext(), fluid.getDataContext() ), + InputError, getDataContext(), fluid.getDataContext() ); + +Parallelism +=========== + +RAJA Kernels / CHAI Memory +--------------------------- + +Use RAJA Reductions +^^^^^^^^^^^^^^^^^^^ + +Use RAJA reductions for summation, min, max, and other parallel operations: + +.. dropdown:: Example: RAJA reductions + :icon: code + + .. code-block:: c++ + + RAJA::ReduceSum totalSum( 0.0 ); + RAJA::ReduceMin minValue( std::numeric_limits::max() ); + RAJA::ReduceMax maxValue( std::numeric_limits::lowest() ); + + forAll>( n, [=] GEOS_HOST_DEVICE ( localIndex i ) + { + totalSum += values[i]; + minValue.min( values[i] ); + maxValue.max( values[i] ); + }); + + real64 sum = totalSum.get(); + real64 min = minValue.get(); + real64 max = maxValue.get(); + +Do Not Launch Kernels from Within Kernels +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +**Never launch RAJA kernels from within other RAJA kernels.** This causes portability issues and is not supported on GPU. + +.. dropdown:: Example: Nested kernel anti-pattern + :icon: code + + .. code-block:: c++ + + // BAD - Do not nest kernels + forAll>( n, [=] GEOS_HOST_DEVICE ( localIndex i ) + { + forAll>( m, [=] GEOS_HOST_DEVICE ( localIndex j ) + { + // ... + }); + }); + + // GOOD - Use nested loops or restructure algorithm + forAll>( n, [=] GEOS_HOST_DEVICE ( localIndex i ) + { + for( localIndex j = 0; j < m; ++j ) + { + // ... + } + }); + +CHAI Memory Management +^^^^^^^^^^^^^^^^^^^^^^ + +GEOS uses CHAI for **automatic kernel memory migration**. Follow these rules: + +- **Use CHAI-managed arrays:** ``array1d``, ``array2d``, etc. +- **Call** ``.move(space)`` **to explicitly control memory location** +- **Never use manual CUDA calls** - RAJA has wrappers for device operations + +.. dropdown:: Example: CHAI memory management + :icon: code + + .. code-block:: c++ + + array1d data( 1000 ); + + // Explicitly move to device before kernel + // (would have been done automatically when giving an LvArray view to a RAJA kernel) + data.move( parallelDeviceMemorySpace ); + + // Use in kernel + auto dataView = data.toView(); + forAll>( numElems, [=] GEOS_HOST_DEVICE ( localIndex i ) + { + dataView[i] = computeValue( i ); + }); + + // Move back to host if needed for CPU operations (costly) + data.move( hostMemorySpace ); + +MPI Communication +----------------- + +Ensure All Ranks Participate in Collective Operations +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + + +**All MPI ranks must participate in collective operations.** Beware of branch deadlocks. + +.. dropdown:: Example: MPI collective operations + :icon: code + + .. code-block:: c++ + + // BAD - only rank 0 calls MPI_AllReduce() - DEADLOCK! + if( MpiWrapper::commRank() == 0 ) + { + MpiWrapper::allReduce( ... ); + } + + // BAD - some ranks may have this condition to false - not every ranks call MPI_Allreduce - DEADLOCK! + if( localWellElementsCount > 0 ) + { + MpiWrapper::allReduce( ... ); // WRONG - other ranks waiting! + } + + // GOOD - All ranks participate, even if they may not have influence in the result + real64 localSum = computeLocalSum(); + real64 globalSum = MpiWrapper::allReduce( localSum ); + + // GOOD - Conditional execution AFTER collective operation + real64 globalSum = MpiWrapper::allReduce( localSum ); + if( MpiWrapper::commRank() == 0 ) + { + GEOS_LOG_RANK_0( GEOS_FMT("Global sum = {}", globalSum) ); + } + +Always Use MpiWrapper +^^^^^^^^^^^^^^^^^^^^^ + +**Always use ``MpiWrapper`` instead of raw MPI calls.** +Refer to ``MpiWrapper.hpp`` for more. + +.. dropdown:: Example: Using MpiWrapper + :icon: code + + .. code-block:: c++ + + // BAD - Raw MPI calls + int rank; + MPI_Comm_rank( MPI_COMM_WORLD, &rank ); + + // GOOD - Use MpiWrapper + int rank = MpiWrapper::commRank(); + int size = MpiWrapper::commSize(); + MpiWrapper::barrier(); + + // Common patterns + real64 globalMax = MpiWrapper::max( localValue ); + real64 globalSum = MpiWrapper::sum( localValue ); + MpiWrapper::allGather( sendBuf, recvBuf, count ); + +Validation / Precision +====================== + +Use Tolerance-Based Comparisons +--------------------------------- + +**Always consider proper tolerance** for floating-point numbers comparisons, taking into account rounding errors, even for extreme values. + +.. dropdown:: Example: Correct float comparison + :icon: code + + .. code-block:: c++ + + // BAD - Direct comparison + if( value == 1.0 ) { ... } + if( a == b ) { ... } + + // GOOD - Absolute tolerance + real64 const absTol = 1.0e-12; + if( std::abs(value) < absTol ) { ... } + + // GOOD - Relative tolerance + real64 const relTol = 1.0e-8; + real64 const scale = std::max( std::abs(a), std::abs(b) ); + if( std::abs(a - b) < relTol * scale ) { ... } + +Use LvArray Math Utilities +^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The ``LvArray::math`` namespace provides safe numerical utilities, avoid LvArray math utilities instead of platform specific math function. + +.. dropdown:: Example: LvArray math utilities + :icon: code + + .. code-block:: c++ + + #include "LvArray/src/math.hpp" + + real64 const absValue = LvArray::math::abs( value ); + real64 const maxValue = LvArray::math::max( a, b ); + real64 const minValue = LvArray::math::min( a, b ); + +Validate in postInputInitialization +------------------------------------ + +**Perform cross-parameter validation as soon as possible**, and when possible **only once **. +This can be done in the ``postInputInitialization()`` method, or later if needed, (i.e. in the +``registerDataOnMesh()`` or during the simulation). + +At the ``postInputInitialization()`` stage: + +- All XML input has been parsed +- Default values have been applied +- Cross-references can be resolved +- User can receive clear error messages before simulation starts. + +.. dropdown:: Example: Post-input validation + :icon: code + + .. code-block:: c++ + + void MyClass::postInputInitialization() + { + // Call parent first + Base::postInputInitialization(); + + // Validate individual parameters + GEOS_THROW_IF_LT( m_timeStep, 0.0, + "timeStep must be positive", + InputError, getDataContext() ); + + // Validate cross-parameter constraints + GEOS_THROW_IF( m_endTime <= m_startTime, + "endTime must be greater than startTime", + InputError, getDataContext() ); + + // Validate against other components + GEOS_THROW_IF( !hasRequiredSolver(), + "Required solver not found in configuration", + InputError, getDataContext() ); + } + +Additional Validation Guidelines +------------------------------------ + +- Use ``setInputFlag()`` to mark parameters as ``REQUIRED`` or ``OPTIONAL`` +- Use ``setApplyDefaultValue()`` to provide sensible defaults +- Use ``setRTTypeName()`` for runtime type validation (e.g., ``rtTypes::CustomTypes::positive``) +- Document valid ranges in ``setDescription()`` + +Performance +=========== + +Profile When Affecting Kernel behaviour +--------------------------------------- + +When modifying performance-critical code (kernels, assembly loops, solvers): + +1. **Use Caliper profiling** on representative test cases before changes +2. **Measure performance** after changes +3. **Document performance impact** in pull request description + +.. dropdown:: Example: Performance profiling workflow + :icon: code + + .. code-block:: bash + + # Profile before changes + geos -i test_case.xml --caliper-output=baseline.cali + + # Make code changes + + # Profile after changes + geos -i test_case.xml --caliper-output=modified.cali + + # Compare results + cali-query -q "SELECT time.duration WHERE annotation='myKernel'" baseline.cali + cali-query -q "SELECT time.duration WHERE annotation='myKernel'" modified.cali + +Caliper Integration +------------------- + +Use ``GEOS_MARK_FUNCTION`` and ``Timer`` for performance tracking for the main computation functions / scopes. + +.. dropdown:: Example: Performance instrumentation + :icon: code + + .. code-block:: c++ + + void PhysicsSolverBase::solverStep( ... ) + { + GEOS_MARK_FUNCTION; // Caliper tracking for entire function + + { + GEOS_CALIPER_MARK_SCOPE("assemble system"); + Timer timer( m_timers.get_inserted( "assembly" ) ); + assembleSystem(); + } + + { + GEOS_CALIPER_MARK_SCOPE("linear solve"); + Timer timer( m_timers.get_inserted( "linear solver" ) ); + linearSolve(); + } + } + +Be Const / Constexpr +-------------------- + +**Use ``const`` and ``constexpr`` when possible**, enabling: + +- **Compiler optimization,** enables better code generation by giving constraints to the code, +- **Code safety,** prevents accidental modification for constant contexts, +- **Show code intention,** make code more readable. + +Memory Management Best Practices +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- **Mark methods const if the method is not designed to modify** the object state, +- **Minimize dynamic memory allocation** as much as possible, particularly in performance-critical code, +- For frequently allocated/deallocated objects, **consider memory pools**, +- For containers with known size, **reserve capacity upfront**. + +Value / Const Function Arguments +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +**Pass large and non-trivially copyable objects by const reference**: +- A "large" size can be defined as "more that 16 bytes, ", which is 2 pointers, 2 integer / index, or 2 double, +- An object is non-trivially copyable objects when it needs to perform sub-allocation when being copied, + +.. dropdown:: Value / Const Reference function parameters practices + :icon: code + + .. code-block:: c++ + + // 16 bytes (string_view are equivalent to 2 constant pointers): PASS BY VALUE + static constexpr string_view str0 = "Hello"; + string const str1 = getTimeStr(); // constant string are passed as string_view, beware of lifetime! + void func( string_view param ); + + arrayView2d< int > myDataView; + + // 16 bytes: PASS BY VALUE + stdArray< double, 2 > vec2D; + void func( vec2D param ); + + // 12 to 16 bytes (globalIndex is 8, integer is 4 to 8): PASS BY VALUE + struct SmallStruct + { + globalIndex id; + integer count; + }; + void func( SmallStruct param ); + + // 16 to 20 bytes (globalIndex is 8, integer is 4 to 8): PASS BY REFERENCE + struct BigStruct + { + globalIndex id; + integer countA; + integer countB; + }; + void func( BigStruct const & param ); + + // Does dynamic allocation: PASS BY REFERENCE + map< int, long > myMap; + stdVector< int > myList { 123 }; + void func( map< int, long > const & myMap, + stdVector< int > const & myList ); + + // LvArray types parameter practices depends on the intent + array2d< int > myConstantData; + array2d< int > myMutableData;, + // passing as non-view means that we will potencially resize the arrays: PASS BY REFERENCE + void arrayResizer( array2d< const int > & myConstantData + array2d< int > & myMutableData ); + // passing as view means that we may just affect the data, only if non-const template: PASS BY VALUE + void arrayProcess( arrayView2d< const int > myConstantData + arrayView2d< int > myMutableData ); + + // 16 bytes superficially... but does dynamic allocation: PASS BY REFERENCE + struct DynAllocStruct + { + string const a; + }; + void func( DynAllocStruct const & param ); + + // Any Group subclass is large: PASS BY REFERENCE + void modify( DomainPartition & domain ); + +Pointer / Reference Function Arguments +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +When passing nullable objects as parameter, pointers can be passed. + +**Pointer parameters must always be null-checked when dereferenced**. +If applying this rule produces repetitive code, consider using a const reference. + +This rule imply the following: +- **A reference may never be null**, +- In a given method, **``this`` may never be null**. + +Passing a pointer or reference parameter show the intention of modifying the underlying object. +The difference of the two practices is that passing a pointer means that we consider that the +underlying object can be null. + +Constexpr for Compile-Time Constants +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +**Use ``constexpr`` for values known at compile time**. + +The rule is not absolute: when the impact is not significative, and the code needs is getting really unclear, the rule can be ignored. + +.. dropdown:: Example: Constexpr usage + :icon: code + + .. code-block:: c++ + + // Compile-time constants + // A rule of thumb is that oftenly any time a value is constexpr, it can also be static. + static constexpr localIndex numDimensions = 3; + static constexpr real64 tolerance = 1.0e-10; + + // Constexpr functions (evaluated at compile time when possible) + constexpr localIndex computeArraySize( localIndex n ) + { return n * n + 2 * n + 1; } + +Provide Views to Arrays +------------------------ + +Provide Views to Arrays +^^^^^^^^^^^^^^^^^^^^^^^ + +- When possible, prefer provide views to arrays (to const data if possible). +- This **must** be done especially for RAJA kernels, and views must be **captured by value**. + +The rule is generalizable to ``string_view`` for strings, but not applicable in GPU context. + +Why Use Views? + +- **No memory allocation:** Views are lightweight references +- **Mutability correctness:** Can provide ``const`` read-only access to inner data +- **GPU compatibility:** Views work seamlessly on device + +.. dropdown:: Example: Views for arrays + :icon: code + + .. code-block:: c++ + + // BAD - Creates a copy + array1d< real64 > copy = originalArray; + + // GOOD - Creates a view (lightweight) + arrayView1d< real64 > view = originalArray.toView(); + + // GOOD - Const view for read-only access + arrayView1d< real64 const > constView = originalArray.toViewConst(); + +General Architecture +==================== + +Avoid Coupling +-------------- + +Minimize coupling between components when possible **without compromising memory efficiency or execution speed**. + +Principles +^^^^^^^^^^ +BESOIN DE VERIFIER CES REGLES +- **Loose coupling:** Components should depend on interfaces, not concrete implementations, +- **No circular dependancies:** Consider the existing dependancies, do not make components co-dependant (headers inclusion, packages referencing in ``CMakeLists.txt``) +- **Dependency injection:** Pass dependencies explicitly rather than accessing globals, relying on **lambda**, **templates** and **minimal interfaces**, +- **Performance exceptions:** Tight coupling is acceptable when required for performance + +.. dropdown:: Example: Reducing coupling + :icon: code + + .. code-block:: c++ + + // BAD - Tight coupling to concrete class / implementation + class SolverA + { + void registerSomeMeshData( VTKMesh & mesh ); + void solveHypre( HypreInterface * solver, HypreMatrix * matrix ); + }; + + // GOOD - Depends on interface + class SolverB + { + void registerSomeMeshData( MeshLevel & mesh ); // More general interface + void solve( LAInterface * solver, MatrixBase * matrix ); + }; + + // Performance-critical tight coupling + template< typename FIELD_T > + class Kernel + { + // Direct access to specific data layout for performance + void compute( FIELD_T const & data ); + }; + + template<> Kernel::compute( arrayView1d< double > & field ) + { /* template specialization */ } + + template<> Kernel::compute( arrayView2d< double > & field ); + { /* template specialization */ } + +Avoid Globals and Mutable State +-------------------------------- + +**Minimize mutable global states when possible.** +Prefer passing context explicitly. + +Why to avoid globals: + +- **Thread safety:** Globals can cause race conditions in parallel code +- **Testability:** Hard to test code that depends on global state +- **Predictability:** Global state makes code behaviour harder to understand +- **Encapsulation:** Violates modularity principles + +.. dropdown:: Example: Avoiding global state + :icon: code + + .. code-block:: c++ + + // BAD - Global mutable state + static real64 g_tolerance = 1.0e-6; + + void solve() + { + if( residual < g_tolerance ) { ... } // Depends on global + } + + // GOOD - Pass as parameter + void solve( real64 tolerance ) + { + if( residual < tolerance ) { ... } + } + + // GOOD - Member variable + class Solver + { + real64 const m_tolerance; + public: + Solver( real64 tolerance ) : m_tolerance( 1.0e-6 ) {} + + void solve() + { + if( residual < m_tolerance ) { ... } + } + }; + +Acceptable Global Usage: + +- **Library wrappers** (MPI wrapper, system-level resources, dependancies interface) +- **Read-only configuration** (immutable after initialization) +- **Performance counters** (for profiling) From eea776b98505246470e69baa5b4d352813de7ae1 Mon Sep 17 00:00:00 2001 From: MelReyCG Date: Wed, 26 Nov 2025 09:34:52 +0100 Subject: [PATCH 02/13] =?UTF-8?q?=F0=9F=93=9D=20LLM=20test?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../Contributing/CodeRules2.rst | 1346 +++++++++++++++++ 1 file changed, 1346 insertions(+) create mode 100644 src/docs/sphinx/developerGuide/Contributing/CodeRules2.rst diff --git a/src/docs/sphinx/developerGuide/Contributing/CodeRules2.rst b/src/docs/sphinx/developerGuide/Contributing/CodeRules2.rst new file mode 100644 index 00000000000..e183ad8e767 --- /dev/null +++ b/src/docs/sphinx/developerGuide/Contributing/CodeRules2.rst @@ -0,0 +1,1346 @@ +############################################################################### +Code Rules +############################################################################### + +.. contents:: Table of Contents + :depth: 3 + :local: + :backlinks: none + +============ +Introduction +============ + +This document outlines the coding rules and best practices that must be followed when contributing to GEOS. +These rules ensure code consistency, performance, maintainability, and reliability across the entire codebase, +especially considering our HPC context with CPU/GPU execution environments. + +All code must comply with these rules before being merged into the develop branch. Violations should be caught +during code review, and CI tests enforce many of these rules automatically. + +======= +Logging +======= + +Output Stream Rules +------------------- + +**Rule: No Direct Output Streams** + Never use ``std::cout``, ``std::cerr``, ``printf``, or any direct output functions in production code. + Always use the GEOS logging macros defined in ``Logger.hpp``. + + **Rationale:** Direct output streams bypass MPI rank control and log level filtering, leading to + uncontrolled output in parallel runs. + + .. code-block:: cpp + + // WRONG + std::cout << "Value is " << value << std::endl; + printf("Error: %s\n", message); + std::cerr << "Debug: " << variable << "\n"; + + // CORRECT + GEOS_LOG( "Value is " << value ); + GEOS_LOG_RANK_0( "Value is " << value ); // Only on rank 0 + GEOS_LOG_LEVEL( logInfo::Convergence, "Value is " << value ); + +**Rule: Avoid Log Spam** + Use appropriate log levels and conditional logging to avoid flooding output, especially in loops + and frequently called functions. + + **Rationale:** Excessive logging degrades performance and makes important messages hard to find. + + .. code-block:: cpp + + // WRONG - logs every iteration + for( localIndex iter = 0; iter < 1000; ++iter ) + { + GEOS_LOG( "Iteration " << iter ); + // ... + } + + // CORRECT - conditional logging + for( localIndex iter = 0; iter < 1000; ++iter ) + { + GEOS_LOG_RANK_0_IF( iter % 100 == 0, + GEOS_FMT( "Iteration {}: residual = {:.2e}", iter, residual ) ); + // ... + } + +GPU Kernel Logging Restrictions +-------------------------------- + +**Rule: No Logging in Device Code** + Never use logging macros in RAJA kernels or any code that may execute on GPU in production builds. + + **Rationale:** Logging in device code causes cache reservations, produces suboptimal code, and + will fail compilation on most GPU architectures. + + .. code-block:: cpp + + // WRONG - will fail on GPU + forAll< policy >( size, [=] GEOS_HOST_DEVICE ( localIndex const i ) + { + GEOS_LOG( "Processing element " << i ); // COMPILATION ERROR on GPU! + array[i] *= 2.0; + }); + + // CORRECT - log before/after kernel + GEOS_LOG_RANK_0( "Starting kernel execution for " << size << " elements" ); + forAll< policy >( size, [=] GEOS_HOST_DEVICE ( localIndex const i ) + { + array[i] *= 2.0; + }); + GEOS_LOG_RANK_0( "Kernel execution complete" ); + + // OK for temporary debugging only (remove before merge) + #ifdef DEBUG_KERNEL + forAll< policy >( size, [=] GEOS_HOST_DEVICE ( localIndex const i ) + { + if( i == debugIndex ) + printf( "Debug: element %d value %f\n", i, array[i] ); + array[i] *= 2.0; + }); + #endif + +Structured Data Output +---------------------- + +**Rule: Use Proper Formatting Tools** + For structured data output, use specialized formatting tools instead of manual formatting. + + **Rationale:** Ensures consistent, maintainable, and correctly aligned output. + + .. code-block:: cpp + + // WRONG - manual formatting + GEOS_LOG( "Iter Residual Time" ); + GEOS_LOG( iter << " " << residual << " " << time ); + + // CORRECT - use TableFormatter + TableFormatter table; + table.addColumn( "Iteration" ).addColumn( "Residual" ).addColumn( "Time [s]" ); + table.addRow( iter, residual, elapsedTime ); + GEOS_LOG( table.toString() ); + +================ +Error Management +================ + +Error Handling Hierarchy +------------------------ + +**Rule: Use Appropriate Error Macros** + Choose the correct error macro based on severity and recoverability. + + **Rationale:** Proper error handling enables graceful degradation and helps debugging. + + .. admonition:: Error Macro Selection Guide + + - ``GEOS_ERROR``: Unrecoverable errors requiring immediate termination + - ``GEOS_THROW``: Errors that should propagate with context + - ``GEOS_WARNING``: Non-critical issues that don't stop execution + - ``GEOS_ASSERT``: Debug-only checks (compiled out in release) + +Fatal Errors +------------ + +**Rule: Use GEOS_ERROR for Unrecoverable States** + Use ``GEOS_ERROR`` when the program cannot continue safely. + + **Rationale:** Immediate termination prevents data corruption and cascading failures. + + .. code-block:: cpp + + // CORRECT - unrecoverable dimension mismatch + GEOS_ERROR_IF( matrix.numRows() != rhs.size(), + GEOS_FMT( "Matrix-vector dimension mismatch: {} != {}", + matrix.numRows(), rhs.size() ) ); + + // CORRECT - using comparison macros + GEOS_ERROR_IF_NE_MSG( status, SUCCESS, + "Solver failed with status code " << status ); + + // CORRECT - null pointer check + GEOS_ERROR_IF( ptr == nullptr, + "Critical pointer is null in " << __FUNCTION__ ); + +Recoverable Errors +------------------ + +**Rule: Use GEOS_THROW for Propagating Errors** + Use ``GEOS_THROW`` when parent code might handle the error. + + **Rationale:** Allows error context to be added at different levels. + + .. code-block:: cpp + + // CORRECT - input validation + GEOS_THROW_IF( dt <= 0.0, + GEOS_FMT( "Invalid timestep: {} <= 0", dt ), + InputError ); + + // CORRECT - with data context + GEOS_THROW_IF( !materialModel, + "Material model not found", + std::runtime_error, + getWrapperDataContext( viewKeyStruct::constitutiveNameString() ) ); + + .. warning:: + **Never catch and recover from exceptions in normal flow.** Exceptions should only be used + for error propagation due to performance implications and GPU incompatibility. + +Non-Critical Issues +------------------- + +**Rule: Use GEOS_WARNING for Non-Blocking Issues** + Use ``GEOS_WARNING`` for issues that should be reported but don't require stopping execution. + + **Rationale:** Allows users to identify potential problems without disrupting computation. + + .. code-block:: cpp + + // CORRECT - performance warning + GEOS_WARNING_IF( numIterations > 0.8 * maxIterations, + GEOS_FMT( "Solver approaching iteration limit: {}/{}", + numIterations, maxIterations ) ); + + // CORRECT - deprecation notice + GEOS_WARNING( "Method computeOld() is deprecated, use computeNew() instead" ); + +===================== +Parallel/MPI Rules +===================== + +Collective Operation Consistency +--------------------------------- + +**Rule: All Ranks Must Call Collective Operations** + Every MPI collective operation must be called by all ranks in the communicator with compatible parameters. + + **Rationale:** Mismatched collective calls cause deadlocks or undefined behavior. + + .. code-block:: cpp + + // WRONG - conditional collective causes DEADLOCK + void processData( real64 const localValue ) + { + if( MpiWrapper::commRank() == 0 ) + { + real64 globalSum = MpiWrapper::sum( localValue ); // DEADLOCK! Only rank 0 calls + } + } + + // WRONG - different operations on different ranks + void inconsistentOperation( real64 const value ) + { + if( MpiWrapper::commRank() == 0 ) + real64 result = MpiWrapper::sum( value ); // sum on rank 0 + else + real64 result = MpiWrapper::max( value ); // max on others - DEADLOCK! + } + + // CORRECT - all ranks participate + void processDataCorrectly( real64 const localValue ) + { + // All ranks call sum + real64 globalSum = MpiWrapper::sum( localValue ); + + // Only rank 0 uses the result + if( MpiWrapper::commRank() == 0 ) + { + GEOS_LOG( "Global sum: " << globalSum ); + } + } + + // CORRECT - consistent operations with different data + void consistentOperation( array1d< real64 > const & data ) + { + real64 localMax = data.empty() ? 0.0 : LvArray::math::max( data.begin(), data.end() ); + real64 globalMax = MpiWrapper::max( localMax ); // All ranks call with their local value + } + +Communication Batching +---------------------- + +**Rule: Batch Small Communications** + Combine multiple small messages into larger ones to reduce communication overhead. + + **Rationale:** MPI has significant per-message overhead; batching improves performance. + + .. code-block:: cpp + + // WRONG - many small communications + void sendDataInefficiently( arrayView1d< real64 const > const & data, + int const destination ) + { + for( localIndex i = 0; i < data.size(); ++i ) + { + MpiWrapper::send( &data[i], 1, destination, tag + i ); // N messages! + } + } + + // CORRECT - single batched communication + void sendDataEfficiently( arrayView1d< real64 const > const & data, + int const destination ) + { + MpiWrapper::send( data.data(), data.size(), destination, tag ); // 1 message + } + + // CORRECT - pack multiple arrays + void sendMultipleArrays( arrayView1d< real64 const > const & array1, + arrayView1d< real64 const > const & array2 ) + { + buffer_type buffer; + buffer.pack( array1.size() ).pack( array1 ).pack( array2.size() ).pack( array2 ); + MpiWrapper::send( buffer.data(), buffer.size(), destination, tag ); + } + +Non-Blocking Communication +-------------------------- + +**Rule: Overlap Communication with Computation** + Use non-blocking communication to hide latency when possible. + + **Rationale:** Enables computation-communication overlap, improving parallel efficiency. + + .. code-block:: cpp + + // WRONG - sequential communication and computation + void sequentialGhostUpdate( DomainPartition & domain ) + { + sendReceiveGhosts( domain ); // Blocking communication + computeInteriorCells( domain ); // Computation waits + computeBoundaryCells( domain ); + } + + // CORRECT - overlapped communication and computation + void overlappedGhostUpdate( DomainPartition & domain ) + { + // Start non-blocking ghost exchange + MPI_Request ghostRequests[2]; + startGhostExchange( domain, ghostRequests ); + + // Compute interior while ghosts are in flight + computeInteriorCells( domain ); + + // Wait for ghosts before boundary computation + MpiWrapper::waitAll( 2, ghostRequests ); + + // Now safe to compute boundary cells + computeBoundaryCells( domain ); + } + +====================================== +Precision and Numerical Considerations +====================================== + +Floating-Point Comparisons +-------------------------- + +**Rule: Never Use Direct Equality for Floating-Point** + Always use tolerance-based comparisons for floating-point numbers. + + **Rationale:** Floating-point arithmetic has inherent rounding errors. + + .. code-block:: cpp + + // WRONG - exact comparison + if( value == 0.0 ) { } + if( a == b ) { } + + // CORRECT - absolute tolerance for zero comparison + constexpr real64 machinePrecision = std::numeric_limits< real64 >::epsilon(); + constexpr real64 tolerance = 100.0 * machinePrecision; + + if( LvArray::math::abs( value ) < tolerance ) + { + // value is effectively zero + } + + // CORRECT - relative tolerance for general comparison + bool areEqual( real64 const a, real64 const b ) + { + real64 const diff = LvArray::math::abs( a - b ); + real64 const magnitude = LvArray::math::max( LvArray::math::abs( a ), + LvArray::math::abs( b ) ); + return diff <= tolerance * magnitude; + } + + // CORRECT - combined absolute and relative tolerance + bool robustEqual( real64 const a, real64 const b, + real64 const relTol = 1e-10, + real64 const absTol = 1e-14 ) + { + return LvArray::math::abs( a - b ) <= relTol * LvArray::math::max( + LvArray::math::abs( a ), LvArray::math::abs( b ) ) + absTol; + } + +Division Safety +--------------- + +**Rule: Check Denominators Before Division** + Always verify denominators are not zero or near-zero before division. + + **Rationale:** Prevents NaN/Inf propagation and numerical instabilities. + + .. code-block:: cpp + + // WRONG - unprotected division + real64 const normalizedResidual = residual / initialResidual; + real64 const strainRate = velocityGradient / thickness; + + // CORRECT - protected division + real64 computeNormalizedResidual( real64 const residual, + real64 const initialResidual ) + { + if( initialResidual > machinePrecision ) + return residual / initialResidual; + else + return residual; // or return a flag indicating special case + } + + // CORRECT - with error reporting + real64 safeDivide( real64 const numerator, + real64 const denominator, + string const & context ) + { + GEOS_ERROR_IF( LvArray::math::abs( denominator ) < machinePrecision, + GEOS_FMT( "Division by zero in {}: denominator = {:.2e}", + context, denominator ) ); + return numerator / denominator; + } + +Overflow Prevention +------------------- + +**Rule: Check for Potential Overflow** + Validate that operations won't exceed type limits, especially for index calculations. + + **Rationale:** Overflow leads to undefined behavior and memory corruption. + + .. code-block:: cpp + + // WRONG - unchecked multiplication + globalIndex totalDofs = numElements * numNodesPerElement * numDofsPerNode; + + // CORRECT - overflow checking + globalIndex computeTotalDofs( globalIndex const numElements, + localIndex const numNodesPerElement, + localIndex const numDofsPerNode ) + { + globalIndex const maxValue = std::numeric_limits< globalIndex >::max(); + + // Check first multiplication + GEOS_ERROR_IF( numElements > maxValue / numNodesPerElement, + "DOF calculation would overflow: elements * nodes" ); + + globalIndex const temp = numElements * numNodesPerElement; + + // Check second multiplication + GEOS_ERROR_IF( temp > maxValue / numDofsPerNode, + "DOF calculation would overflow: total DOFs" ); + + return temp * numDofsPerNode; + } + + // CORRECT - capped increment for iterative methods + real64 applyIncrement( real64 const current, real64 const increment ) + { + constexpr real64 maxRelativeChange = 10.0; + + if( LvArray::math::abs( increment ) > maxRelativeChange * LvArray::math::abs( current ) ) + { + GEOS_WARNING( GEOS_FMT( "Large increment capped: {:.2e} -> {:.2e}", + increment, maxRelativeChange * current ) ); + return current + LvArray::math::signum( increment ) * maxRelativeChange * + LvArray::math::abs( current ); + } + return current + increment; + } + +Conditioning Monitoring +----------------------- + +**Rule: Monitor Numerical Conditioning** + Check and report conditioning issues in linear systems and algorithms. + + **Rationale:** Poor conditioning leads to inaccurate results and convergence failures. + + .. code-block:: cpp + + // CORRECT - condition number monitoring + class LinearSystemSolver + { + void checkConditioning( Matrix const & matrix ) + { + real64 const conditionNumber = estimateConditionNumber( matrix ); + + GEOS_WARNING_IF( conditionNumber > 1e10, + GEOS_FMT( "Poor conditioning detected: κ = {:.2e}", conditionNumber ) ); + + GEOS_ERROR_IF( conditionNumber > 1e16, + "Matrix is numerically singular" ); + + // Adjust solver based on conditioning + if( conditionNumber > 1e8 ) + { + m_solver.useRobustPreconditioner(); + m_solver.increaseTolerance( 10.0 ); + } + } + }; + +====================== +Performance Guidelines +====================== + +Loop Optimization +----------------- + +**Rule: Hoist Loop Invariants** + Move computations that don't change during iterations outside the loop. + + **Rationale:** Reduces redundant computations, especially critical in GPU kernels with limited registers. + + .. code-block:: cpp + + // WRONG - recomputing invariants + forAll< policy >( size, [=] GEOS_HOST_DEVICE ( localIndex const i ) + { + real64 const factor = viscosity * deltaTime / density; // Computed size times! + velocity[i] += factor * acceleration[i]; + }); + + // CORRECT - precompute invariants + real64 const factor = viscosity * deltaTime / density; // Computed once + forAll< policy >( size, [=] GEOS_HOST_DEVICE ( localIndex const i ) + { + velocity[i] += factor * acceleration[i]; + }); + + // CORRECT - complex invariant hoisting + void computeStress( arrayView3d< real64 > const & stress, + arrayView3d< real64 const > const & strain ) + { + // Hoist material property calculations + real64 const lambda = bulkModulus - 2.0 * shearModulus / 3.0; + real64 const twoMu = 2.0 * shearModulus; + + forAll< policy >( numElements, [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + // Only element-specific calculations in loop + real64 const volumetricStrain = strain[k][0][0] + strain[k][1][1] + strain[k][2][2]; + stress[k][0][0] = lambda * volumetricStrain + twoMu * strain[k][0][0]; + stress[k][1][1] = lambda * volumetricStrain + twoMu * strain[k][1][1]; + stress[k][2][2] = lambda * volumetricStrain + twoMu * strain[k][2][2]; + }); + } + +Memory Access Patterns +---------------------- + +**Rule: Optimize Memory Access for Cache and Coalescing** + Access memory sequentially and ensure coalesced access on GPUs. + + **Rationale:** Memory bandwidth is often the bottleneck; proper access patterns can improve performance 10x or more. + + .. code-block:: cpp + + // WRONG - strided access (poor cache usage) + void transposeNaive( arrayView2d< real64 > const & output, + arrayView2d< real64 const > const & input ) + { + forAll< policy >( input.size( 0 ), [=] GEOS_HOST_DEVICE ( localIndex const i ) + { + for( localIndex j = 0; j < input.size( 1 ); ++j ) + { + output[j][i] = input[i][j]; // Strided write! + } + }); + } + + // CORRECT - tiled transpose for better cache usage + void transposeTiled( arrayView2d< real64 > const & output, + arrayView2d< real64 const > const & input ) + { + constexpr localIndex tileSize = 32; + + forAll< policy >( input.size( 0 ) / tileSize, + [=] GEOS_HOST_DEVICE ( localIndex const tileI ) + { + for( localIndex tileJ = 0; tileJ < input.size( 1 ) / tileSize; ++tileJ ) + { + // Process tile with better locality + for( localIndex i = 0; i < tileSize; ++i ) + { + for( localIndex j = 0; j < tileSize; ++j ) + { + localIndex const globalI = tileI * tileSize + i; + localIndex const globalJ = tileJ * tileSize + j; + output[globalJ][globalI] = input[globalI][globalJ]; + } + } + } + }); + } + +Kernel Fusion +------------- + +**Rule: Combine Related Kernels** + Fuse multiple simple kernels to reduce memory traffic and launch overhead. + + **Rationale:** Memory bandwidth is precious; each kernel launch has overhead. + + .. code-block:: cpp + + // WRONG - multiple passes over data + void updateParticlesSeparate( arrayView1d< R1Tensor > const & velocity, + arrayView1d< R1Tensor > const & position, + arrayView1d< R1Tensor const > const & acceleration, + real64 const dt ) + { + // First kernel - update velocity + forAll< policy >( velocity.size(), [=] GEOS_HOST_DEVICE ( localIndex const i ) + { + velocity[i] += acceleration[i] * dt; + }); + + // Second kernel - update position (requires updated velocity) + forAll< policy >( position.size(), [=] GEOS_HOST_DEVICE ( localIndex const i ) + { + position[i] += velocity[i] * dt; + }); + } + + // CORRECT - single fused kernel + void updateParticlesFused( arrayView1d< R1Tensor > const & velocity, + arrayView1d< R1Tensor > const & position, + arrayView1d< R1Tensor const > const & acceleration, + real64 const dt ) + { + forAll< policy >( velocity.size(), [=] GEOS_HOST_DEVICE ( localIndex const i ) + { + velocity[i] += acceleration[i] * dt; + position[i] += velocity[i] * dt; // Single pass through memory + }); + } + +========== +Data Types +========== + +Type Consistency +---------------- + +**Rule: Use GEOS Type Aliases** + Always use GEOS-defined type aliases instead of native C++ types. + + **Rationale:** Ensures portability, consistency, and allows global type changes. + + .. list-table:: Type Correspondence + :header-rows: 1 + :widths: 25 25 50 + + * - GEOS Type + - Native Type + - Usage + * - ``real64`` + - ``double`` + - Default floating-point + * - ``real32`` + - ``float`` + - Single precision when memory-limited + * - ``integer`` + - ``int`` + - General integer values + * - ``localIndex`` + - ``GEOS_LOCALINDEX_TYPE`` + - Indexing within MPI rank + * - ``globalIndex`` + - ``GEOS_GLOBALINDEX_TYPE`` + - Global indexing across ranks + * - ``size_t`` + - ``std::size_t`` + - Size and count types + * - ``string`` + - ``std::string`` + - String data + * - ``string_view`` + - ``std::string_view`` + - Non-owning string references + + .. code-block:: cpp + + // WRONG - using native types + double pressure = 101325.0; + int numIterations = 0; + std::vector< double > values; + + // CORRECT - using GEOS types + real64 pressure = 101325.0; + integer numIterations = 0; + array1d< real64 > values; + +Container Usage +--------------- + +**Rule: Use GEOS Array Classes** + Never use STL containers for numerical data. Use GEOS array classes for GPU portability. + + **Rationale:** STL containers are not GPU-portable and lack necessary memory management features. + + .. code-block:: cpp + + // WRONG - STL containers + std::vector< double > forces; + double matrix[100][100]; + std::array< double, 3 > vector; + + // CORRECT - GEOS arrays + array1d< real64 > forces; + array2d< real64 > matrix( 100, 100 ); + R1Tensor vector = { 0.0, 0.0, 0.0 }; + + // CORRECT - views for kernel arguments + arrayView1d< real64 > const forcesView = forces.toView(); + arrayView2d< real64 > const matrixView = matrix.toView(); + + // CORRECT - stack arrays for small, fixed-size data + stackArray1d< real64, 10 > localValues; + +====================== +Physics-Specific Rules +====================== + +Unit Consistency +---------------- + +**Rule: Document and Verify Physical Units** + Always document units in function signatures and verify dimensional consistency. + + **Rationale:** Unit errors are a common source of bugs in physics simulations. + + .. code-block:: cpp + + /** + * @brief Compute pressure from equation of state + * @param[in] density fluid density [kg/m³] + * @param[in] bulkModulus bulk modulus [Pa] + * @param[in] referenceDensity reference density [kg/m³] + * @return pressure [Pa] + */ + GEOS_HOST_DEVICE + real64 computePressure( real64 const density, // [kg/m³] + real64 const bulkModulus, // [Pa] + real64 const referenceDensity ) // [kg/m³] + { + // Dimensional analysis: [Pa] = [Pa] * [dimensionless] + return bulkModulus * ( density / referenceDensity - 1.0 ); + } + + // CORRECT - unit conversion utilities + namespace units + { + constexpr real64 convertPressure( real64 const value, + PressureUnit const from, + PressureUnit const to ) + { + real64 const pascals = ( from == PressureUnit::PSI ) ? value * 6894.76 : value; + return ( to == PressureUnit::PSI ) ? pascals / 6894.76 : pascals; + } + } + +Physical Bounds Checking +------------------------ + +**Rule: Enforce Physical Constraints** + Validate that state variables remain within physically meaningful bounds. + + **Rationale:** Unphysical values indicate errors and can cause solver failures. + + .. code-block:: cpp + + class ConstitutiveModel + { + void updateState( real64 const pressure, + real64 const temperature, + real64 const porosity ) + { + // Check physical bounds + GEOS_ERROR_IF( pressure < -1e5, + GEOS_FMT( "Pressure {:.2e} Pa below cavitation limit", pressure ) ); + + GEOS_ERROR_IF( temperature < 0.0, + GEOS_FMT( "Absolute temperature {:.2f} K is negative", temperature ) ); + + GEOS_ERROR_IF( porosity < 0.0 || porosity > 1.0, + GEOS_FMT( "Porosity {:.3f} outside [0,1]", porosity ) ); + + // Warn about unusual but valid values + GEOS_WARNING_IF( temperature > 1000.0, + GEOS_FMT( "High temperature {:.0f} K may be outside model validity", + temperature ) ); + } + }; + +Conservation Verification +------------------------- + +**Rule: Verify Conservation Laws** + Check that mass, momentum, and energy are conserved within tolerance. + + **Rationale:** Conservation violations indicate bugs or numerical issues. + + .. code-block:: cpp + + class FlowSolver + { + void verifyMassConservation( DomainPartition const & domain ) + { + real64 const totalMassInitial = computeTotalMass( domain, 0 ); + real64 const totalMassCurrent = computeTotalMass( domain, currentTime ); + real64 const netFlux = computeCumulativeFlux( domain ); + + real64 const imbalance = LvArray::math::abs( + ( totalMassCurrent - totalMassInitial ) - netFlux ); + + real64 const relativeTolerance = 1e-10; + real64 const scale = LvArray::math::max( + LvArray::math::abs( totalMassInitial ), + LvArray::math::abs( netFlux ) ); + + GEOS_ERROR_IF( imbalance > relativeTolerance * scale, + GEOS_FMT( "Mass conservation violated: imbalance = {:.2e} kg ({:.2e}%)", + imbalance, 100.0 * imbalance / scale ) ); + } + }; + +===================== +Solver-Specific Rules +===================== + +Convergence Criteria +-------------------- + +**Rule: Implement Robust Convergence Checks** + Use both absolute and relative tolerances with appropriate scaling. + + **Rationale:** Single tolerance types fail for problems with varying scales. + + .. code-block:: cpp + + class ConvergenceCriteria + { + bool checkConvergence( real64 const residualNorm, + real64 const initialResidualNorm, + real64 const solutionNorm, + real64 const updateNorm ) const + { + // Absolute tolerance for near-zero problems + if( residualNorm < m_absoluteTolerance ) + { + GEOS_LOG_LEVEL( logInfo::Convergence, + "Converged: absolute tolerance" ); + return true; + } + + // Relative reduction from initial + if( initialResidualNorm > 0 && + residualNorm < m_relativeTolerance * initialResidualNorm ) + { + GEOS_LOG_LEVEL( logInfo::Convergence, + "Converged: relative reduction" ); + return true; + } + + // Solution-scaled tolerance for nonlinear problems + if( solutionNorm > 0 && + residualNorm < m_relativeTolerance * solutionNorm ) + { + GEOS_LOG_LEVEL( logInfo::Convergence, + "Converged: solution-scaled" ); + return true; + } + + // Update-based convergence for Newton methods + if( updateNorm < m_updateTolerance * solutionNorm ) + { + GEOS_LOG_LEVEL( logInfo::Convergence, + "Converged: small update" ); + return true; + } + + return false; + } + }; + +Linear Solver Configuration +--------------------------- + +**Rule: Adapt Solver to Problem Characteristics** + Choose solver and preconditioner based on matrix properties. + + **Rationale:** Optimal solver choice can improve performance by orders of magnitude. + + .. code-block:: cpp + + void configureLinearSolver( LinearSolver & solver, + localIndex const matrixSize, + real64 const estimatedConditionNumber, + bool const isSymmetric ) + { + // Small problems: direct solver + if( matrixSize < 1000 ) + { + solver.parameters().setDirect(); + GEOS_LOG_LEVEL( logInfo::LinearSolver, "Using direct solver for small system" ); + } + // Symmetric positive definite: CG + else if( isSymmetric && estimatedConditionNumber < 1e6 ) + { + solver.parameters().setSolver( LinearSolverParameters::SolverType::CG ); + solver.parameters().setPreconditioner( LinearSolverParameters::PreconditionerType::AMG ); + GEOS_LOG_LEVEL( logInfo::LinearSolver, "Using CG+AMG for SPD system" ); + } + // Ill-conditioned or nonsymmetric: GMRES + else + { + solver.parameters().setSolver( LinearSolverParameters::SolverType::GMRES ); + solver.parameters().setKrylovSubspaceSize( 50 ); + + if( estimatedConditionNumber > 1e10 ) + { + solver.parameters().setPreconditioner( LinearSolverParameters::PreconditionerType::ILU ); + solver.parameters().setILUFill( 2 ); + GEOS_LOG_LEVEL( logInfo::LinearSolver, "Using GMRES+ILU(2) for ill-conditioned system" ); + } + else + { + solver.parameters().setPreconditioner( LinearSolverParameters::PreconditionerType::AMG ); + GEOS_LOG_LEVEL( logInfo::LinearSolver, "Using GMRES+AMG for general system" ); + } + } + } + +======================== +Memory Management Rules +======================== + +View Lifetime Management +------------------------ + +**Rule: Never Outlive Parent Arrays** + Ensure array views don't outlive their parent arrays. + + **Rationale:** Dangling views cause segmentation faults and undefined behavior. + + .. code-block:: cpp + + // WRONG - returning view of local array + arrayView1d< real64 > getDanglingView() + { + array1d< real64 > tempArray( 100 ); + return tempArray.toView(); // DANGER: tempArray destroyed, view dangles! + } + + // WRONG - storing view of temporary + class DataHolder + { + arrayView1d< real64 > m_view; + + void setData() + { + array1d< real64 > temp( 100 ); + m_view = temp.toView(); // DANGER: temp destroyed at end of scope! + } + }; + + // CORRECT - return the array + array1d< real64 > getSafeArray() + { + array1d< real64 > result( 100 ); + // populate... + return result; // Move semantics ensure safety + } + + // CORRECT - store array, create view when needed + class SafeDataHolder + { + array1d< real64 > m_data; + + arrayView1d< real64 > getView() { return m_data.toView(); } + arrayView1d< real64 const > getViewConst() const { return m_data.toViewConst(); } + }; + +GPU Memory Management +--------------------- + +**Rule: Minimize Host-Device Transfers** + Keep data on the appropriate memory space and minimize transfers. + + **Rationale:** PCIe bandwidth is often the bottleneck in GPU applications. + + .. code-block:: cpp + + class GPUOptimizedSolver + { + // WRONG - repeated transfers + void inefficientCompute( array1d< real64 > & data ) + { + for( integer iter = 0; iter < 100; ++iter ) + { + data.move( LvArray::MemorySpace::cuda ); // Transfer to GPU + computeOnGPU( data.toView() ); + data.move( LvArray::MemorySpace::host ); // Transfer back + real64 norm = computeNorm( data ); // On CPU + } + } + + // CORRECT - minimize transfers + void efficientCompute( array1d< real64 > & data ) + { + // Single transfer to GPU + data.move( LvArray::MemorySpace::cuda ); + + arrayView1d< real64 > dataView = data.toView(); + array1d< real64 > gpuNorms( 100 ); + arrayView1d< real64 > normsView = gpuNorms.toView(); + + for( integer iter = 0; iter < 100; ++iter ) + { + computeOnGPU( dataView ); + computeNormOnGPU( dataView, normsView[iter] ); // Stay on GPU + } + + // Single transfer back + gpuNorms.move( LvArray::MemorySpace::host ); + // Process norms on host + } + }; + +==================== +Testing Requirements +==================== + +Numerical Testing +----------------- + +**Rule: Test Edge Cases and Known Solutions** + Include tests for boundary conditions, special cases, and problems with analytical solutions. + + **Rationale:** Ensures correctness and catches regressions. + + .. code-block:: cpp + + TEST( MatrixOperations, Inversion ) + { + // Test identity matrix + R2SymTensor identity = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 }; + R2SymTensor inverse; + real64 det = identity.invert( inverse ); + + EXPECT_NEAR( det, 1.0, 1e-14 ); + for( integer i = 0; i < 6; ++i ) + { + EXPECT_NEAR( inverse[i], identity[i], 1e-14 ); + } + + // Test singular matrix detection + R2SymTensor singular = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 }; + det = singular.invert( inverse ); + EXPECT_NEAR( det, 0.0, 1e-14 ); + + // Test numerical stability + R2SymTensor nearSingular = { 1.0, 1.0, 1.0 + 1e-10, 0.0, 0.0, 0.0 }; + det = nearSingular.invert( inverse ); + EXPECT_GT( LvArray::math::abs( det ), 0.0 ); + + // Verify A * A^-1 = I + R2SymTensor shouldBeIdentity; + shouldBeIdentity.Rij_eq_AikBkj( nearSingular, inverse ); + for( integer i = 0; i < 3; ++i ) + { + EXPECT_NEAR( shouldBeIdentity[i], 1.0, 1e-8 ); // Relaxed tolerance for ill-conditioned + } + } + +GPU Testing +----------- + +**Rule: Test CPU and GPU Paths Identically** + Ensure both execution paths produce identical results within tolerance. + + **Rationale:** Catches GPU-specific bugs and ensures portability. + + .. code-block:: cpp + + template< typename POLICY > + class KernelTest : public ::testing::Test + { + protected: + void testDotProduct() + { + constexpr localIndex size = 10000; + array1d< real64 > a( size ); + array1d< real64 > b( size ); + + // Initialize with deterministic values + for( localIndex i = 0; i < size; ++i ) + { + a[i] = std::sin( i * 0.1 ); + b[i] = std::cos( i * 0.1 ); + } + + // Compute on device + a.move( LvArray::MemorySpace::cuda, false ); + b.move( LvArray::MemorySpace::cuda, false ); + + arrayView1d< real64 const > aView = a.toViewConst(); + arrayView1d< real64 const > bView = b.toViewConst(); + + RAJA::ReduceSum< typename POLICY::ReducePolicy, real64 > dotProduct( 0.0 ); + + forAll< POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const i ) + { + dotProduct += aView[i] * bView[i]; + }); + + real64 const result = dotProduct.get(); + + // Verify against CPU reference + real64 reference = 0.0; + for( localIndex i = 0; i < size; ++i ) + { + reference += a[i] * b[i]; + } + + // GPU reduction may have different rounding + EXPECT_NEAR( result, reference, 1e-10 * LvArray::math::abs( reference ) ); + } + }; + + using TestTypes = ::testing::Types< serialPolicy, + #ifdef GEOS_USE_CUDA + cuda_policy, + #endif + #ifdef GEOS_USE_HIP + hip_policy, + #endif + openmp_policy >; + + TYPED_TEST_SUITE( KernelTest, TestTypes ); + + TYPED_TEST( KernelTest, DotProduct ) + { + this->testDotProduct(); + } + +========================= +Contribution Requirements +========================= + +Documentation Standards +----------------------- + +**Rule: Document All Public APIs** + Provide complete Doxygen documentation for all public interfaces. + + **Rationale:** Good documentation is essential for maintainability and usability. + + .. code-block:: cpp + + /** + * @brief Solve a nonlinear system using Newton's method + * @tparam MATRIX_TYPE type of the system matrix + * @param[in] jacobian the Jacobian matrix of the system + * @param[in,out] solution initial guess on input, solution on output + * @param[in] residualFunction functor to compute residual + * @param[in] parameters solver configuration parameters + * @return convergence information including iterations and final residual + * + * @note The solution vector is modified in place + * @warning The jacobian must be non-singular + * + * This function implements a damped Newton method with line search. + * The convergence criteria are specified in the parameters struct. + */ + template< typename MATRIX_TYPE > + ConvergenceInfo solveNonlinearSystem( MATRIX_TYPE const & jacobian, + arrayView1d< real64 > const & solution, + std::function< void( arrayView1d< real64 const >, + arrayView1d< real64 > ) > const & residualFunction, + NonlinearSolverParameters const & parameters ); + +Unit Test Coverage +------------------ + +**Rule: Test All Code Paths** + Ensure unit tests cover normal operation, edge cases, and error conditions. + + **Rationale:** Comprehensive testing prevents regressions and ensures reliability. + + .. code-block:: cpp + + class SolverTest : public ::testing::Test + { + protected: + void SetUp() override + { + // Initialize test fixtures + m_solver = std::make_unique< LinearSolver >(); + setupTestMatrix( m_matrix ); + setupTestRHS( m_rhs ); + } + + std::unique_ptr< LinearSolver > m_solver; + CRSMatrix< real64 > m_matrix; + array1d< real64 > m_rhs; + }; + + TEST_F( SolverTest, SolveIdentitySystem ) + { + // Test trivial case + m_matrix.setIdentity(); + array1d< real64 > solution = m_rhs; + + EXPECT_TRUE( m_solver->solve( m_matrix, solution ) ); + + for( localIndex i = 0; i < m_rhs.size(); ++i ) + { + EXPECT_NEAR( solution[i], m_rhs[i], 1e-12 ); + } + } + + TEST_F( SolverTest, HandleSingularMatrix ) + { + // Test error handling + m_matrix.zero(); // Singular matrix + array1d< real64 > solution = m_rhs; + + EXPECT_FALSE( m_solver->solve( m_matrix, solution ) ); + EXPECT_GT( m_solver->lastResidualNorm(), m_solver->tolerance() ); + } + +========================= +Additional Best Practices +========================= + +Const Correctness +----------------- + +**Rule: Use Const Wherever Possible** + Mark variables, parameters, and member functions const when they don't modify state. + + **Rationale:** Prevents bugs, enables optimizations, and documents intent. + + .. code-block:: cpp + + class DataProcessor + { + public: + // Const member function + real64 computeAverage() const // Promises not to modify object + { + return m_sum / m_count; + } + + // Const parameters for input + void processData( arrayView1d< real64 const > const & input, // Can't modify input + arrayView1d< real64 > const & output ) // Can modify output + { + // Const local variables + real64 const scale = computeScale( input ); + localIndex const size = input.size(); + + for( localIndex i = 0; i < size; ++i ) + { + output[i] = scale * input[i]; + } + } + + private: + real64 m_sum; + mutable real64 m_cachedAverage; // Can be modified in const functions + localIndex m_count; + }; + +Thread Safety +------------- + +**Rule: Avoid Global Mutable State** + Minimize shared mutable state and use proper synchronization when necessary. + + **Rationale:** Prevents data races and ensures correctness in parallel execution. + + .. code-block:: cpp + + // WRONG - global mutable state + namespace globals + { + real64 convergenceTolerance = 1e-6; // DANGER: race condition! + } + + // CORRECT - pass parameters explicitly + class Solver + { + struct Parameters + { + real64 convergenceTolerance = 1e-6; + integer maxIterations = 100; + }; + + void solve( Parameters const & params ); // Thread-safe + }; + + // CORRECT - use thread-local storage when necessary + namespace stats + { + thread_local integer numKernelCalls = 0; // Per-thread counter + + integer getTotalKernelCalls() + { + // Need synchronization to aggregate + return MpiWrapper::sum( numKernelCalls ); + } + } + +Header Dependencies +------------------- + +**Rule: Minimize Header Inclusions** + Use forward declarations when possible and include only necessary headers. + + **Rationale:** Reduces compilation time and prevents circular dependencies. + + .. code-block:: cpp + + // header file: Solver.hpp + #pragma once + + // WRONG - including full headers unnecessarily + #include "Matrix.hpp" // Full class definition + #include "Vector.hpp" // Full class definition + #include "Preconditioner.hpp" // Full class definition + + // CORRECT - forward declarations when possible + namespace geos + { + class Matrix; // Forward declaration sufficient + class Vector; // Forward declaration sufficient + class Preconditioner; // Forward declaration sufficient + + class Solver + { + public: + void solve( Matrix const & A, Vector & x ); + void setPreconditioner( Preconditioner * precond ); + + private: + Preconditioner * m_preconditioner = nullptr; // Pointer - forward decl OK + }; + } + + // Include full headers only in implementation file + From 25bf77914b584419059f7183fd3e1b570b606199 Mon Sep 17 00:00:00 2001 From: MelReyCG Date: Wed, 26 Nov 2025 16:35:15 +0100 Subject: [PATCH 03/13] =?UTF-8?q?=F0=9F=93=9D=20new=20draft=20version?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../developerGuide/Contributing/CodeRules.rst | 296 +++- .../Contributing/CodeRules2.rst | 1346 ----------------- 2 files changed, 244 insertions(+), 1398 deletions(-) delete mode 100644 src/docs/sphinx/developerGuide/Contributing/CodeRules2.rst diff --git a/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst b/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst index 84aed982b7e..3815c437d96 100644 --- a/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst +++ b/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst @@ -140,7 +140,7 @@ Use GEOS tensor types for geometric and mechanical calculations: ------------------------ It is possible to use ``std`` library components for data on host memory, for doing so, the rule is -to **only use GEOS container wrappers** instead of direct standard library containers. +to **only use GEOS ``std`` container wrappers** instead of direct standard library containers. * - Standard C++ Type - GEOS Type @@ -313,7 +313,6 @@ Unit testing is not about testing every single line of code (an unrealistic goal return result; } - Code Coverage ^^^^^^^^^^^^^ @@ -336,10 +335,11 @@ Logging Rules ----- -Use GEOS Logging Macros -^^^^^^^^^^^^^^^^^^^^^^^ +Use GEOS Logging Macros (``GEOS_LOG*``) +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -**Never use** ``std::cout``, ``std::cerr``, or ``printf`` for logging. Always use GEOS logging macros defined in ``Logger.hpp``: +**Use of ``std::cout``, ``std::cerr``, or ``printf`` for logging must never appear** on the develop branch. +Always use GEOS logging macros defined in ``Logger.hpp``: .. dropdown:: Example: Correct logging usage :icon: code @@ -353,6 +353,7 @@ Use GEOS Logging Macros // GOOD - Use GEOS macros GEOS_LOG("Starting simulation"); GEOS_LOG_RANK_0("Master rank initialization complete"); + **Available logging macros:** - ``GEOS_LOG(...)`` - Log on all ranks @@ -701,31 +702,6 @@ Parallelism RAJA Kernels / CHAI Memory --------------------------- -Use RAJA Reductions -^^^^^^^^^^^^^^^^^^^ - -Use RAJA reductions for summation, min, max, and other parallel operations: - -.. dropdown:: Example: RAJA reductions - :icon: code - - .. code-block:: c++ - - RAJA::ReduceSum totalSum( 0.0 ); - RAJA::ReduceMin minValue( std::numeric_limits::max() ); - RAJA::ReduceMax maxValue( std::numeric_limits::lowest() ); - - forAll>( n, [=] GEOS_HOST_DEVICE ( localIndex i ) - { - totalSum += values[i]; - minValue.min( values[i] ); - maxValue.max( values[i] ); - }); - - real64 sum = totalSum.get(); - real64 min = minValue.get(); - real64 max = maxValue.get(); - Do Not Launch Kernels from Within Kernels ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -784,13 +760,37 @@ GEOS uses CHAI for **automatic kernel memory migration**. Follow these rules: // Move back to host if needed for CPU operations (costly) data.move( hostMemorySpace ); +Use RAJA Reductions +^^^^^^^^^^^^^^^^^^^ + +Use RAJA reductions for summation, min, max, and other parallel operations: + +.. dropdown:: Example: RAJA reductions + :icon: code + + .. code-block:: c++ + + RAJA::ReduceSum totalSum( 0.0 ); + RAJA::ReduceMin minValue( std::numeric_limits::max() ); + RAJA::ReduceMax maxValue( std::numeric_limits::lowest() ); + + forAll>( n, [=] GEOS_HOST_DEVICE ( localIndex i ) + { + totalSum += values[i]; + minValue.min( values[i] ); + maxValue.max( values[i] ); + }); + + real64 sum = totalSum.get(); + real64 min = minValue.get(); + real64 max = maxValue.get(); + MPI Communication ----------------- Ensure All Ranks Participate in Collective Operations ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - **All MPI ranks must participate in collective operations.** Beware of branch deadlocks. .. dropdown:: Example: MPI collective operations @@ -846,11 +846,42 @@ Refer to ``MpiWrapper.hpp`` for more. real64 globalSum = MpiWrapper::sum( localValue ); MpiWrapper::allGather( sendBuf, recvBuf, count ); +Communication Batching +^^^^^^^^^^^^^^^^^^^^^^^ + +**Rule: Batch Small Communications** + Combine multiple small messages into larger ones to reduce communication overhead. + + **Rationale:** MPI has significant per-message overhead; batching improves performance. + + .. code-block:: cpp + + // WRONG - many small communications + void sendDataInefficiently( arrayView1d< real64 const > const & data, + int const destination ) + { + for( localIndex i = 0; i < data.size(); ++i ) + { + MpiWrapper::send( &data[i], 1, destination, tag + i ); // N messages! + } + } + + // CORRECT - single batched communication + void sendDataEfficiently( arrayView1d< real64 const > const & data, + int const destination ) + { + MpiWrapper::send( data.data(), data.size(), destination, tag ); // 1 message + } + + Validation / Precision ====================== +Computation validation +----------------------- + Use Tolerance-Based Comparisons ---------------------------------- +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ **Always consider proper tolerance** for floating-point numbers comparisons, taking into account rounding errors, even for extreme values. @@ -872,6 +903,50 @@ Use Tolerance-Based Comparisons real64 const scale = std::max( std::abs(a), std::abs(b) ); if( std::abs(a - b) < relTol * scale ) { ... } +Division Safety, NaN/Inf values +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +- **Always verify denominators are not zero or near-zero before division**. + +- In General, **we should not make any computation that could result in NaN/Inf values**. + +.. dropdown:: Example: Division Safety + :icon: code + + .. code-block:: cpp + + // WRONG - unprotected division + real64 const normalizedResidual = residual / initialResidual; + real64 const strainRate = velocityGradient / thickness; + + // CORRECT - protected division + real64 computeNormalizedResidual( real64 const residual, + real64 const initialResidual ) + { + if( initialResidual > machinePrecision ) + return residual / initialResidual; + else + return residual; // or return a flag indicating special case + } + + // CORRECT - with error reporting + real64 safeDivide( real64 const numerator, + real64 const denominator, + string const & context ) + { + GEOS_ERROR_IF( LvArray::math::abs( denominator ) < machinePrecision, + GEOS_FMT( "Division by zero in {}: denominator = {:.2e}", + context, denominator ) ); + return numerator / denominator; + } + +Overflow Prevention +^^^^^^^^^^^^^^^^^^^^ + +Overflow leads to undefined behavior and memory corruption. + +**Validate that operations won't exceed type limits**, especially for index calculations. + Use LvArray Math Utilities ^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -887,9 +962,12 @@ The ``LvArray::math`` namespace provides safe numerical utilities, avoid LvArray real64 const absValue = LvArray::math::abs( value ); real64 const maxValue = LvArray::math::max( a, b ); real64 const minValue = LvArray::math::min( a, b ); + +Data repository validation +------------------------------- Validate in postInputInitialization ------------------------------------- +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ **Perform cross-parameter validation as soon as possible**, and when possible **only once **. This can be done in the ``postInputInitialization()`` method, or later if needed, (i.e. in the @@ -929,7 +1007,7 @@ At the ``postInputInitialization()`` stage: } Additional Validation Guidelines ------------------------------------- +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - Use ``setInputFlag()`` to mark parameters as ``REQUIRED`` or ``OPTIONAL`` - Use ``setApplyDefaultValue()`` to provide sensible defaults @@ -939,8 +1017,11 @@ Additional Validation Guidelines Performance =========== +Profiling +--------- + Profile When Affecting Kernel behaviour ---------------------------------------- +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ When modifying performance-critical code (kernels, assembly loops, solvers): @@ -966,7 +1047,7 @@ When modifying performance-critical code (kernels, assembly loops, solvers): cali-query -q "SELECT time.duration WHERE annotation='myKernel'" modified.cali Caliper Integration -------------------- +^^^^^^^^^^^^^^^^^^^^ Use ``GEOS_MARK_FUNCTION`` and ``Timer`` for performance tracking for the main computation functions / scopes. @@ -992,8 +1073,26 @@ Use ``GEOS_MARK_FUNCTION`` and ``Timer`` for performance tracking for the main c } } +Speed Optimization Rules +----------------------------- + +- **Hoist Loop Invariants**, move computations that don't change during iterations outside the loop. +- When it does not critically affect the code architecture and clarity, **Fuse multiple related kernels to reduce memory traffic and launch overhead** (i.e., statistics kernels process all physics field at once). +- **Optimize Memory Access for Cache and Coalescing**. Access memory sequentially and ensure coalesced access, especially on GPUs. +- **Minimize Host-Device Transfers**. Keep data on the appropriate memory space and minimize transfers. + +Memory Management Rules +--------------------------------- + +Generalities +^^^^^^^^^^^^^^^^^^^^ + +- **Minimize dynamic memory allocation** as much as possible, particularly in performance-critical code, +- For frequently allocated/deallocated objects, **consider memory pools**, +- For containers with known size, **reserve capacity upfront**. + Be Const / Constexpr --------------------- +^^^^^^^^^^^^^^^^^^^^ **Use ``const`` and ``constexpr`` when possible**, enabling: @@ -1001,13 +1100,7 @@ Be Const / Constexpr - **Code safety,** prevents accidental modification for constant contexts, - **Show code intention,** make code more readable. -Memory Management Best Practices -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -- **Mark methods const if the method is not designed to modify** the object state, -- **Minimize dynamic memory allocation** as much as possible, particularly in performance-critical code, -- For frequently allocated/deallocated objects, **consider memory pools**, -- For containers with known size, **reserve capacity upfront**. +Also, **Mark methods const if the method is not designed to modify** the object state. Value / Const Function Arguments ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -1113,10 +1206,7 @@ The rule is not absolute: when the impact is not significative, and the code nee { return n * n + 2 * n + 1; } Provide Views to Arrays ------------------------- - -Provide Views to Arrays -^^^^^^^^^^^^^^^^^^^^^^^ +^^^^^^^^^^^^^^^^^^^^^^^^ - When possible, prefer provide views to arrays (to const data if possible). - This **must** be done especially for RAJA kernels, and views must be **captured by value**. @@ -1143,6 +1233,54 @@ Why Use Views? // GOOD - Const view for read-only access arrayView1d< real64 const > constView = originalArray.toViewConst(); +View Lifetime Management +------------------------ + +**Never Outlive Parent Arrays** +Dangling views cause segmentation faults and undefined behavior, that can be particularly hard to diagnose. +The rule is applicable to ``arrayView*`` and ``string_view``. + +.. dropdown:: Example: Lifetime Management + :icon: code + + .. code-block:: cpp + + // WRONG - returning view of local array + arrayView1d< real64 > getDanglingView() + { + array1d< real64 > tempArray( 100 ); + return tempArray.toView(); // DANGER: tempArray destroyed, view dangles! + } + + // WRONG - storing view of temporary + class DataHolder + { + arrayView1d< real64 > m_view; + + void setData() + { + array1d< real64 > temp( 100 ); + m_view = temp.toView(); // DANGER: temp destroyed at end of scope! + } + }; + + // CORRECT - return the array + array1d< real64 > getSafeArray() + { + array1d< real64 > result( 100 ); + // populate... + return result; // Move semantics ensure safety + } + + // CORRECT - store array, create view when needed + class SafeDataHolder + { + array1d< real64 > m_data; + + arrayView1d< real64 > getView() { return m_data.toView(); } + arrayView1d< real64 const > getViewConst() const { return m_data.toViewConst(); } + }; + General Architecture ==================== @@ -1153,11 +1291,12 @@ Minimize coupling between components when possible **without compromising memory Principles ^^^^^^^^^^ -BESOIN DE VERIFIER CES REGLES + - **Loose coupling:** Components should depend on interfaces, not concrete implementations, -- **No circular dependancies:** Consider the existing dependancies, do not make components co-dependant (headers inclusion, packages referencing in ``CMakeLists.txt``) -- **Dependency injection:** Pass dependencies explicitly rather than accessing globals, relying on **lambda**, **templates** and **minimal interfaces**, -- **Performance exceptions:** Tight coupling is acceptable when required for performance +- **No circular dependancies:** Consider the existing GEOS dependancies to not make components co-dependant (headers inclusion, packages referencing in ``CMakeLists.txt``, avoid tightly coupled objects), +- **Dependency injection:** Public components should receive their dependencies from external sources. Pass required dependencies using intermediate types instead of direct implementation types, using lambda, templates, and minimal interfaces, relying on **lambda**, **templates** and **minimal interfaces** (loose coupling, testability), +- **Performance exceptions:** Tight coupling is acceptable when required for performance, +- **Minimize header inclusions and dependancies**. .. dropdown:: Example: Reducing coupling :icon: code @@ -1242,3 +1381,56 @@ Acceptable Global Usage: - **Library wrappers** (MPI wrapper, system-level resources, dependancies interface) - **Read-only configuration** (immutable after initialization) - **Performance counters** (for profiling) + +Magic values, ``Group`` paths +------------------------------- + +**Avoid magic values**: +- **arbitrary values should not be written more than once**, define constants, consider using or extending ``PhysicsConstants.hpp`` / ``Units.hpp``, +- Because of architecture considerations, **Avoid using full raw ``Group`` path**, especially when the path has parts unrelated to the component consulting it -> prefer as relative paths as possible, +- ``Group`` / ``Wrapper`` name in the data repository is considered as a magic value, **use the ``getCatalogName()`` / ``getName()`` getters**, +- **Prefer to let appear the calculus of constants** rather than writing its value directly without explaination (constexpr has no runtime cost). + +====================== +Physics-Specific Rules +====================== + +Unit Consistency +---------------- + +**Verify physical units**, and **document any non-S.I. units** (variable name, comments). +**Non-S.I. units are only for internal computation use.** + +Physical Bounds Checking +------------------------ + +Unphysical values indicate errors and can cause solver failures. +**Validate that state variables remain within physically meaningful bounds**. + +If a value is not strictly disallowed but does not seem possible for the model, **show a warning to the user that he can disable**. + + + .. code-block:: cpp + + class ConstitutiveModel + { + void updateState( real64 const pressure, + real64 const temperature, + real64 const porosity ) + { + // Check physical bounds + GEOS_ERROR_IF( pressure < 1e-5, + GEOS_FMT( "Pressure {:.2e} Pa below cavitation limit", pressure ) ); + + GEOS_ERROR_IF( temperature < 0.0, + GEOS_FMT( "Absolute temperature {:.2f} K is negative", temperature ) ); + + GEOS_ERROR_IF( porosity < 0.0 || porosity > 1.0, + GEOS_FMT( "Porosity {:.3f} outside [0,1]", porosity ) ); + + // Warn about unusual but valid values, allow user to disable warning + GEOS_WARNING_IF( isHighTemperatureWarningEnabled && temperature > 1000.0, + GEOS_FMT( "High temperature {:.0f} K may be outside model validity", + temperature ) ); + } + }; diff --git a/src/docs/sphinx/developerGuide/Contributing/CodeRules2.rst b/src/docs/sphinx/developerGuide/Contributing/CodeRules2.rst deleted file mode 100644 index e183ad8e767..00000000000 --- a/src/docs/sphinx/developerGuide/Contributing/CodeRules2.rst +++ /dev/null @@ -1,1346 +0,0 @@ -############################################################################### -Code Rules -############################################################################### - -.. contents:: Table of Contents - :depth: 3 - :local: - :backlinks: none - -============ -Introduction -============ - -This document outlines the coding rules and best practices that must be followed when contributing to GEOS. -These rules ensure code consistency, performance, maintainability, and reliability across the entire codebase, -especially considering our HPC context with CPU/GPU execution environments. - -All code must comply with these rules before being merged into the develop branch. Violations should be caught -during code review, and CI tests enforce many of these rules automatically. - -======= -Logging -======= - -Output Stream Rules -------------------- - -**Rule: No Direct Output Streams** - Never use ``std::cout``, ``std::cerr``, ``printf``, or any direct output functions in production code. - Always use the GEOS logging macros defined in ``Logger.hpp``. - - **Rationale:** Direct output streams bypass MPI rank control and log level filtering, leading to - uncontrolled output in parallel runs. - - .. code-block:: cpp - - // WRONG - std::cout << "Value is " << value << std::endl; - printf("Error: %s\n", message); - std::cerr << "Debug: " << variable << "\n"; - - // CORRECT - GEOS_LOG( "Value is " << value ); - GEOS_LOG_RANK_0( "Value is " << value ); // Only on rank 0 - GEOS_LOG_LEVEL( logInfo::Convergence, "Value is " << value ); - -**Rule: Avoid Log Spam** - Use appropriate log levels and conditional logging to avoid flooding output, especially in loops - and frequently called functions. - - **Rationale:** Excessive logging degrades performance and makes important messages hard to find. - - .. code-block:: cpp - - // WRONG - logs every iteration - for( localIndex iter = 0; iter < 1000; ++iter ) - { - GEOS_LOG( "Iteration " << iter ); - // ... - } - - // CORRECT - conditional logging - for( localIndex iter = 0; iter < 1000; ++iter ) - { - GEOS_LOG_RANK_0_IF( iter % 100 == 0, - GEOS_FMT( "Iteration {}: residual = {:.2e}", iter, residual ) ); - // ... - } - -GPU Kernel Logging Restrictions --------------------------------- - -**Rule: No Logging in Device Code** - Never use logging macros in RAJA kernels or any code that may execute on GPU in production builds. - - **Rationale:** Logging in device code causes cache reservations, produces suboptimal code, and - will fail compilation on most GPU architectures. - - .. code-block:: cpp - - // WRONG - will fail on GPU - forAll< policy >( size, [=] GEOS_HOST_DEVICE ( localIndex const i ) - { - GEOS_LOG( "Processing element " << i ); // COMPILATION ERROR on GPU! - array[i] *= 2.0; - }); - - // CORRECT - log before/after kernel - GEOS_LOG_RANK_0( "Starting kernel execution for " << size << " elements" ); - forAll< policy >( size, [=] GEOS_HOST_DEVICE ( localIndex const i ) - { - array[i] *= 2.0; - }); - GEOS_LOG_RANK_0( "Kernel execution complete" ); - - // OK for temporary debugging only (remove before merge) - #ifdef DEBUG_KERNEL - forAll< policy >( size, [=] GEOS_HOST_DEVICE ( localIndex const i ) - { - if( i == debugIndex ) - printf( "Debug: element %d value %f\n", i, array[i] ); - array[i] *= 2.0; - }); - #endif - -Structured Data Output ----------------------- - -**Rule: Use Proper Formatting Tools** - For structured data output, use specialized formatting tools instead of manual formatting. - - **Rationale:** Ensures consistent, maintainable, and correctly aligned output. - - .. code-block:: cpp - - // WRONG - manual formatting - GEOS_LOG( "Iter Residual Time" ); - GEOS_LOG( iter << " " << residual << " " << time ); - - // CORRECT - use TableFormatter - TableFormatter table; - table.addColumn( "Iteration" ).addColumn( "Residual" ).addColumn( "Time [s]" ); - table.addRow( iter, residual, elapsedTime ); - GEOS_LOG( table.toString() ); - -================ -Error Management -================ - -Error Handling Hierarchy ------------------------- - -**Rule: Use Appropriate Error Macros** - Choose the correct error macro based on severity and recoverability. - - **Rationale:** Proper error handling enables graceful degradation and helps debugging. - - .. admonition:: Error Macro Selection Guide - - - ``GEOS_ERROR``: Unrecoverable errors requiring immediate termination - - ``GEOS_THROW``: Errors that should propagate with context - - ``GEOS_WARNING``: Non-critical issues that don't stop execution - - ``GEOS_ASSERT``: Debug-only checks (compiled out in release) - -Fatal Errors ------------- - -**Rule: Use GEOS_ERROR for Unrecoverable States** - Use ``GEOS_ERROR`` when the program cannot continue safely. - - **Rationale:** Immediate termination prevents data corruption and cascading failures. - - .. code-block:: cpp - - // CORRECT - unrecoverable dimension mismatch - GEOS_ERROR_IF( matrix.numRows() != rhs.size(), - GEOS_FMT( "Matrix-vector dimension mismatch: {} != {}", - matrix.numRows(), rhs.size() ) ); - - // CORRECT - using comparison macros - GEOS_ERROR_IF_NE_MSG( status, SUCCESS, - "Solver failed with status code " << status ); - - // CORRECT - null pointer check - GEOS_ERROR_IF( ptr == nullptr, - "Critical pointer is null in " << __FUNCTION__ ); - -Recoverable Errors ------------------- - -**Rule: Use GEOS_THROW for Propagating Errors** - Use ``GEOS_THROW`` when parent code might handle the error. - - **Rationale:** Allows error context to be added at different levels. - - .. code-block:: cpp - - // CORRECT - input validation - GEOS_THROW_IF( dt <= 0.0, - GEOS_FMT( "Invalid timestep: {} <= 0", dt ), - InputError ); - - // CORRECT - with data context - GEOS_THROW_IF( !materialModel, - "Material model not found", - std::runtime_error, - getWrapperDataContext( viewKeyStruct::constitutiveNameString() ) ); - - .. warning:: - **Never catch and recover from exceptions in normal flow.** Exceptions should only be used - for error propagation due to performance implications and GPU incompatibility. - -Non-Critical Issues -------------------- - -**Rule: Use GEOS_WARNING for Non-Blocking Issues** - Use ``GEOS_WARNING`` for issues that should be reported but don't require stopping execution. - - **Rationale:** Allows users to identify potential problems without disrupting computation. - - .. code-block:: cpp - - // CORRECT - performance warning - GEOS_WARNING_IF( numIterations > 0.8 * maxIterations, - GEOS_FMT( "Solver approaching iteration limit: {}/{}", - numIterations, maxIterations ) ); - - // CORRECT - deprecation notice - GEOS_WARNING( "Method computeOld() is deprecated, use computeNew() instead" ); - -===================== -Parallel/MPI Rules -===================== - -Collective Operation Consistency ---------------------------------- - -**Rule: All Ranks Must Call Collective Operations** - Every MPI collective operation must be called by all ranks in the communicator with compatible parameters. - - **Rationale:** Mismatched collective calls cause deadlocks or undefined behavior. - - .. code-block:: cpp - - // WRONG - conditional collective causes DEADLOCK - void processData( real64 const localValue ) - { - if( MpiWrapper::commRank() == 0 ) - { - real64 globalSum = MpiWrapper::sum( localValue ); // DEADLOCK! Only rank 0 calls - } - } - - // WRONG - different operations on different ranks - void inconsistentOperation( real64 const value ) - { - if( MpiWrapper::commRank() == 0 ) - real64 result = MpiWrapper::sum( value ); // sum on rank 0 - else - real64 result = MpiWrapper::max( value ); // max on others - DEADLOCK! - } - - // CORRECT - all ranks participate - void processDataCorrectly( real64 const localValue ) - { - // All ranks call sum - real64 globalSum = MpiWrapper::sum( localValue ); - - // Only rank 0 uses the result - if( MpiWrapper::commRank() == 0 ) - { - GEOS_LOG( "Global sum: " << globalSum ); - } - } - - // CORRECT - consistent operations with different data - void consistentOperation( array1d< real64 > const & data ) - { - real64 localMax = data.empty() ? 0.0 : LvArray::math::max( data.begin(), data.end() ); - real64 globalMax = MpiWrapper::max( localMax ); // All ranks call with their local value - } - -Communication Batching ----------------------- - -**Rule: Batch Small Communications** - Combine multiple small messages into larger ones to reduce communication overhead. - - **Rationale:** MPI has significant per-message overhead; batching improves performance. - - .. code-block:: cpp - - // WRONG - many small communications - void sendDataInefficiently( arrayView1d< real64 const > const & data, - int const destination ) - { - for( localIndex i = 0; i < data.size(); ++i ) - { - MpiWrapper::send( &data[i], 1, destination, tag + i ); // N messages! - } - } - - // CORRECT - single batched communication - void sendDataEfficiently( arrayView1d< real64 const > const & data, - int const destination ) - { - MpiWrapper::send( data.data(), data.size(), destination, tag ); // 1 message - } - - // CORRECT - pack multiple arrays - void sendMultipleArrays( arrayView1d< real64 const > const & array1, - arrayView1d< real64 const > const & array2 ) - { - buffer_type buffer; - buffer.pack( array1.size() ).pack( array1 ).pack( array2.size() ).pack( array2 ); - MpiWrapper::send( buffer.data(), buffer.size(), destination, tag ); - } - -Non-Blocking Communication --------------------------- - -**Rule: Overlap Communication with Computation** - Use non-blocking communication to hide latency when possible. - - **Rationale:** Enables computation-communication overlap, improving parallel efficiency. - - .. code-block:: cpp - - // WRONG - sequential communication and computation - void sequentialGhostUpdate( DomainPartition & domain ) - { - sendReceiveGhosts( domain ); // Blocking communication - computeInteriorCells( domain ); // Computation waits - computeBoundaryCells( domain ); - } - - // CORRECT - overlapped communication and computation - void overlappedGhostUpdate( DomainPartition & domain ) - { - // Start non-blocking ghost exchange - MPI_Request ghostRequests[2]; - startGhostExchange( domain, ghostRequests ); - - // Compute interior while ghosts are in flight - computeInteriorCells( domain ); - - // Wait for ghosts before boundary computation - MpiWrapper::waitAll( 2, ghostRequests ); - - // Now safe to compute boundary cells - computeBoundaryCells( domain ); - } - -====================================== -Precision and Numerical Considerations -====================================== - -Floating-Point Comparisons --------------------------- - -**Rule: Never Use Direct Equality for Floating-Point** - Always use tolerance-based comparisons for floating-point numbers. - - **Rationale:** Floating-point arithmetic has inherent rounding errors. - - .. code-block:: cpp - - // WRONG - exact comparison - if( value == 0.0 ) { } - if( a == b ) { } - - // CORRECT - absolute tolerance for zero comparison - constexpr real64 machinePrecision = std::numeric_limits< real64 >::epsilon(); - constexpr real64 tolerance = 100.0 * machinePrecision; - - if( LvArray::math::abs( value ) < tolerance ) - { - // value is effectively zero - } - - // CORRECT - relative tolerance for general comparison - bool areEqual( real64 const a, real64 const b ) - { - real64 const diff = LvArray::math::abs( a - b ); - real64 const magnitude = LvArray::math::max( LvArray::math::abs( a ), - LvArray::math::abs( b ) ); - return diff <= tolerance * magnitude; - } - - // CORRECT - combined absolute and relative tolerance - bool robustEqual( real64 const a, real64 const b, - real64 const relTol = 1e-10, - real64 const absTol = 1e-14 ) - { - return LvArray::math::abs( a - b ) <= relTol * LvArray::math::max( - LvArray::math::abs( a ), LvArray::math::abs( b ) ) + absTol; - } - -Division Safety ---------------- - -**Rule: Check Denominators Before Division** - Always verify denominators are not zero or near-zero before division. - - **Rationale:** Prevents NaN/Inf propagation and numerical instabilities. - - .. code-block:: cpp - - // WRONG - unprotected division - real64 const normalizedResidual = residual / initialResidual; - real64 const strainRate = velocityGradient / thickness; - - // CORRECT - protected division - real64 computeNormalizedResidual( real64 const residual, - real64 const initialResidual ) - { - if( initialResidual > machinePrecision ) - return residual / initialResidual; - else - return residual; // or return a flag indicating special case - } - - // CORRECT - with error reporting - real64 safeDivide( real64 const numerator, - real64 const denominator, - string const & context ) - { - GEOS_ERROR_IF( LvArray::math::abs( denominator ) < machinePrecision, - GEOS_FMT( "Division by zero in {}: denominator = {:.2e}", - context, denominator ) ); - return numerator / denominator; - } - -Overflow Prevention -------------------- - -**Rule: Check for Potential Overflow** - Validate that operations won't exceed type limits, especially for index calculations. - - **Rationale:** Overflow leads to undefined behavior and memory corruption. - - .. code-block:: cpp - - // WRONG - unchecked multiplication - globalIndex totalDofs = numElements * numNodesPerElement * numDofsPerNode; - - // CORRECT - overflow checking - globalIndex computeTotalDofs( globalIndex const numElements, - localIndex const numNodesPerElement, - localIndex const numDofsPerNode ) - { - globalIndex const maxValue = std::numeric_limits< globalIndex >::max(); - - // Check first multiplication - GEOS_ERROR_IF( numElements > maxValue / numNodesPerElement, - "DOF calculation would overflow: elements * nodes" ); - - globalIndex const temp = numElements * numNodesPerElement; - - // Check second multiplication - GEOS_ERROR_IF( temp > maxValue / numDofsPerNode, - "DOF calculation would overflow: total DOFs" ); - - return temp * numDofsPerNode; - } - - // CORRECT - capped increment for iterative methods - real64 applyIncrement( real64 const current, real64 const increment ) - { - constexpr real64 maxRelativeChange = 10.0; - - if( LvArray::math::abs( increment ) > maxRelativeChange * LvArray::math::abs( current ) ) - { - GEOS_WARNING( GEOS_FMT( "Large increment capped: {:.2e} -> {:.2e}", - increment, maxRelativeChange * current ) ); - return current + LvArray::math::signum( increment ) * maxRelativeChange * - LvArray::math::abs( current ); - } - return current + increment; - } - -Conditioning Monitoring ------------------------ - -**Rule: Monitor Numerical Conditioning** - Check and report conditioning issues in linear systems and algorithms. - - **Rationale:** Poor conditioning leads to inaccurate results and convergence failures. - - .. code-block:: cpp - - // CORRECT - condition number monitoring - class LinearSystemSolver - { - void checkConditioning( Matrix const & matrix ) - { - real64 const conditionNumber = estimateConditionNumber( matrix ); - - GEOS_WARNING_IF( conditionNumber > 1e10, - GEOS_FMT( "Poor conditioning detected: κ = {:.2e}", conditionNumber ) ); - - GEOS_ERROR_IF( conditionNumber > 1e16, - "Matrix is numerically singular" ); - - // Adjust solver based on conditioning - if( conditionNumber > 1e8 ) - { - m_solver.useRobustPreconditioner(); - m_solver.increaseTolerance( 10.0 ); - } - } - }; - -====================== -Performance Guidelines -====================== - -Loop Optimization ------------------ - -**Rule: Hoist Loop Invariants** - Move computations that don't change during iterations outside the loop. - - **Rationale:** Reduces redundant computations, especially critical in GPU kernels with limited registers. - - .. code-block:: cpp - - // WRONG - recomputing invariants - forAll< policy >( size, [=] GEOS_HOST_DEVICE ( localIndex const i ) - { - real64 const factor = viscosity * deltaTime / density; // Computed size times! - velocity[i] += factor * acceleration[i]; - }); - - // CORRECT - precompute invariants - real64 const factor = viscosity * deltaTime / density; // Computed once - forAll< policy >( size, [=] GEOS_HOST_DEVICE ( localIndex const i ) - { - velocity[i] += factor * acceleration[i]; - }); - - // CORRECT - complex invariant hoisting - void computeStress( arrayView3d< real64 > const & stress, - arrayView3d< real64 const > const & strain ) - { - // Hoist material property calculations - real64 const lambda = bulkModulus - 2.0 * shearModulus / 3.0; - real64 const twoMu = 2.0 * shearModulus; - - forAll< policy >( numElements, [=] GEOS_HOST_DEVICE ( localIndex const k ) - { - // Only element-specific calculations in loop - real64 const volumetricStrain = strain[k][0][0] + strain[k][1][1] + strain[k][2][2]; - stress[k][0][0] = lambda * volumetricStrain + twoMu * strain[k][0][0]; - stress[k][1][1] = lambda * volumetricStrain + twoMu * strain[k][1][1]; - stress[k][2][2] = lambda * volumetricStrain + twoMu * strain[k][2][2]; - }); - } - -Memory Access Patterns ----------------------- - -**Rule: Optimize Memory Access for Cache and Coalescing** - Access memory sequentially and ensure coalesced access on GPUs. - - **Rationale:** Memory bandwidth is often the bottleneck; proper access patterns can improve performance 10x or more. - - .. code-block:: cpp - - // WRONG - strided access (poor cache usage) - void transposeNaive( arrayView2d< real64 > const & output, - arrayView2d< real64 const > const & input ) - { - forAll< policy >( input.size( 0 ), [=] GEOS_HOST_DEVICE ( localIndex const i ) - { - for( localIndex j = 0; j < input.size( 1 ); ++j ) - { - output[j][i] = input[i][j]; // Strided write! - } - }); - } - - // CORRECT - tiled transpose for better cache usage - void transposeTiled( arrayView2d< real64 > const & output, - arrayView2d< real64 const > const & input ) - { - constexpr localIndex tileSize = 32; - - forAll< policy >( input.size( 0 ) / tileSize, - [=] GEOS_HOST_DEVICE ( localIndex const tileI ) - { - for( localIndex tileJ = 0; tileJ < input.size( 1 ) / tileSize; ++tileJ ) - { - // Process tile with better locality - for( localIndex i = 0; i < tileSize; ++i ) - { - for( localIndex j = 0; j < tileSize; ++j ) - { - localIndex const globalI = tileI * tileSize + i; - localIndex const globalJ = tileJ * tileSize + j; - output[globalJ][globalI] = input[globalI][globalJ]; - } - } - } - }); - } - -Kernel Fusion -------------- - -**Rule: Combine Related Kernels** - Fuse multiple simple kernels to reduce memory traffic and launch overhead. - - **Rationale:** Memory bandwidth is precious; each kernel launch has overhead. - - .. code-block:: cpp - - // WRONG - multiple passes over data - void updateParticlesSeparate( arrayView1d< R1Tensor > const & velocity, - arrayView1d< R1Tensor > const & position, - arrayView1d< R1Tensor const > const & acceleration, - real64 const dt ) - { - // First kernel - update velocity - forAll< policy >( velocity.size(), [=] GEOS_HOST_DEVICE ( localIndex const i ) - { - velocity[i] += acceleration[i] * dt; - }); - - // Second kernel - update position (requires updated velocity) - forAll< policy >( position.size(), [=] GEOS_HOST_DEVICE ( localIndex const i ) - { - position[i] += velocity[i] * dt; - }); - } - - // CORRECT - single fused kernel - void updateParticlesFused( arrayView1d< R1Tensor > const & velocity, - arrayView1d< R1Tensor > const & position, - arrayView1d< R1Tensor const > const & acceleration, - real64 const dt ) - { - forAll< policy >( velocity.size(), [=] GEOS_HOST_DEVICE ( localIndex const i ) - { - velocity[i] += acceleration[i] * dt; - position[i] += velocity[i] * dt; // Single pass through memory - }); - } - -========== -Data Types -========== - -Type Consistency ----------------- - -**Rule: Use GEOS Type Aliases** - Always use GEOS-defined type aliases instead of native C++ types. - - **Rationale:** Ensures portability, consistency, and allows global type changes. - - .. list-table:: Type Correspondence - :header-rows: 1 - :widths: 25 25 50 - - * - GEOS Type - - Native Type - - Usage - * - ``real64`` - - ``double`` - - Default floating-point - * - ``real32`` - - ``float`` - - Single precision when memory-limited - * - ``integer`` - - ``int`` - - General integer values - * - ``localIndex`` - - ``GEOS_LOCALINDEX_TYPE`` - - Indexing within MPI rank - * - ``globalIndex`` - - ``GEOS_GLOBALINDEX_TYPE`` - - Global indexing across ranks - * - ``size_t`` - - ``std::size_t`` - - Size and count types - * - ``string`` - - ``std::string`` - - String data - * - ``string_view`` - - ``std::string_view`` - - Non-owning string references - - .. code-block:: cpp - - // WRONG - using native types - double pressure = 101325.0; - int numIterations = 0; - std::vector< double > values; - - // CORRECT - using GEOS types - real64 pressure = 101325.0; - integer numIterations = 0; - array1d< real64 > values; - -Container Usage ---------------- - -**Rule: Use GEOS Array Classes** - Never use STL containers for numerical data. Use GEOS array classes for GPU portability. - - **Rationale:** STL containers are not GPU-portable and lack necessary memory management features. - - .. code-block:: cpp - - // WRONG - STL containers - std::vector< double > forces; - double matrix[100][100]; - std::array< double, 3 > vector; - - // CORRECT - GEOS arrays - array1d< real64 > forces; - array2d< real64 > matrix( 100, 100 ); - R1Tensor vector = { 0.0, 0.0, 0.0 }; - - // CORRECT - views for kernel arguments - arrayView1d< real64 > const forcesView = forces.toView(); - arrayView2d< real64 > const matrixView = matrix.toView(); - - // CORRECT - stack arrays for small, fixed-size data - stackArray1d< real64, 10 > localValues; - -====================== -Physics-Specific Rules -====================== - -Unit Consistency ----------------- - -**Rule: Document and Verify Physical Units** - Always document units in function signatures and verify dimensional consistency. - - **Rationale:** Unit errors are a common source of bugs in physics simulations. - - .. code-block:: cpp - - /** - * @brief Compute pressure from equation of state - * @param[in] density fluid density [kg/m³] - * @param[in] bulkModulus bulk modulus [Pa] - * @param[in] referenceDensity reference density [kg/m³] - * @return pressure [Pa] - */ - GEOS_HOST_DEVICE - real64 computePressure( real64 const density, // [kg/m³] - real64 const bulkModulus, // [Pa] - real64 const referenceDensity ) // [kg/m³] - { - // Dimensional analysis: [Pa] = [Pa] * [dimensionless] - return bulkModulus * ( density / referenceDensity - 1.0 ); - } - - // CORRECT - unit conversion utilities - namespace units - { - constexpr real64 convertPressure( real64 const value, - PressureUnit const from, - PressureUnit const to ) - { - real64 const pascals = ( from == PressureUnit::PSI ) ? value * 6894.76 : value; - return ( to == PressureUnit::PSI ) ? pascals / 6894.76 : pascals; - } - } - -Physical Bounds Checking ------------------------- - -**Rule: Enforce Physical Constraints** - Validate that state variables remain within physically meaningful bounds. - - **Rationale:** Unphysical values indicate errors and can cause solver failures. - - .. code-block:: cpp - - class ConstitutiveModel - { - void updateState( real64 const pressure, - real64 const temperature, - real64 const porosity ) - { - // Check physical bounds - GEOS_ERROR_IF( pressure < -1e5, - GEOS_FMT( "Pressure {:.2e} Pa below cavitation limit", pressure ) ); - - GEOS_ERROR_IF( temperature < 0.0, - GEOS_FMT( "Absolute temperature {:.2f} K is negative", temperature ) ); - - GEOS_ERROR_IF( porosity < 0.0 || porosity > 1.0, - GEOS_FMT( "Porosity {:.3f} outside [0,1]", porosity ) ); - - // Warn about unusual but valid values - GEOS_WARNING_IF( temperature > 1000.0, - GEOS_FMT( "High temperature {:.0f} K may be outside model validity", - temperature ) ); - } - }; - -Conservation Verification -------------------------- - -**Rule: Verify Conservation Laws** - Check that mass, momentum, and energy are conserved within tolerance. - - **Rationale:** Conservation violations indicate bugs or numerical issues. - - .. code-block:: cpp - - class FlowSolver - { - void verifyMassConservation( DomainPartition const & domain ) - { - real64 const totalMassInitial = computeTotalMass( domain, 0 ); - real64 const totalMassCurrent = computeTotalMass( domain, currentTime ); - real64 const netFlux = computeCumulativeFlux( domain ); - - real64 const imbalance = LvArray::math::abs( - ( totalMassCurrent - totalMassInitial ) - netFlux ); - - real64 const relativeTolerance = 1e-10; - real64 const scale = LvArray::math::max( - LvArray::math::abs( totalMassInitial ), - LvArray::math::abs( netFlux ) ); - - GEOS_ERROR_IF( imbalance > relativeTolerance * scale, - GEOS_FMT( "Mass conservation violated: imbalance = {:.2e} kg ({:.2e}%)", - imbalance, 100.0 * imbalance / scale ) ); - } - }; - -===================== -Solver-Specific Rules -===================== - -Convergence Criteria --------------------- - -**Rule: Implement Robust Convergence Checks** - Use both absolute and relative tolerances with appropriate scaling. - - **Rationale:** Single tolerance types fail for problems with varying scales. - - .. code-block:: cpp - - class ConvergenceCriteria - { - bool checkConvergence( real64 const residualNorm, - real64 const initialResidualNorm, - real64 const solutionNorm, - real64 const updateNorm ) const - { - // Absolute tolerance for near-zero problems - if( residualNorm < m_absoluteTolerance ) - { - GEOS_LOG_LEVEL( logInfo::Convergence, - "Converged: absolute tolerance" ); - return true; - } - - // Relative reduction from initial - if( initialResidualNorm > 0 && - residualNorm < m_relativeTolerance * initialResidualNorm ) - { - GEOS_LOG_LEVEL( logInfo::Convergence, - "Converged: relative reduction" ); - return true; - } - - // Solution-scaled tolerance for nonlinear problems - if( solutionNorm > 0 && - residualNorm < m_relativeTolerance * solutionNorm ) - { - GEOS_LOG_LEVEL( logInfo::Convergence, - "Converged: solution-scaled" ); - return true; - } - - // Update-based convergence for Newton methods - if( updateNorm < m_updateTolerance * solutionNorm ) - { - GEOS_LOG_LEVEL( logInfo::Convergence, - "Converged: small update" ); - return true; - } - - return false; - } - }; - -Linear Solver Configuration ---------------------------- - -**Rule: Adapt Solver to Problem Characteristics** - Choose solver and preconditioner based on matrix properties. - - **Rationale:** Optimal solver choice can improve performance by orders of magnitude. - - .. code-block:: cpp - - void configureLinearSolver( LinearSolver & solver, - localIndex const matrixSize, - real64 const estimatedConditionNumber, - bool const isSymmetric ) - { - // Small problems: direct solver - if( matrixSize < 1000 ) - { - solver.parameters().setDirect(); - GEOS_LOG_LEVEL( logInfo::LinearSolver, "Using direct solver for small system" ); - } - // Symmetric positive definite: CG - else if( isSymmetric && estimatedConditionNumber < 1e6 ) - { - solver.parameters().setSolver( LinearSolverParameters::SolverType::CG ); - solver.parameters().setPreconditioner( LinearSolverParameters::PreconditionerType::AMG ); - GEOS_LOG_LEVEL( logInfo::LinearSolver, "Using CG+AMG for SPD system" ); - } - // Ill-conditioned or nonsymmetric: GMRES - else - { - solver.parameters().setSolver( LinearSolverParameters::SolverType::GMRES ); - solver.parameters().setKrylovSubspaceSize( 50 ); - - if( estimatedConditionNumber > 1e10 ) - { - solver.parameters().setPreconditioner( LinearSolverParameters::PreconditionerType::ILU ); - solver.parameters().setILUFill( 2 ); - GEOS_LOG_LEVEL( logInfo::LinearSolver, "Using GMRES+ILU(2) for ill-conditioned system" ); - } - else - { - solver.parameters().setPreconditioner( LinearSolverParameters::PreconditionerType::AMG ); - GEOS_LOG_LEVEL( logInfo::LinearSolver, "Using GMRES+AMG for general system" ); - } - } - } - -======================== -Memory Management Rules -======================== - -View Lifetime Management ------------------------- - -**Rule: Never Outlive Parent Arrays** - Ensure array views don't outlive their parent arrays. - - **Rationale:** Dangling views cause segmentation faults and undefined behavior. - - .. code-block:: cpp - - // WRONG - returning view of local array - arrayView1d< real64 > getDanglingView() - { - array1d< real64 > tempArray( 100 ); - return tempArray.toView(); // DANGER: tempArray destroyed, view dangles! - } - - // WRONG - storing view of temporary - class DataHolder - { - arrayView1d< real64 > m_view; - - void setData() - { - array1d< real64 > temp( 100 ); - m_view = temp.toView(); // DANGER: temp destroyed at end of scope! - } - }; - - // CORRECT - return the array - array1d< real64 > getSafeArray() - { - array1d< real64 > result( 100 ); - // populate... - return result; // Move semantics ensure safety - } - - // CORRECT - store array, create view when needed - class SafeDataHolder - { - array1d< real64 > m_data; - - arrayView1d< real64 > getView() { return m_data.toView(); } - arrayView1d< real64 const > getViewConst() const { return m_data.toViewConst(); } - }; - -GPU Memory Management ---------------------- - -**Rule: Minimize Host-Device Transfers** - Keep data on the appropriate memory space and minimize transfers. - - **Rationale:** PCIe bandwidth is often the bottleneck in GPU applications. - - .. code-block:: cpp - - class GPUOptimizedSolver - { - // WRONG - repeated transfers - void inefficientCompute( array1d< real64 > & data ) - { - for( integer iter = 0; iter < 100; ++iter ) - { - data.move( LvArray::MemorySpace::cuda ); // Transfer to GPU - computeOnGPU( data.toView() ); - data.move( LvArray::MemorySpace::host ); // Transfer back - real64 norm = computeNorm( data ); // On CPU - } - } - - // CORRECT - minimize transfers - void efficientCompute( array1d< real64 > & data ) - { - // Single transfer to GPU - data.move( LvArray::MemorySpace::cuda ); - - arrayView1d< real64 > dataView = data.toView(); - array1d< real64 > gpuNorms( 100 ); - arrayView1d< real64 > normsView = gpuNorms.toView(); - - for( integer iter = 0; iter < 100; ++iter ) - { - computeOnGPU( dataView ); - computeNormOnGPU( dataView, normsView[iter] ); // Stay on GPU - } - - // Single transfer back - gpuNorms.move( LvArray::MemorySpace::host ); - // Process norms on host - } - }; - -==================== -Testing Requirements -==================== - -Numerical Testing ------------------ - -**Rule: Test Edge Cases and Known Solutions** - Include tests for boundary conditions, special cases, and problems with analytical solutions. - - **Rationale:** Ensures correctness and catches regressions. - - .. code-block:: cpp - - TEST( MatrixOperations, Inversion ) - { - // Test identity matrix - R2SymTensor identity = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 }; - R2SymTensor inverse; - real64 det = identity.invert( inverse ); - - EXPECT_NEAR( det, 1.0, 1e-14 ); - for( integer i = 0; i < 6; ++i ) - { - EXPECT_NEAR( inverse[i], identity[i], 1e-14 ); - } - - // Test singular matrix detection - R2SymTensor singular = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 }; - det = singular.invert( inverse ); - EXPECT_NEAR( det, 0.0, 1e-14 ); - - // Test numerical stability - R2SymTensor nearSingular = { 1.0, 1.0, 1.0 + 1e-10, 0.0, 0.0, 0.0 }; - det = nearSingular.invert( inverse ); - EXPECT_GT( LvArray::math::abs( det ), 0.0 ); - - // Verify A * A^-1 = I - R2SymTensor shouldBeIdentity; - shouldBeIdentity.Rij_eq_AikBkj( nearSingular, inverse ); - for( integer i = 0; i < 3; ++i ) - { - EXPECT_NEAR( shouldBeIdentity[i], 1.0, 1e-8 ); // Relaxed tolerance for ill-conditioned - } - } - -GPU Testing ------------ - -**Rule: Test CPU and GPU Paths Identically** - Ensure both execution paths produce identical results within tolerance. - - **Rationale:** Catches GPU-specific bugs and ensures portability. - - .. code-block:: cpp - - template< typename POLICY > - class KernelTest : public ::testing::Test - { - protected: - void testDotProduct() - { - constexpr localIndex size = 10000; - array1d< real64 > a( size ); - array1d< real64 > b( size ); - - // Initialize with deterministic values - for( localIndex i = 0; i < size; ++i ) - { - a[i] = std::sin( i * 0.1 ); - b[i] = std::cos( i * 0.1 ); - } - - // Compute on device - a.move( LvArray::MemorySpace::cuda, false ); - b.move( LvArray::MemorySpace::cuda, false ); - - arrayView1d< real64 const > aView = a.toViewConst(); - arrayView1d< real64 const > bView = b.toViewConst(); - - RAJA::ReduceSum< typename POLICY::ReducePolicy, real64 > dotProduct( 0.0 ); - - forAll< POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const i ) - { - dotProduct += aView[i] * bView[i]; - }); - - real64 const result = dotProduct.get(); - - // Verify against CPU reference - real64 reference = 0.0; - for( localIndex i = 0; i < size; ++i ) - { - reference += a[i] * b[i]; - } - - // GPU reduction may have different rounding - EXPECT_NEAR( result, reference, 1e-10 * LvArray::math::abs( reference ) ); - } - }; - - using TestTypes = ::testing::Types< serialPolicy, - #ifdef GEOS_USE_CUDA - cuda_policy, - #endif - #ifdef GEOS_USE_HIP - hip_policy, - #endif - openmp_policy >; - - TYPED_TEST_SUITE( KernelTest, TestTypes ); - - TYPED_TEST( KernelTest, DotProduct ) - { - this->testDotProduct(); - } - -========================= -Contribution Requirements -========================= - -Documentation Standards ------------------------ - -**Rule: Document All Public APIs** - Provide complete Doxygen documentation for all public interfaces. - - **Rationale:** Good documentation is essential for maintainability and usability. - - .. code-block:: cpp - - /** - * @brief Solve a nonlinear system using Newton's method - * @tparam MATRIX_TYPE type of the system matrix - * @param[in] jacobian the Jacobian matrix of the system - * @param[in,out] solution initial guess on input, solution on output - * @param[in] residualFunction functor to compute residual - * @param[in] parameters solver configuration parameters - * @return convergence information including iterations and final residual - * - * @note The solution vector is modified in place - * @warning The jacobian must be non-singular - * - * This function implements a damped Newton method with line search. - * The convergence criteria are specified in the parameters struct. - */ - template< typename MATRIX_TYPE > - ConvergenceInfo solveNonlinearSystem( MATRIX_TYPE const & jacobian, - arrayView1d< real64 > const & solution, - std::function< void( arrayView1d< real64 const >, - arrayView1d< real64 > ) > const & residualFunction, - NonlinearSolverParameters const & parameters ); - -Unit Test Coverage ------------------- - -**Rule: Test All Code Paths** - Ensure unit tests cover normal operation, edge cases, and error conditions. - - **Rationale:** Comprehensive testing prevents regressions and ensures reliability. - - .. code-block:: cpp - - class SolverTest : public ::testing::Test - { - protected: - void SetUp() override - { - // Initialize test fixtures - m_solver = std::make_unique< LinearSolver >(); - setupTestMatrix( m_matrix ); - setupTestRHS( m_rhs ); - } - - std::unique_ptr< LinearSolver > m_solver; - CRSMatrix< real64 > m_matrix; - array1d< real64 > m_rhs; - }; - - TEST_F( SolverTest, SolveIdentitySystem ) - { - // Test trivial case - m_matrix.setIdentity(); - array1d< real64 > solution = m_rhs; - - EXPECT_TRUE( m_solver->solve( m_matrix, solution ) ); - - for( localIndex i = 0; i < m_rhs.size(); ++i ) - { - EXPECT_NEAR( solution[i], m_rhs[i], 1e-12 ); - } - } - - TEST_F( SolverTest, HandleSingularMatrix ) - { - // Test error handling - m_matrix.zero(); // Singular matrix - array1d< real64 > solution = m_rhs; - - EXPECT_FALSE( m_solver->solve( m_matrix, solution ) ); - EXPECT_GT( m_solver->lastResidualNorm(), m_solver->tolerance() ); - } - -========================= -Additional Best Practices -========================= - -Const Correctness ------------------ - -**Rule: Use Const Wherever Possible** - Mark variables, parameters, and member functions const when they don't modify state. - - **Rationale:** Prevents bugs, enables optimizations, and documents intent. - - .. code-block:: cpp - - class DataProcessor - { - public: - // Const member function - real64 computeAverage() const // Promises not to modify object - { - return m_sum / m_count; - } - - // Const parameters for input - void processData( arrayView1d< real64 const > const & input, // Can't modify input - arrayView1d< real64 > const & output ) // Can modify output - { - // Const local variables - real64 const scale = computeScale( input ); - localIndex const size = input.size(); - - for( localIndex i = 0; i < size; ++i ) - { - output[i] = scale * input[i]; - } - } - - private: - real64 m_sum; - mutable real64 m_cachedAverage; // Can be modified in const functions - localIndex m_count; - }; - -Thread Safety -------------- - -**Rule: Avoid Global Mutable State** - Minimize shared mutable state and use proper synchronization when necessary. - - **Rationale:** Prevents data races and ensures correctness in parallel execution. - - .. code-block:: cpp - - // WRONG - global mutable state - namespace globals - { - real64 convergenceTolerance = 1e-6; // DANGER: race condition! - } - - // CORRECT - pass parameters explicitly - class Solver - { - struct Parameters - { - real64 convergenceTolerance = 1e-6; - integer maxIterations = 100; - }; - - void solve( Parameters const & params ); // Thread-safe - }; - - // CORRECT - use thread-local storage when necessary - namespace stats - { - thread_local integer numKernelCalls = 0; // Per-thread counter - - integer getTotalKernelCalls() - { - // Need synchronization to aggregate - return MpiWrapper::sum( numKernelCalls ); - } - } - -Header Dependencies -------------------- - -**Rule: Minimize Header Inclusions** - Use forward declarations when possible and include only necessary headers. - - **Rationale:** Reduces compilation time and prevents circular dependencies. - - .. code-block:: cpp - - // header file: Solver.hpp - #pragma once - - // WRONG - including full headers unnecessarily - #include "Matrix.hpp" // Full class definition - #include "Vector.hpp" // Full class definition - #include "Preconditioner.hpp" // Full class definition - - // CORRECT - forward declarations when possible - namespace geos - { - class Matrix; // Forward declaration sufficient - class Vector; // Forward declaration sufficient - class Preconditioner; // Forward declaration sufficient - - class Solver - { - public: - void solve( Matrix const & A, Vector & x ); - void setPreconditioner( Preconditioner * precond ); - - private: - Preconditioner * m_preconditioner = nullptr; // Pointer - forward decl OK - }; - } - - // Include full headers only in implementation file - From d95a006affef1591a42e1fa60228ac32b1c0a187 Mon Sep 17 00:00:00 2001 From: MelReyCG Date: Fri, 28 Nov 2025 10:11:24 +0100 Subject: [PATCH 04/13] =?UTF-8?q?=F0=9F=93=9Dprecised=20usage=20of=20conta?= =?UTF-8?q?iner=20(Arnaud=20review)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../sphinx/developerGuide/Contributing/CodeRules.rst | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst b/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst index 3815c437d96..9e1a69a59ab 100644 --- a/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst +++ b/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst @@ -110,13 +110,15 @@ Refer to the rule of :ref:`CHAI Memory Management`. .. code-block:: c++ - // BAD - Do not use std::vector + // BAD - Do not use std::vector for data that can be packed or kernel addressable std::vector values; + std::array fixedData; - // GOOD - Use GEOS arrays + // GOOD - Use CHAI managed arrays for kernels array1d values; arrayView1d valuesView = values.toView(); arrayView1d constView = values.toViewConst(); + stackArray1d fixedData; Tensor Types ^^^^^^^^^^^^ @@ -142,6 +144,8 @@ Use GEOS tensor types for geometric and mechanical calculations: It is possible to use ``std`` library components for data on host memory, for doing so, the rule is to **only use GEOS ``std`` container wrappers** instead of direct standard library containers. +This rule allow us to control bounds checking depending on ``GEOS_USE_BOUNDS_CHECK`` macro / ``GEOS_ENABLE_BOUNDS_CHECK`` cmake option.. + * - Standard C++ Type - GEOS Type - Description @@ -169,8 +173,8 @@ to **only use GEOS ``std`` container wrappers** instead of direct standard libra std::map orderedMap; // GOOD - GEOS wrappers - array1d data; - stackArray1d fixedData; + stdVector data; + stdArray fixedData; map orderedMap; The following standard library components are allowed: From f63e53d1ae3e03d1e0d26ab41717da64c8f592d9 Mon Sep 17 00:00:00 2001 From: Herve Gross <40979822+herve-gross@users.noreply.github.com> Date: Wed, 10 Dec 2025 18:01:37 -0800 Subject: [PATCH 05/13] Apply suggestions from code review --- src/docs/sphinx/developerGuide/Contributing/CodeRules.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst b/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst index 9e1a69a59ab..1bcd914d716 100644 --- a/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst +++ b/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst @@ -234,7 +234,7 @@ Keep Doxygen comments and naming as clear as possible. Keep in mind that it is t :icon: code .. code-block:: c++ - // GOOD - Documentation that gives useful information about what the function *does*, without + // GOOD - The documentation gives useful information about what the function *does*, without // going into theorical knowledge nor implementation details. /** * @brief Computes the residual for the nonlinear solver. From 205c603cc8c75e3d362cc8676c355f6df4f2d44f Mon Sep 17 00:00:00 2001 From: MelReyCG Date: Wed, 17 Dec 2025 20:42:26 +0100 Subject: [PATCH 06/13] =?UTF-8?q?=F0=9F=93=9D=20first=20pass=20of=20review?= =?UTF-8?q?s?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../developerGuide/Contributing/CodeRules.rst | 64 +++++++++++-------- 1 file changed, 39 insertions(+), 25 deletions(-) diff --git a/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst b/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst index 1bcd914d716..d077b8492a0 100644 --- a/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst +++ b/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst @@ -178,7 +178,7 @@ This rule allow us to control bounds checking depending on ``GEOS_USE_BOUNDS_CHE map orderedMap; The following standard library components are allowed: -- ``std::pair``, ``std::tuple`` for temporary structures, but beware of memory allocation, sometimes struct are to prefer. +- ``std::pair``, ``std::tuple`` for temporary structures, `sometimes struct are to prefer `_. - ``std::function`` for callbacks - ``std::optional`` for optional return values - ``std::move``, ``std::forward`` for move semantics @@ -228,7 +228,7 @@ User-facing features must be documented in the Sphinx RST documentation (``src/d Doxygen & Naming (Developer-Oriented) """""""""""""""""""""""""""""""""""""" -Keep Doxygen comments and naming as clear as possible. Keep in mind that it is target to developers from other domains. +Keep Doxygen comments and naming as clear as possible. Variable names and inline comments should provide help for developers even if they are not experts in the domain. .. dropdown:: Example: Doxygen documentation :icon: code @@ -253,7 +253,7 @@ Keep Doxygen comments and naming as clear as possible. Keep in mind that it is t real64 const dt, arrayView1d const & residual ); - // BAD - Documentation does not tell anything about the point of this function. + // BAD - documentation does not provide useful information beyond the function name. /** * @brief Compute CFL numbers * @param[in] time current time @@ -271,8 +271,9 @@ Testing Unit Tests ^^^^^^^^^^ -**Always provide unit tests** to ensure code behaves according to expectations. -Unit testing is not about testing every single line of code (an unrealistic goal), but about *verifying assumptions and preventing regressions*. Don't assume your code works as intended—prove it with tests. +**Always provide unit tests** to ensure every function behaves according to expectations. +Unit testing is not about testing every single line of code (an unrealistic goal), but about *verifying assumptions and preventing regressions*. Do not assume that your code will always work as intended — protect it with unit tests. + - Focus on **critical logic** (core functionality and algorithms), **edge cases** (extreme, erroneous, empty values). - **Use GoogleTest framework**. - In GEOS, the goal of unit test is never to take position on the validity of physical results (this is the point of integrated-tests). @@ -307,6 +308,19 @@ Unit testing is not about testing every single line of code (an unrealistic goal EXPECT_THROW( component.compute( -1.0 ), std::invalid_argument ); } + // Test the `testMyFunc()` function for GPU (geos::parallelDevicePolicy, if applicable) and + // CPU scenario (serialPolicy) + #ifdef GEOS_USE_DEVICE + TEST( MyFuncTests, testMyFuncDevice ) + { + testMyFunc< geos::parallelDevicePolicy< > >(); + } + #endif + TEST( MyFuncTests, testMyFuncHost ) + { + testMyFunc< serialPolicy >(); + } + // entry point of the test binary int main( int argc, char * * argv ) { @@ -321,14 +335,14 @@ Code Coverage ^^^^^^^^^^^^^ **Code coverage should never decrease.** New code contributions must maintain or improve overall code coverage. -Use code-cov reports to identify untested code paths. +Use code-cov to report untested code paths. -Integrated Tests & examples +Integrated Tests & Examples ^^^^^^^^^^^^^^^^^^^^^^^^^^^ -**Every XML-usable & serializable component must have examples (at least, have an integrated-test case):** -- Each different usage pattern requires its own example, -- Examples must be ran as integrated tests to guaranty they keep valid, +**Always provide an example or an integrated test for every XML-usable & serializable component:** +- Each usage pattern should have its own dedicated example, +- Examples must run as part of the integrated tests to protect features against regressions, - Place examples in ``examples/`` directory, and integrated tests in ``inputFiles/``, For more info, please refer to :ref:`IntegratedTests.rst`. @@ -339,7 +353,7 @@ Logging Rules ----- -Use GEOS Logging Macros (``GEOS_LOG*``) +GEOS Logging Macros (``GEOS_LOG*``) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ **Use of ``std::cout``, ``std::cerr``, or ``printf`` for logging must never appear** on the develop branch. @@ -368,7 +382,7 @@ Always use GEOS logging macros defined in ``Logger.hpp``: - ``GEOS_LOG_RANK(msg)`` - Log to rank-specific output stream - ``GEOS_LOG_RANK_VAR(var)`` - Log variable on rank stream -Use Log Levels +Using Log Levels ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ **Logging should be meaningful and controlled**: @@ -412,7 +426,7 @@ See :doc:`LogLevels` for using / adding log level (e.g., ``logInfo::Convergence` Logging in RAJA Kernels ^^^^^^^^^^^^^^^^^^^^^^^^ -**Logging inside RAJA kernels only for DEBUG purposes** and must never appear on the develop branch. +**Logging inside RAJA kernels should only be done for debugging purposes** and must never appear on the develop branch. **Rationale:** Logging on GPU can: - Reserve cache lines and degrade performance, @@ -497,7 +511,7 @@ Error Management ================ - **Developpers must provide runtime validations** of their code behaviour, -- for runtime validation fails that can be encountered by users, **meaningful messages and context must be provided**. +- For runtime validation fails that can be triggered by users, **meaningful messages and context must be provided** to help troubleshooting. Errors (GEOS_ERROR*) -------------------- @@ -586,7 +600,7 @@ do not serve as an alternative code path to consult a behaviour state. - **Performance:** Exception handling is expensive, especially with stack unwinding, - **Predictability:** Catching and continuing makes behaviour unpredictable and hard to debug, -- **Safety:** The error state may have corrupted simulation data, broke the code path at an unexpected place, and/or left the system in an invalid state. +- **Safety:** The error state may have corrupted simulation data, broken the code path at an unexpected place, and/or left the system in an invalid state. .. dropdown:: Example: Exception handling anti-pattern :icon: code @@ -622,7 +636,7 @@ Use ``GEOS_WARNING*`` macros when: - An issue is detected but **simulation can continue** - Simulation may be **affected but not completly invalid** -- User should be notified of a **potential problem** +- The user should be notified of a **potential problem** .. dropdown:: Example: Warning usage :icon: code @@ -638,10 +652,10 @@ Use ``GEOS_WARNING*`` macros when: Errors & Warning in RAJA Kernels -------------------------------- -Same rule as logging: **Error and warning output inside RAJA kernels only for DEBUG purposes** and must never appear on the develop branch. +Similarly to Logging, **error and warning output inside RAJA kernels is only for debugging purposes** and must never appear on the develop branch. - GPU kernel errors cause immediate kernel termination, adding potentially costly branches, -- Error handling on device has significant performance & cache impact, +- Error handling on device has significant performance and cache impact, - GPU error handling may be unsupported depending on the platform / device, - Can cause deadlocks in parallel execution. @@ -881,7 +895,7 @@ Communication Batching Validation / Precision ====================== -Computation validation +Computation Validation ----------------------- Use Tolerance-Based Comparisons @@ -907,10 +921,10 @@ Use Tolerance-Based Comparisons real64 const scale = std::max( std::abs(a), std::abs(b) ); if( std::abs(a - b) < relTol * scale ) { ... } -Division Safety, NaN/Inf values +Division Safety, NaN/Inf Values ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -- **Always verify denominators are not zero or near-zero before division**. +- **Always verify that denominators are not zero or near-zero before a division**. - In General, **we should not make any computation that could result in NaN/Inf values**. @@ -1024,7 +1038,7 @@ Performance Profiling --------- -Profile When Affecting Kernel behaviour +Profile When Affecting Kernel Behaviour ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ When modifying performance-critical code (kernels, assembly loops, solvers): @@ -1081,7 +1095,7 @@ Speed Optimization Rules ----------------------------- - **Hoist Loop Invariants**, move computations that don't change during iterations outside the loop. -- When it does not critically affect the code architecture and clarity, **Fuse multiple related kernels to reduce memory traffic and launch overhead** (i.e., statistics kernels process all physics field at once). +- When it does not critically affect the code architecture and clarity, **fuse multiple related kernels to reduce memory traffic and launch overhead** (i.e., statistics kernels process all physics field at once). - **Optimize Memory Access for Cache and Coalescing**. Access memory sequentially and ensure coalesced access, especially on GPUs. - **Minimize Host-Device Transfers**. Keep data on the appropriate memory space and minimize transfers. @@ -1104,7 +1118,7 @@ Be Const / Constexpr - **Code safety,** prevents accidental modification for constant contexts, - **Show code intention,** make code more readable. -Also, **Mark methods const if the method is not designed to modify** the object state. +Also, **mark methods const if the method is not designed to modify** the object state. Value / Const Function Arguments ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -1297,7 +1311,7 @@ Principles ^^^^^^^^^^ - **Loose coupling:** Components should depend on interfaces, not concrete implementations, -- **No circular dependancies:** Consider the existing GEOS dependancies to not make components co-dependant (headers inclusion, packages referencing in ``CMakeLists.txt``, avoid tightly coupled objects), +- **No circular dependencies:** Consider the existing GEOS dependencies to not make components co-dependent (headers inclusion, packages referencing in ``CMakeLists.txt``, avoid tightly coupled objects), - **Dependency injection:** Public components should receive their dependencies from external sources. Pass required dependencies using intermediate types instead of direct implementation types, using lambda, templates, and minimal interfaces, relying on **lambda**, **templates** and **minimal interfaces** (loose coupling, testability), - **Performance exceptions:** Tight coupling is acceptable when required for performance, - **Minimize header inclusions and dependancies**. From e379ac3892b3ec7432b6aac87bf951e3d132a578 Mon Sep 17 00:00:00 2001 From: MelReyCG Date: Thu, 18 Dec 2025 15:28:37 +0100 Subject: [PATCH 07/13] =?UTF-8?q?=F0=9F=93=9D=20Continuing=20on=20reviews,?= =?UTF-8?q?=20removed=20"Rationale:"=20in=20favor=20of=20dropdowns?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../developerGuide/Contributing/CodeRules.rst | 94 ++++++++++++------- 1 file changed, 58 insertions(+), 36 deletions(-) diff --git a/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst b/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst index d077b8492a0..af91419e75c 100644 --- a/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst +++ b/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst @@ -132,11 +132,11 @@ Use GEOS tensor types for geometric and mechanical calculations: * - Type - Description * - ``R1Tensor`` - - 3D vector (real64) + - 3D vector (real64, default choice over 32bit counterpart) * - ``R1Tensor32`` - 3D vector (real32) * - ``R2SymTensor`` - - Symmetric 6-component tensor (Voigt notation) + - Symmetric 6-component tensor (`Voigt notation `_) ``std`` Standard Library ------------------------ @@ -178,7 +178,7 @@ This rule allow us to control bounds checking depending on ``GEOS_USE_BOUNDS_CHE map orderedMap; The following standard library components are allowed: -- ``std::pair``, ``std::tuple`` for temporary structures, `sometimes struct are to prefer `_. +- ``std::pair``, ``std::tuple`` for temporary structures, `but sometimes struct are to prefer `_. - ``std::function`` for callbacks - ``std::optional`` for optional return values - ``std::move``, ``std::forward`` for move semantics @@ -340,7 +340,8 @@ Use code-cov to report untested code paths. Integrated Tests & Examples ^^^^^^^^^^^^^^^^^^^^^^^^^^^ -**Always provide an example or an integrated test for every XML-usable & serializable component:** +**Always provide an example or an integrated test for every XML-usable & serializable component:** + - Each usage pattern should have its own dedicated example, - Examples must run as part of the integrated tests to protect features against regressions, - Place examples in ``examples/`` directory, and integrated tests in ``inputFiles/``, @@ -386,6 +387,7 @@ Using Log Levels ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ **Logging should be meaningful and controlled**: + - When possible, **use ``GEOS_LOG_LEVEL_RANK*`` macros with appropriate ``logInfo``** to allow output filtering, - **Avoid unnecessary logging**, prefer rank 0 logging for global information to avoid redundant output, - **Consider logs performance impact** (output flush frequency, formatting). @@ -428,11 +430,15 @@ Logging in RAJA Kernels **Logging inside RAJA kernels should only be done for debugging purposes** and must never appear on the develop branch. -**Rationale:** Logging on GPU can: -- Reserve cache lines and degrade performance, -- Cause race conditions in output, -- Generate massive amounts of output, -- Produce unpredictable behaviour (e.g. runtime crashes on AMD). +.. dropdown:: Why logging on GPU is not allowed on production code? + :icon: info + + Logging on GPU can: + + - Reserve cache lines and degrade performance, + - Cause race conditions in output, + - Generate massive amounts of output, + - Produce unpredictable behaviour (e.g. runtime crashes on AMD). .. dropdown:: Example: Debug-only kernel logging :icon: code @@ -596,7 +602,7 @@ Never Catch and Continue **Do not catch exceptions and continue execution**. Exceptions in GEOS indicate serious problems, they do not serve as an alternative code path to consult a behaviour state. -**Rationale:** +Why not to use exceptions for anything else than passing information up to the call stack: - **Performance:** Exception handling is expensive, especially with stack unwinding, - **Predictability:** Catching and continuing makes behaviour unpredictable and hard to debug, @@ -654,10 +660,13 @@ Errors & Warning in RAJA Kernels Similarly to Logging, **error and warning output inside RAJA kernels is only for debugging purposes** and must never appear on the develop branch. -- GPU kernel errors cause immediate kernel termination, adding potentially costly branches, -- Error handling on device has significant performance and cache impact, -- GPU error handling may be unsupported depending on the platform / device, -- Can cause deadlocks in parallel execution. +.. dropdown:: Why errors are forbidden on GPU in production code? + :icon: info + + - GPU kernel errors cause immediate kernel termination, adding potentially costly branches, + - Error handling on device has significant performance and cache impact, + - GPU error handling may be unsupported depending on the platform / device, + - Can cause deadlocks in parallel execution. .. dropdown:: Example: Kernel error handling (debug only) :icon: code @@ -843,6 +852,7 @@ Always Use MpiWrapper ^^^^^^^^^^^^^^^^^^^^^ **Always use ``MpiWrapper`` instead of raw MPI calls.** + Refer to ``MpiWrapper.hpp`` for more. .. dropdown:: Example: Using MpiWrapper @@ -867,10 +877,15 @@ Refer to ``MpiWrapper.hpp`` for more. Communication Batching ^^^^^^^^^^^^^^^^^^^^^^^ -**Rule: Batch Small Communications** - Combine multiple small messages into larger ones to reduce communication overhead. +**Batch small communications as much as possible**. + +.. dropdown:: Why? + :icon: info - **Rationale:** MPI has significant per-message overhead; batching improves performance. + As MPI has significant per-message overhead, combining multiple small messages into larger ones can reduce communication overhead. + +.. dropdown:: Example: Batching MPI communication + :icon: code .. code-block:: cpp @@ -914,7 +929,7 @@ Use Tolerance-Based Comparisons // GOOD - Absolute tolerance real64 const absTol = 1.0e-12; - if( std::abs(value) < absTol ) { ... } + if( LvArray::math::abs(value) < absTol ) { ... } // GOOD - Relative tolerance real64 const relTol = 1.0e-8; @@ -925,7 +940,6 @@ Division Safety, NaN/Inf Values ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - **Always verify that denominators are not zero or near-zero before a division**. - - In General, **we should not make any computation that could result in NaN/Inf values**. .. dropdown:: Example: Division Safety @@ -961,9 +975,7 @@ Division Safety, NaN/Inf Values Overflow Prevention ^^^^^^^^^^^^^^^^^^^^ -Overflow leads to undefined behavior and memory corruption. - -**Validate that operations won't exceed type limits**, especially for index calculations. +Overflow leads to undefined behavior and memory corruption. **Validate that operations won't exceed type / container limits**, especially for index calculations. Use LvArray Math Utilities ^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -1124,6 +1136,7 @@ Value / Const Function Arguments ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ **Pass large and non-trivially copyable objects by const reference**: + - A "large" size can be defined as "more that 16 bytes, ", which is 2 pointers, 2 integer / index, or 2 double, - An object is non-trivially copyable objects when it needs to perform sub-allocation when being copied, @@ -1198,9 +1211,12 @@ This rule imply the following: - **A reference may never be null**, - In a given method, **``this`` may never be null**. -Passing a pointer or reference parameter show the intention of modifying the underlying object. -The difference of the two practices is that passing a pointer means that we consider that the -underlying object can be null. +.. dropdown:: Why? + :icon: info + + Passing a pointer or reference parameter show the intention of modifying the underlying object. + The difference of the two practices is that passing a pointer means that we consider that the + underlying object can be null. Constexpr for Compile-Time Constants ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -1231,11 +1247,12 @@ Provide Views to Arrays The rule is generalizable to ``string_view`` for strings, but not applicable in GPU context. -Why Use Views? +.. dropdown:: Why Use Views? + :icon: info -- **No memory allocation:** Views are lightweight references -- **Mutability correctness:** Can provide ``const`` read-only access to inner data -- **GPU compatibility:** Views work seamlessly on device + - **No memory allocation:** Views are lightweight references + - **Mutability correctness:** Can provide ``const`` read-only access to inner data + - **GPU compatibility:** Views work seamlessly on device .. dropdown:: Example: Views for arrays :icon: code @@ -1293,10 +1310,11 @@ The rule is applicable to ``arrayView*`` and ``string_view``. // CORRECT - store array, create view when needed class SafeDataHolder { - array1d< real64 > m_data; - + public: arrayView1d< real64 > getView() { return m_data.toView(); } arrayView1d< real64 const > getViewConst() const { return m_data.toViewConst(); } + private: + array1d< real64 > m_data; }; General Architecture @@ -1355,12 +1373,13 @@ Avoid Globals and Mutable State **Minimize mutable global states when possible.** Prefer passing context explicitly. -Why to avoid globals: +.. dropdown:: Why to avoid globals? + :icon: info -- **Thread safety:** Globals can cause race conditions in parallel code -- **Testability:** Hard to test code that depends on global state -- **Predictability:** Global state makes code behaviour harder to understand -- **Encapsulation:** Violates modularity principles + - **Thread safety:** Globals can cause race conditions in parallel code + - **Testability:** Hard to test code that depends on global state + - **Predictability:** Global state makes code behaviour harder to understand + - **Encapsulation:** Violates modularity principles .. dropdown:: Example: Avoiding global state :icon: code @@ -1404,6 +1423,7 @@ Magic values, ``Group`` paths ------------------------------- **Avoid magic values**: + - **arbitrary values should not be written more than once**, define constants, consider using or extending ``PhysicsConstants.hpp`` / ``Units.hpp``, - Because of architecture considerations, **Avoid using full raw ``Group`` path**, especially when the path has parts unrelated to the component consulting it -> prefer as relative paths as possible, - ``Group`` / ``Wrapper`` name in the data repository is considered as a magic value, **use the ``getCatalogName()`` / ``getName()`` getters**, @@ -1427,6 +1447,8 @@ Unphysical values indicate errors and can cause solver failures. If a value is not strictly disallowed but does not seem possible for the model, **show a warning to the user that he can disable**. +.. dropdown:: Example: Avoiding global state + :icon: code .. code-block:: cpp From 61f807ce4375bea5432c0052632754cdaacbb52f Mon Sep 17 00:00:00 2001 From: MelReyCG Date: Tue, 30 Dec 2025 16:36:08 +0100 Subject: [PATCH 08/13] =?UTF-8?q?=F0=9F=93=9D=20Added=20GEOS=5FCALIPER=5FM?= =?UTF-8?q?ARK=5FSCOPE=20for=20profiling=20docs?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/docs/sphinx/developerGuide/Contributing/CodeRules.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst b/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst index af91419e75c..5cdcefd5c7e 100644 --- a/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst +++ b/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst @@ -1079,7 +1079,7 @@ When modifying performance-critical code (kernels, assembly loops, solvers): Caliper Integration ^^^^^^^^^^^^^^^^^^^^ -Use ``GEOS_MARK_FUNCTION`` and ``Timer`` for performance tracking for the main computation functions / scopes. +Use ``GEOS_MARK_FUNCTION``, ``GEOS_CALIPER_MARK_SCOPE`` and ``Timer`` for performance tracking for the main computation functions / scopes. .. dropdown:: Example: Performance instrumentation :icon: code From 0f2dbf8fde04150dd2896080da29955caca77e19 Mon Sep 17 00:00:00 2001 From: MelReyCG Date: Tue, 30 Dec 2025 16:47:47 +0100 Subject: [PATCH 09/13] =?UTF-8?q?=F0=9F=93=9Drewrite=20of=20data=20context?= =?UTF-8?q?=20docs?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../sphinx/developerGuide/Contributing/CodeRules.rst | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst b/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst index 5cdcefd5c7e..82457d3c63d 100644 --- a/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst +++ b/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst @@ -691,13 +691,14 @@ information on the source of the data, which can be: - the file & line of the simulation XML file, - if no source exists, ``Group`` / ``Wrapper`` path in the data repository. -**Add ``DataContext`` when** Errors occur related to a **data repository component**: +**Add ``DataContext`` when** an error occurs because of: -- a ``Wrapper`` when validation fails for a given user input setting, -- a ``Group`` when the error is implied by the instance state, +- a ``Wrapper`` when validation fails for a given user input setting (``getWrapperDataContext()``), +- a ``Group`` when the error is implied by the instance state (``getDataContext()``), - multiple data-contexts when the error is due to multiple objects / settings (first = more important). -Use ``getDataContext()`` or ``getWrapperDataContext()`` as last parameters of the logging macros to give context. +**Use ``getDataContext()`` or ``getWrapperDataContext()`` as last parameters of the logging macros to pass the data context metadata.** +It will be provided properly in the error message depending on the format (log / YAML). .. dropdown:: Example: Using DataContext :icon: code @@ -718,6 +719,7 @@ Use ``getDataContext()`` or ``getWrapperDataContext()`` as last parameters of th getDataContext() ); // Exemple 3, using multiple DataContext because the error is related with multiple objects + // The first ones are considered the highest priority. GEOS_THROW_IF( m_isThermal && !isFluidModelThermal, GEOS_FMT( "CompositionalMultiphaseBase {}: the thermal option is enabled in the solver, but the fluid model {} is incompatible with the thermal option", getDataContext(), fluid.getDataContext() ), From 62ba67a2f286af722729974ee29b38a12362ce74 Mon Sep 17 00:00:00 2001 From: MelReyCG Date: Tue, 30 Dec 2025 17:06:15 +0100 Subject: [PATCH 10/13] =?UTF-8?q?=F0=9F=93=9Dadded=20ref=20for=20cuda=20tr?= =?UTF-8?q?ap=20issue?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/docs/sphinx/developerGuide/Contributing/CodeRules.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst b/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst index 82457d3c63d..da4ed819e5b 100644 --- a/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst +++ b/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst @@ -666,7 +666,7 @@ Similarly to Logging, **error and warning output inside RAJA kernels is only for - GPU kernel errors cause immediate kernel termination, adding potentially costly branches, - Error handling on device has significant performance and cache impact, - GPU error handling may be unsupported depending on the platform / device, - - Can cause deadlocks in parallel execution. + - Can cause deadlocks in parallel execution (`as discussed here`_, a rank can die, then the other ranks will wait for it at next MPI call). .. dropdown:: Example: Kernel error handling (debug only) :icon: code From 2d6af957e2e93f2f57ceb5e73c1de947637efa8b Mon Sep 17 00:00:00 2001 From: MelReyCG Date: Tue, 30 Dec 2025 17:37:07 +0100 Subject: [PATCH 11/13] =?UTF-8?q?=F0=9F=93=9D=20small=20precision=20added?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/docs/sphinx/developerGuide/Contributing/CodeRules.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst b/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst index da4ed819e5b..f54a457232f 100644 --- a/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst +++ b/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst @@ -1254,7 +1254,7 @@ The rule is generalizable to ``string_view`` for strings, but not applicable in - **No memory allocation:** Views are lightweight references - **Mutability correctness:** Can provide ``const`` read-only access to inner data - - **GPU compatibility:** Views work seamlessly on device + - **GPU compatibility:** LvArray views work seamlessly on device .. dropdown:: Example: Views for arrays :icon: code From 3eaa3520e225f50c9ed79ab91c56dccc82ef8dbf Mon Sep 17 00:00:00 2001 From: MelReyCG Date: Wed, 31 Dec 2025 10:40:36 +0100 Subject: [PATCH 12/13] =?UTF-8?q?=F0=9F=93=9Dchange=20docs=20title=20level?= =?UTF-8?q?s?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../sphinx/developerGuide/Contributing/CodeRules.rst | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst b/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst index f54a457232f..7b4f5deafd2 100644 --- a/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst +++ b/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst @@ -191,13 +191,10 @@ Contribution Documentation ------------- -Always Provide Comprehensive Documentation -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -All new features must include appropriate documentation at multiple levels: +**All new features must include appropriate documentation, for developpers and user, at the following levels:** Wrapper Documentation (User-Oriented) -"""""""""""""""""""""""""""""""""""""" +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ All data repository wrappers must be documented with ``setDescription()``. As much as possible, valid values rules should be provided, and default values should be provided with ``setApplyDefaultValue()``. @@ -215,7 +212,7 @@ As much as possible, valid values rules should be provided, and default values s "Valid range: (0, 1e-3]. See User Guide Chapter 5." ); RST Documentation (User Guide) -""""""""""""""""""""""""""""""" +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ User-facing features must be documented in the Sphinx RST documentation (``src/docs/sphinx/``): @@ -226,7 +223,7 @@ User-facing features must be documented in the Sphinx RST documentation (``src/d - Add **references** to related features Doxygen & Naming (Developer-Oriented) -"""""""""""""""""""""""""""""""""""""" +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Keep Doxygen comments and naming as clear as possible. Variable names and inline comments should provide help for developers even if they are not experts in the domain. From db6f8507fc41fb763bf3de89936a3248d4e5c136 Mon Sep 17 00:00:00 2001 From: MelReyCG Date: Wed, 31 Dec 2025 16:20:39 +0100 Subject: [PATCH 13/13] =?UTF-8?q?=F0=9F=93=9D=20splitting=20code=20rules?= =?UTF-8?q?=20document=20in=202=20+=20adding=20indexation?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../Contributing/CodeGeosFeatures.rst | 1100 ++++++++++++ .../developerGuide/Contributing/CodeRules.rst | 1470 +---------------- .../Contributing/CodingPractices.rst | 415 +++++ .../Contributing/index_contributing.rst | 2 + 4 files changed, 1526 insertions(+), 1461 deletions(-) create mode 100644 src/docs/sphinx/developerGuide/Contributing/CodeGeosFeatures.rst create mode 100644 src/docs/sphinx/developerGuide/Contributing/CodingPractices.rst diff --git a/src/docs/sphinx/developerGuide/Contributing/CodeGeosFeatures.rst b/src/docs/sphinx/developerGuide/Contributing/CodeGeosFeatures.rst new file mode 100644 index 00000000000..3513e0dbb33 --- /dev/null +++ b/src/docs/sphinx/developerGuide/Contributing/CodeGeosFeatures.rst @@ -0,0 +1,1100 @@ +.. _CodeGeosFeatures: + + +############################################################################### +GEOS Coding Features Rules +############################################################################### + +This section covers the coding rules regarding the features available in the GEOS C++ infrastructure. + +All contributors should familiarize themselves with these rules to maintain code quality, consistency, and reliability across the GEOS codebase. + +.. note:: + This document uses collapsible code blocks. Click on "Show/Hide Code" to expand examples. + +Contribution Requirements +========================= + +This section contains the rules when adding features to GEOS. + +Documentation +------------- + +**All new features must include appropriate documentation, for developpers and user, at the following levels:** + +Wrapper Documentation (User-Oriented) +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +All data repository wrappers must be documented with ``setDescription()``. +As much as possible, valid values rules should be provided, and default values should be provided with ``setApplyDefaultValue()``. + +.. dropdown:: Example: Wrapper documentation + :icon: code + + .. code-block:: c++ + + registerWrapper( "timeStep", &m_timeStep ). + setInputFlag( InputFlags::REQUIRED ). + setApplyDefaultValue( 1.0e-6 ). + setDescription( "Initial time step value. This parameter controls " + "the starting time increment for the simulation. " + "Valid range: (0, 1e-3]. See User Guide Chapter 5." ); + +RST Documentation (User Guide) +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +User-facing features must be documented in the Sphinx RST documentation (``src/docs/sphinx/``): + +- Explain **what** the feature does +- Provide **usage examples** +- Document **input parameters** +- Include **expected output** or behaviour +- Add **references** to related features + +Doxygen & Naming (Developer-Oriented) +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Keep Doxygen comments and naming as clear as possible. Variable names and inline comments should provide help for developers even if they are not experts in the domain. + +.. dropdown:: Example: Doxygen documentation + :icon: code + + .. code-block:: c++ + // GOOD - The documentation gives useful information about what the function *does*, without + // going into theorical knowledge nor implementation details. + /** + * @brief Computes the residual for the nonlinear solver. + * @param[in] domain The domain partition containing mesh and field data + * @param[in] time Current simulation time + * @param[in] dt Time step size + * @param[out] residual The computed residual vector + * @return Residual L2 norm + * + * This method assembles the discrete residual for the physics equations + * by looping over all elements in the target regions. The residual is + * stored in the provided vector and its L2 norm is returned. + */ + real64 computeResidual( DomainPartition & domain, + real64 const time, + real64 const dt, + arrayView1d const & residual ); + + // BAD - documentation does not provide useful information beyond the function name. + /** + * @brief Compute CFL numbers + * @param[in] time current time + * @param[in] dt the time step size + * @param[in] domain the domain partition + */ + void computeCFLNumbers( real64 const time, + real64 const dt, + DomainPartition & domain ) const; + + +Testing +------- + +Unit Tests +^^^^^^^^^^ + +**Always provide unit tests** to ensure every function behaves according to expectations. +Unit testing is not about testing every single line of code (an unrealistic goal), but about *verifying assumptions and preventing regressions*. Do not assume that your code will always work as intended — protect it with unit tests. + +- Focus on **critical logic** (core functionality and algorithms), **edge cases** (extreme, erroneous, empty values). +- **Use GoogleTest framework**. +- In GEOS, the goal of unit test is never to take position on the validity of physical results (this is the point of integrated-tests). +- Place tests in ``src/coreComponents//tests/``. +- Test GPU code paths if applicable (use ``#ifdef GEOS_USE_DEVICE``). + +.. dropdown:: Example: Unit test structure + :icon: code + + .. code-block:: c++ + + #include + + // Test core functionality + TEST( MyComponentTest, ComputesCorrectResult ) + { + MyComponent component; + EXPECT_DOUBLE_EQ( component.compute( 2.0 ), 4.0 ); + } + + // Test edge case + TEST( MyComponentTest, HandlesZeroInput ) + { + MyComponent component; + EXPECT_DOUBLE_EQ( component.compute( 0.0 ), 0.0 ); + } + + // Test error detection + TEST( MyComponentTest, RejectsInvalidInput ) + { + MyComponent component; + EXPECT_THROW( component.compute( -1.0 ), std::invalid_argument ); + } + + // Test the `testMyFunc()` function for GPU (geos::parallelDevicePolicy, if applicable) and + // CPU scenario (serialPolicy) + #ifdef GEOS_USE_DEVICE + TEST( MyFuncTests, testMyFuncDevice ) + { + testMyFunc< geos::parallelDevicePolicy< > >(); + } + #endif + TEST( MyFuncTests, testMyFuncHost ) + { + testMyFunc< serialPolicy >(); + } + + // entry point of the test binary + int main( int argc, char * * argv ) + { + ::testing::InitGoogleTest( &argc, argv ); + g_commandLineOptions = *geos::basicSetup( argc, argv ); + int const result = RUN_ALL_TESTS(); + geos::basicCleanup(); + return result; + } + +Code Coverage +^^^^^^^^^^^^^ + +**Code coverage should never decrease.** New code contributions must maintain or improve overall code coverage. +Use code-cov to report untested code paths. + +Integrated Tests & Examples +^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +**Always provide an example or an integrated test for every XML-usable & serializable component:** + +- Each usage pattern should have its own dedicated example, +- Examples must run as part of the integrated tests to protect features against regressions, +- Place examples in ``examples/`` directory, and integrated tests in ``inputFiles/``, + +For more info, please refer to :ref:`IntegratedTests.rst`. + +Code Infrastructure +==================== + +**Always use GEOS types / functions instead of standard / external counterparts.** + +Typing +------------------------- + +Basic Types +^^^^^^^^^^^ + +.. list-table:: Standard C++ vs GEOS Types + :header-rows: 1 + :widths: 30 30 40 + + * - Standard C++ Type + - GEOS Type + - Description + * - ``int`` + - ``integer`` + - Signed integer + * - ``std::size_t`` + - ``localIndex`` + - Local array indexing (MPI partition) + * - N/A + - ``globalIndex`` + - Global indexing (across MPI) + * - ``float`` + - ``real32`` + - 32-bit floating point + * - ``double`` + - ``real64`` + - 64-bit floating point + * - ``std::string`` + - ``string`` + - String type + * - ``std::string_view`` + - ``string_view`` + - Non-owning string view + +.. dropdown:: Example: Using GEOS basic types + :icon: code + + .. code-block:: c++ + + // BAD - Do not use standard types + int count = 0; + double value = 1.5; + + // GOOD - Use GEOS types + integer count = 0; + real64 value = 1.5; + localIndex elementIndex = 42; + +Array Types +^^^^^^^^^^^ + +GEOS provides multi-dimensional CHAI array types for managed host-device data. + +Refer to the rule of :ref:`CHAI Memory Management`. + +.. list-table:: GEOS Array Types + :header-rows: 1 + :widths: 35 65 + + * - GEOS Type + - Description + * - ``array1d`` + - 1D owning array + * - ``arrayView1d`` + - 1D non-owning view (modifiable) + * - ``arraySlice1d`` + - 1D slice (from multi-dim array) + * - ``array2d`` + - 2D owning array + * - ``arrayView2d`` + - 2D non-owning view + * - ``arraySlice2d`` + - 2D slice + * - ``array3d``, ``array4d``, ``array5d`` + - 3D/4D/5D owning arrays + * - ``arrayView3d``, etc. + - 3D/4D/5D views + * - ``SortedArray`` + - Sorted array with search capabilities + * - ``ArrayOfArrays`` + - Jagged 2D array (variable row sizes) + * - ``ArrayOfSets`` + - Array of sets (unique elements per row) + * - ``stackArray1d`` + - Stack-allocated array (max size N) + +.. dropdown:: Example: Using GEOS arrays + :icon: code + + .. code-block:: c++ + + // BAD - Do not use std::vector for data that can be packed or kernel addressable + std::vector values; + std::array fixedData; + + // GOOD - Use CHAI managed arrays for kernels + array1d values; + arrayView1d valuesView = values.toView(); + arrayView1d constView = values.toViewConst(); + stackArray1d fixedData; + +Tensor Types +^^^^^^^^^^^^ + +Use GEOS tensor types for geometric and mechanical calculations: + +.. list-table:: GEOS Tensor Types + :header-rows: 1 + :widths: 30 70 + + * - Type + - Description + * - ``R1Tensor`` + - 3D vector (real64, default choice over 32bit counterpart) + * - ``R1Tensor32`` + - 3D vector (real32) + * - ``R2SymTensor`` + - Symmetric 6-component tensor (`Voigt notation `_) + +External Dependencies +---------------------- + +``std`` Standard Library Usage +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +It is possible to use ``std`` library components for data on host memory, for doing so, the rule is +to **only use GEOS ``std`` container wrappers** instead of direct standard library containers. + +This rule allow us to control bounds checking depending on ``GEOS_USE_BOUNDS_CHECK`` macro / ``GEOS_ENABLE_BOUNDS_CHECK`` cmake option.. + + * - Standard C++ Type + - GEOS Type + - Description + * - ``std::vector`` + - ``stdVector`` + - ``std`` dynamically allocated array + * - ``std::array`` + - ``stdArray`` + - ``std`` fixed constexpr sized array + * - ``std::map`` + - ``map`` + - ``std`` sorted dictionary + * - ``std::unordered_map`` + - ``unordered_map`` + - ``std`` unsorted dictionary (hash-map) + +.. dropdown:: Example: Container usage + :icon: code + + .. code-block:: c++ + + // BAD - Direct std containers + std::vector data; + std::array fixedData; + std::map orderedMap; + + // GOOD - GEOS wrappers + stdVector data; + stdArray fixedData; + map orderedMap; + +The following standard library components are allowed: +- ``std::pair``, ``std::tuple`` for temporary structures, `but sometimes struct are to prefer `_. +- ``std::function`` for callbacks +- ``std::optional`` for optional return values +- ``std::move``, ``std::forward`` for move semantics +- Algorithms from ```` where appropriate +- ``std::unique_ptr``, ``std::shared_ptr`` for memory management + +LvArray Math Utilities +^^^^^^^^^^^^^^^^^^^^^^^ + +The ``LvArray::math`` namespace provides safe numerical utilities, **use LvArray math utilities instead of platform specific math functions.** + +.. dropdown:: Example: LvArray math utilities + :icon: code + + .. code-block:: c++ + + #include "LvArray/src/math.hpp" + + real64 const absValue = LvArray::math::abs( value ); + real64 const maxValue = LvArray::math::max( a, b ); + real64 const minValue = LvArray::math::min( a, b ); + +Logging +======= + +Rules +----- + +GEOS Logging Macros (``GEOS_LOG*``) +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +**Use of ``std::cout``, ``std::cerr``, or ``printf`` for logging must never appear** on the develop branch. +Always use GEOS logging macros defined in ``Logger.hpp``: + +.. dropdown:: Example: Correct logging usage + :icon: code + + .. code-block:: c++ + + // BAD - Do not use + std::cout << "Starting simulation" << std::endl; + printf("Value = %f\n", value); + + // GOOD - Use GEOS macros + GEOS_LOG("Starting simulation"); + GEOS_LOG_RANK_0("Master rank initialization complete"); + +**Available logging macros:** + +- ``GEOS_LOG(...)`` - Log on all ranks +- ``GEOS_LOG_VAR(...)`` - Log variable name and value +- ``GEOS_LOG_IF(condition, msg)`` - Conditional logging +- ``GEOS_LOG_RANK_0(msg)`` - Log only on rank 0 (MPI master) +- ``GEOS_LOG_RANK_0_IF(condition, msg)`` - Conditional logging on rank 0 +- ``GEOS_LOG_RANK(msg)`` - Log to rank-specific output stream +- ``GEOS_LOG_RANK_VAR(var)`` - Log variable on rank stream + +Using Log Levels +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +**Logging should be meaningful and controlled**: + +- When possible, **use ``GEOS_LOG_LEVEL_RANK*`` macros with appropriate ``logInfo``** to allow output filtering, +- **Avoid unnecessary logging**, prefer rank 0 logging for global information to avoid redundant output, +- **Consider logs performance impact** (output flush frequency, formatting). + +See :doc:`LogLevels` for using / adding log level (e.g., ``logInfo::Convergence``, ``logInfo::TimeStep``). + +.. dropdown:: Example: Using log levels + :icon: code + + .. code-block:: c++ + + // physics solver cpp - effective logging + GEOS_LOG_LEVEL_RANK_0( logInfo::Convergence, + GEOS_FMT( "Newton iteration {}: residual = {:4.2e}", + iter, residual ) ); + + // src/coreComponents/physicsSolvers/LogLevelsInfo.hpp - convergence logging config + struct Convergence + { + static constexpr int getMinLogLevel() { return 1; } + static constexpr std::string_view getDescription() { return "Convergence information"; } + }; + +.. dropdown:: Example: Avoiding log spam + :icon: code + + .. code-block:: c++ + + // BAD - Will spam output for every element on every ranks + forAll< parallelDevicePolicy<> >( numElems, [=] GEOS_HOST_DEVICE ( localIndex const ei ) + { + GEOS_LOG("Processing element " << ei); // DO NOT DO THIS + }); + + // GOOD - Log summary information + GEOS_LOG_RANK_0(GEOS_FMT("Processed {} elements", numElems)); + +Logging in RAJA Kernels +^^^^^^^^^^^^^^^^^^^^^^^^ + +**Logging inside RAJA kernels should only be done for debugging purposes** and must never appear on the develop branch. + +.. dropdown:: Why logging on GPU is not allowed on production code? + :icon: info + + Logging on GPU can: + + - Reserve cache lines and degrade performance, + - Cause race conditions in output, + - Generate massive amounts of output, + - Produce unpredictable behaviour (e.g. runtime crashes on AMD). + +.. dropdown:: Example: Debug-only kernel logging + :icon: code + + .. code-block:: c++ + + forAll< parallelDevicePolicy<> >( numElems, [=] GEOS_HOST_DEVICE ( localIndex const ei ) + { + GEOS_LOG("processing..."); // Compiles, but remove before merge to develop! + GEOS_LOG(GEOS_FMT("processing elem {}/{}", ei, numElems)); // Will crash on CUDA builds (unsupported GEOS_FMT) + printf("DEBUG: ei=%d\n", ei); // Compiles, but remove before merge to develop! + }); + +Structured Data Output +---------------------- + +For tabular or statistical output, use structured logging facilities instead of plain text logging: + +- Use ``TableTextFormatter`` with a ``logLevel`` XML setting for formatted text table output, +- Use ``TableCSVFormatter`` with a ``writeCSV()`` XML setting for structured data output. + +.. dropdown:: Example: Log table output + :icon: code + + .. code-block:: c++ + // log output preparation, layout precomputation can be kept for periodic outputs + TableLayout statsLogLayout( "", { "Flux(es)", "Region", "Element Count", massColumn, rateColumn } ); + m_logLayout = statsLogLayout; + + // filling table data (log format is for current timestep) + m_logData.addRow( fluxName, + elementSetName, + GEOS_FMT( "{}", wrappedStats.stats().m_elementCount ), + GEOS_FMT( "{}", wrappedStats.stats().m_producedMass ), + GEOS_FMT( "{}", wrappedStats.stats().m_productionRate ) ); + + // statistics CSV output (only if enabled and from rank 0) + if( logLevelActive && MpiWrapper::commRank() == 0 ) + { + string const title = GEOS_FMT( "{}, flux statistics for: {}", + getName(), fluxesStr ); + string const multilineTitle = wrapTextToMaxLength( title, 80 ); // enabling automatic line wrapping + m_logLayout.setTitle( multilineTitle ); + TableTextFormatter const tableStatFormatter( m_logLayout ); + GEOS_LOG_RANK( tableStatFormatter.toString( statsData ) ); + } + +.. dropdown:: Example: CSV output + :icon: code + + .. code-block:: c++ + + // CSV output preparation, layout precomputation can be kept for periodic outputs + TableLayout const statsCSVLayout( "", {"Time [s]", "Flux(es)", "Region", "Element Count", massColumn, rateColumn} ); + m_csvLayout = statsCSVLayout; + + // filling table data (CSV format is all timestep series) + m_csvData.addRow( wrappedStats.getStatsPeriodStart(), + fluxName, + elementSetName, + wrappedStats.stats().m_elementCount, + wrappedStats.stats().m_producedMass, + wrappedStats.stats().m_productionRate ); + + // statistics CSV output (only if enabled and from rank 0) + if( m_writeCSV > 0 && MpiWrapper::commRank() == 0 ) + { + std::ofstream outputFile( m_csvFilename ); + TableCSVFormatter const tableStatFormatter( m_csvLayout ); + outputFile << tableStatFormatter.toString( csvData ); + outputFile.close(); + csvData.clear(); + } + +Error Management +================ + +- **Developpers must provide runtime validations** of their code behaviour, +- For runtime validation fails that can be triggered by users, **meaningful messages and context must be provided** to help troubleshooting. + +Errors (GEOS_ERROR*) +-------------------- + +Why Use Errors? +^^^^^^^^^^^^^^^ + +Use ``GEOS_ERROR*`` macros when encountering a **blocking error** where: + +- The simulation result would be certifiably invalid, +- The application state is unrecoverable, +- Continuing execution would be unsafe or meaningless. + +.. dropdown:: Example: Error handling + :icon: code + + .. code-block:: c++ + + // Example: not implemented feature + GEOS_ERROR( "Rounding interpolation with derivatives not implemented" ); + + // Example: invalid mesh topology + GEOS_ERROR_IF( numFaces == 0, + GEOS_FMT("Element {} has zero faces - invalid mesh", elemId) ); + + // Example: comparison-based error + GEOS_ERROR_IF_LT( dt, 0.0, "Time step must be positive" ); + +Exceptions (GEOS_THROW*) +------------------------ + +Why Use Exceptions? +^^^^^^^^^^^^^^^^^^^ + +Use ``GEOS_THROW*`` macros for the same reasons as ``GEOS_ERROR*`` (**unrecoverable state**), and when: + +- You want to **add context** at higher call stack levels, +- You need to **propagate detailed error information** upward +- If you use custom exception class, you need to document them here. + +.. list-table:: Available exception types + :header-rows: 1 + :widths: 35 65 + + * - Error class + - Used for + * - ``std::runtime_error`` + - General runtime errors + * - ``geos::InputError`` + - Errors in user input (XML, data files) + * - ``geos::SimulationError`` + - Errors during simulation execution + * - ``geos::BadTypeError`` + - Type conversion errors + * - ``geos::NotAnError`` + - Non critical, volontary stopping state + +.. dropdown:: Example: Throwing exceptions + :icon: code + + .. code-block:: c++ + + // Example 1: throw with context + GEOS_THROW_IF( !file.is_open(), + GEOS_FMT("Cannot open file {}", filename), + InputError ); + + // Example 2: rethrow with context info (because of unclear exception) + try + { + tolerance = stod( inputParams[8] ); + } + catch( std::invalid_argument const & e ) + { + GEOS_THROW( GEOS_FMT( "{}: invalid model parameter value: {}", functionName, e.what() ), + InputError, getDataContext() ); + } + +Never Catch and Continue +^^^^^^^^^^^^^^^^^^^^^^^^^ + +**Do not catch exceptions and continue execution**. Exceptions in GEOS indicate serious problems, they +do not serve as an alternative code path to consult a behaviour state. + +Why not to use exceptions for anything else than passing information up to the call stack: + +- **Performance:** Exception handling is expensive, especially with stack unwinding, +- **Predictability:** Catching and continuing makes behaviour unpredictable and hard to debug, +- **Safety:** The error state may have corrupted simulation data, broken the code path at an unexpected place, and/or left the system in an invalid state. + +.. dropdown:: Example: Exception handling anti-pattern + :icon: code + + .. code-block:: c++ + + // BAD - Do not catch to continue + // Here, the "solver" should have a "isFailed" flag or a returned resulting state + try + { + solver.solve(); + } + catch( std::exception const & e ) // <- Note: this exception could have been throw for unexpected reason + { + GEOS_LOG("Solver failed, trying alternative"); + alternativeSolver.solve(); + } + + // GOOD - Manage your code paths manually + if(! solver.solve() ) + { + GEOS_LOG("Solver failed, trying alternative"); + alternativeSolver.solve(); + } + +Warnings (GEOS_WARNING*) +------------------------ + +Why Use Warnings? +^^^^^^^^^^^^^^^^^ + +Use ``GEOS_WARNING*`` macros when: + +- An issue is detected but **simulation can continue** +- Simulation may be **affected but not completly invalid** +- The user should be notified of a **potential problem** + +.. dropdown:: Example: Warning usage + :icon: code + + .. code-block:: c++ + + // Example: suboptimal but valid configuration + GEOS_WARNING_IF( dt > recommendedDt, + GEOS_FMT("Time step {} exceeds recommended value {}. This may affect solution accuracy.", + dt, recommendedDt), + getDataContext() ); + +Errors & Warning in RAJA Kernels +-------------------------------- + +Similarly to Logging, **error and warning output inside RAJA kernels is only for debugging purposes** and must never appear on the develop branch. + +.. dropdown:: Why errors are forbidden on GPU in production code? + :icon: info + + - GPU kernel errors cause immediate kernel termination, adding potentially costly branches, + - Error handling on device has significant performance and cache impact, + - GPU error handling may be unsupported depending on the platform / device, + - Can cause deadlocks in parallel execution (`as discussed here`_, a rank can die, then the other ranks will wait for it at next MPI call). + +.. dropdown:: Example: Kernel error handling (debug only) + :icon: code + + .. code-block:: c++ + + // On GPU, errors will cause kernel abort + forAll< parallelDevicePolicy<> >( numElems, [=] GEOS_HOST_DEVICE ( localIndex const ei ) + { + #if !defined( NDEBUG ) + // This will trap on GPU - use sparingly for debugging only + GEOS_ERROR_IF( badCondition, "Bad condition detected" ); + #endif + }); + +Contextualization with DataContext +----------------------------------- + +GEOS provides a abstract contextualization mechanism using ``DataContext`` to provide detailed error +information on the source of the data, which can be: + +- the file & line of the simulation XML file, +- if no source exists, ``Group`` / ``Wrapper`` path in the data repository. + +**Add ``DataContext`` when** an error occurs because of: + +- a ``Wrapper`` when validation fails for a given user input setting (``getWrapperDataContext()``), +- a ``Group`` when the error is implied by the instance state (``getDataContext()``), +- multiple data-contexts when the error is due to multiple objects / settings (first = more important). + +**Use ``getDataContext()`` or ``getWrapperDataContext()`` as last parameters of the logging macros to pass the data context metadata.** +It will be provided properly in the error message depending on the format (log / YAML). + +.. dropdown:: Example: Using DataContext + :icon: code + + .. code-block:: c++ + + // Example 1, using the getWrapperDataContext() because the error is + // related to a precise input setting. + GEOS_ERROR_IF( !meshBodies.hasGroup( meshBodyName ), + GEOS_FMT( "MeshBody ({}) is specified in targetRegions, but does not exist.", + meshBodyName ), + getWrapperDataContext( viewKeyStruct::targetRegionsString() ) ); + + // Exemple 2, using getDataContext() because the error is related with the entire instance + // state (a physics solver for instance) + GEOS_ERROR_IF( convergenceFailed, + "Failed to converge!", + getDataContext() ); + + // Exemple 3, using multiple DataContext because the error is related with multiple objects + // The first ones are considered the highest priority. + GEOS_THROW_IF( m_isThermal && !isFluidModelThermal, + GEOS_FMT( "CompositionalMultiphaseBase {}: the thermal option is enabled in the solver, but the fluid model {} is incompatible with the thermal option", + getDataContext(), fluid.getDataContext() ), + InputError, getDataContext(), fluid.getDataContext() ); + +Parallelism +=========== + +RAJA Kernels / CHAI Memory +--------------------------- + +Do Not Launch Kernels from Within Kernels +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +**Never launch RAJA kernels from within other RAJA kernels.** This causes portability issues and is not supported on GPU. + +.. dropdown:: Example: Nested kernel anti-pattern + :icon: code + + .. code-block:: c++ + + // BAD - Do not nest kernels + forAll>( n, [=] GEOS_HOST_DEVICE ( localIndex i ) + { + forAll>( m, [=] GEOS_HOST_DEVICE ( localIndex j ) + { + // ... + }); + }); + + // GOOD - Use nested loops or restructure algorithm + forAll>( n, [=] GEOS_HOST_DEVICE ( localIndex i ) + { + for( localIndex j = 0; j < m; ++j ) + { + // ... + } + }); + +CHAI Memory Management +^^^^^^^^^^^^^^^^^^^^^^ + +GEOS uses CHAI for **automatic kernel memory migration**. Follow these rules: + +- **Use CHAI-managed arrays for memory shared with kernels:** ``arrayView1d``, ``arraySlice2d``, ``array1d``, ``stackArray2d``, etc. + - ``array*`` types are designed for **host dynamic allocations**, + - ``arrayView*`` / ``arraySlice*`` types are for **non-reallocating scenario** (lambda capture, function passing), + - ``arrayView*`` types also serve for **automatic host-to-device memory move**, + - ``stackArray*`` types are used **to directly pass / store full arrays values**, use with care, +- **Call** ``.move(space)`` **to explicitly control memory location** +- **Never use manual CUDA calls** - RAJA has wrappers for device operations (compatibility layer for all platforms). + +.. dropdown:: Example: CHAI memory management + :icon: code + + .. code-block:: c++ + + array1d data( 1000 ); + + // Explicitly move to device before kernel + // (would have been done automatically when giving an LvArray view to a RAJA kernel) + data.move( parallelDeviceMemorySpace ); + + // Use in kernel + auto dataView = data.toView(); + forAll>( numElems, [=] GEOS_HOST_DEVICE ( localIndex i ) + { + dataView[i] = computeValue( i ); + }); + + // Move back to host if needed for CPU operations (costly) + data.move( hostMemorySpace ); + +Use RAJA Reductions +^^^^^^^^^^^^^^^^^^^ + +Use RAJA reductions for summation, min, max, and other parallel operations: + +.. dropdown:: Example: RAJA reductions + :icon: code + + .. code-block:: c++ + + RAJA::ReduceSum totalSum( 0.0 ); + RAJA::ReduceMin minValue( std::numeric_limits::max() ); + RAJA::ReduceMax maxValue( std::numeric_limits::lowest() ); + + forAll>( n, [=] GEOS_HOST_DEVICE ( localIndex i ) + { + totalSum += values[i]; + minValue.min( values[i] ); + maxValue.max( values[i] ); + }); + + real64 sum = totalSum.get(); + real64 min = minValue.get(); + real64 max = maxValue.get(); + +MPI Communication +----------------- + +Ensure All Ranks Participate in Collective Operations +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +**All MPI ranks must participate in collective operations.** Beware of branch deadlocks. + +.. dropdown:: Example: MPI collective operations + :icon: code + + .. code-block:: c++ + + // BAD - only rank 0 calls MPI_AllReduce() - DEADLOCK! + if( MpiWrapper::commRank() == 0 ) + { + MpiWrapper::allReduce( ... ); + } + + // BAD - some ranks may have this condition to false - not every ranks call MPI_Allreduce - DEADLOCK! + if( localWellElementsCount > 0 ) + { + MpiWrapper::allReduce( ... ); // WRONG - other ranks waiting! + } + + // GOOD - All ranks participate, even if they may not have influence in the result + real64 localSum = computeLocalSum(); + real64 globalSum = MpiWrapper::allReduce( localSum ); + + // GOOD - Conditional execution AFTER collective operation + real64 globalSum = MpiWrapper::allReduce( localSum ); + if( MpiWrapper::commRank() == 0 ) + { + GEOS_LOG_RANK_0( GEOS_FMT("Global sum = {}", globalSum) ); + } + +Always Use MpiWrapper +^^^^^^^^^^^^^^^^^^^^^ + +**Always use ``MpiWrapper`` instead of raw MPI calls.** + +Refer to ``MpiWrapper.hpp`` for more. + +.. dropdown:: Example: Using MpiWrapper + :icon: code + + .. code-block:: c++ + + // BAD - Raw MPI calls + int rank; + MPI_Comm_rank( MPI_COMM_WORLD, &rank ); + + // GOOD - Use MpiWrapper + int rank = MpiWrapper::commRank(); + int size = MpiWrapper::commSize(); + MpiWrapper::barrier(); + + // Common patterns + real64 globalMax = MpiWrapper::max( localValue ); + real64 globalSum = MpiWrapper::sum( localValue ); + MpiWrapper::allGather( sendBuf, recvBuf, count ); + +Communication Batching +^^^^^^^^^^^^^^^^^^^^^^^ + +**Batch small communications as much as possible**. + +.. dropdown:: Why? + :icon: info + + As MPI has significant per-message overhead, combining multiple small messages into larger ones can reduce communication overhead. + +.. dropdown:: Example: Batching MPI communication + :icon: code + + .. code-block:: cpp + + // WRONG - many small communications + void sendDataInefficiently( arrayView1d< real64 const > const & data, + int const destination ) + { + for( localIndex i = 0; i < data.size(); ++i ) + { + MpiWrapper::send( &data[i], 1, destination, tag + i ); // N messages! + } + } + + // CORRECT - single batched communication + void sendDataEfficiently( arrayView1d< real64 const > const & data, + int const destination ) + { + MpiWrapper::send( data.data(), data.size(), destination, tag ); // 1 message + } + +Performance +=========== + +Profiling +--------- + +Profile When Affecting performance-critical code +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +When modifying performance-critical code (kernels, assembly loops, solvers...): + +1. **Use Caliper profiling** on representative test cases before changes +2. **Measure performance** after changes +3. **Document performance impact** in pull request description + +.. dropdown:: Example: Performance profiling workflow + :icon: code + + .. code-block:: bash + + # Profile before changes + geos -i test_case.xml --caliper-output=baseline.cali + + # Make code changes + + # Profile after changes + geos -i test_case.xml --caliper-output=modified.cali + + # Compare results + cali-query -q "SELECT time.duration WHERE annotation='myKernel'" baseline.cali + cali-query -q "SELECT time.duration WHERE annotation='myKernel'" modified.cali + +Caliper Integration +^^^^^^^^^^^^^^^^^^^^ + +Use ``GEOS_MARK_FUNCTION``, ``GEOS_CALIPER_MARK_SCOPE`` and ``Timer`` for performance tracking for the main computation functions / scopes. + +.. dropdown:: Example: Performance instrumentation + :icon: code + + .. code-block:: c++ + + void PhysicsSolverBase::solverStep( ... ) + { + GEOS_MARK_FUNCTION; // Caliper tracking for entire function + + { + GEOS_CALIPER_MARK_SCOPE("assemble system"); + Timer timer( m_timers.get_inserted( "assembly" ) ); + assembleSystem(); + } + + { + GEOS_CALIPER_MARK_SCOPE("linear solve"); + Timer timer( m_timers.get_inserted( "linear solver" ) ); + linearSolve(); + } + } + +Data Repository / Validation +============================ + +Data Repository Paths +--------------------- + +Note that ``Group`` / ``Wrapper`` names & paths in the data repository are considered as magic values: +- **use getters rather than arbitrary names** (use view keys, ``getCatalogName()``, ``getName()`` or create custom ones). +- To avoid coupling / architecture stiffening, **avoid using full raw ``Group`` path**, especially when the path has parts unrelated to the component consulting it, prefer getters, + +Data validation +-------------------------- + +When To Validate +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +**Perform cross-parameter validation as soon as possible**, and when possible **only once**. +This can be done in the ``postInputInitialization()`` method, or later if needed, (i.e. in the +``registerDataOnMesh()`` when mesh data is needed, or during the simulation with lazy-init). + +.. dropdown:: Why doing validations at the ``postInputInitialization()`` stage? + :icon: info + + - All XML input has been parsed + - Default values have been applied + - Cross-references can be resolved + - User can receive clear error messages before simulation starts. + +.. dropdown:: Example: Post-input validation + :icon: code + + .. code-block:: c++ + + void MyClass::postInputInitialization() + { + // Call parent first + Base::postInputInitialization(); + + // Validate individual parameters + GEOS_THROW_IF_LT( m_timeStep, 0.0, + "timeStep must be positive", + InputError, getDataContext() ); + + // Validate cross-parameter constraints + GEOS_THROW_IF( m_endTime <= m_startTime, + "endTime must be greater than startTime", + InputError, getDataContext() ); + + // Validate against other components + GEOS_THROW_IF( !hasRequiredSolver(), + "Required solver not found in configuration", + InputError, getDataContext() ); + } + +Wrapped Data Validation +^^^^^^^^^^^^^^^^^^^^^^^ + +For data that is wrapped with a ``Wrapper`` (input / output of GEOS): + +- Use ``setInputFlag()`` to mark parameters as ``REQUIRED`` or ``OPTIONAL`` +- Use ``setApplyDefaultValue()`` to provide sensible defaults +- Use ``setRTTypeName()`` for runtime type validation (e.g., ``rtTypes::CustomTypes::positive``) +- Document valid ranges in ``setDescription()`` + +Physics-Specific Rules +---------------------- + +Unit Consistency +^^^^^^^^^^^^^^^^ + +- **Verify physical units consistancy** in your computations, +- **By default, values are in S.I. units**, refer to ``Units.hpp`` to see default metrics / conversion functions, +- **document any non-S.I. units** (variable name, comments). +- **Non-S.I. units are only for internal computation use.** + +Physical Bounds Checking +^^^^^^^^^^^^^^^^^^^^^^^^ + +Unphysical values indicate errors and can cause solver failures. +**Validate that state variables remain within physically meaningful bounds**. + +If a value is not strictly disallowed but does not seem possible for the model, **show a warning to the user that he can disable**. + +.. dropdown:: Example: Avoiding global state + :icon: code + + .. code-block:: cpp + + class ConstitutiveModel + { + void updateState( real64 const pressure, + real64 const temperature, + real64 const porosity ) + { + // Check physical bounds + GEOS_ERROR_IF( pressure < 1e-5, + GEOS_FMT( "Pressure {:.2e} Pa below cavitation limit", pressure ) ); + + GEOS_ERROR_IF( temperature < 0.0, + GEOS_FMT( "Absolute temperature {:.2f} K is negative", temperature ) ); + + GEOS_ERROR_IF( porosity < 0.0 || porosity > 1.0, + GEOS_FMT( "Porosity {:.3f} outside [0,1]", porosity ) ); + + // Warn about unusual but valid values, allow user to disable warning + GEOS_WARNING_IF( isHighTemperatureWarningEnabled && temperature > 1000.0, + GEOS_FMT( "High temperature {:.0f} K may be outside model validity", + temperature ) ); + } + }; diff --git a/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst b/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst index 7b4f5deafd2..22887b26dcf 100644 --- a/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst +++ b/src/docs/sphinx/developerGuide/Contributing/CodeRules.rst @@ -1,3 +1,6 @@ +.. _CodeRules: + + ############################################################################### Code Rules ############################################################################### @@ -7,1469 +10,14 @@ Code Rules Introduction ============ -This document provides comprehensive coding rules and best practices for developing with GEOS. These rules complement the :doc:`CodeStyle` guidelines and focus on typing, error handling, parallelism, performance, and general architecture principles. - -All contributors should familiarize themselves with these rules to maintain code quality, consistency, and reliability across the GEOS codebase. - -Typing -====== - -Type Correspondence Table -------------------------- - -GEOS defines its own type system for **portability**, **clarity**, and **performance**. - -**Always use GEOS types** instead of standard C++ types. - -Basic Types -^^^^^^^^^^^ - -.. list-table:: Standard C++ vs GEOS Types - :header-rows: 1 - :widths: 30 30 40 - - * - Standard C++ Type - - GEOS Type - - Description - * - ``int`` - - ``integer`` - - Signed integer - * - ``std::size_t`` - - ``localIndex`` - - Local array indexing (MPI partition) - * - N/A - - ``globalIndex`` - - Global indexing (across MPI) - * - ``float`` - - ``real32`` - - 32-bit floating point - * - ``double`` - - ``real64`` - - 64-bit floating point - * - ``std::string`` - - ``string`` - - String type - * - ``std::string_view`` - - ``string_view`` - - Non-owning string view - -.. dropdown:: Example: Using GEOS basic types - :icon: code - - .. code-block:: c++ - - // BAD - Do not use standard types - int count = 0; - double value = 1.5; - - // GOOD - Use GEOS types - integer count = 0; - real64 value = 1.5; - localIndex elementIndex = 42; - -Array Types -^^^^^^^^^^^ - -GEOS provides multi-dimensional CHAI array types for managed host-device data. - -Refer to the rule of :ref:`CHAI Memory Management`. - -.. list-table:: GEOS Array Types - :header-rows: 1 - :widths: 35 65 - - * - GEOS Type - - Description - * - ``array1d`` - - 1D owning array - * - ``arrayView1d`` - - 1D non-owning view (modifiable) - * - ``arraySlice1d`` - - 1D slice (from multi-dim array) - * - ``array2d`` - - 2D owning array - * - ``arrayView2d`` - - 2D non-owning view - * - ``arraySlice2d`` - - 2D slice - * - ``array3d``, ``array4d``, ``array5d`` - - 3D/4D/5D owning arrays - * - ``arrayView3d``, etc. - - 3D/4D/5D views - * - ``SortedArray`` - - Sorted array with search capabilities - * - ``ArrayOfArrays`` - - Jagged 2D array (variable row sizes) - * - ``ArrayOfSets`` - - Array of sets (unique elements per row) - * - ``stackArray1d`` - - Stack-allocated array (max size N) - -.. dropdown:: Example: Using GEOS arrays - :icon: code - - .. code-block:: c++ - - // BAD - Do not use std::vector for data that can be packed or kernel addressable - std::vector values; - std::array fixedData; - - // GOOD - Use CHAI managed arrays for kernels - array1d values; - arrayView1d valuesView = values.toView(); - arrayView1d constView = values.toViewConst(); - stackArray1d fixedData; - -Tensor Types -^^^^^^^^^^^^ - -Use GEOS tensor types for geometric and mechanical calculations: - -.. list-table:: GEOS Tensor Types - :header-rows: 1 - :widths: 30 70 - - * - Type - - Description - * - ``R1Tensor`` - - 3D vector (real64, default choice over 32bit counterpart) - * - ``R1Tensor32`` - - 3D vector (real32) - * - ``R2SymTensor`` - - Symmetric 6-component tensor (`Voigt notation `_) - -``std`` Standard Library ------------------------- - -It is possible to use ``std`` library components for data on host memory, for doing so, the rule is -to **only use GEOS ``std`` container wrappers** instead of direct standard library containers. - -This rule allow us to control bounds checking depending on ``GEOS_USE_BOUNDS_CHECK`` macro / ``GEOS_ENABLE_BOUNDS_CHECK`` cmake option.. - - * - Standard C++ Type - - GEOS Type - - Description - * - ``std::vector`` - - ``stdVector`` - - ``std`` dynamically allocated array - * - ``std::array`` - - ``stdArray`` - - ``std`` fixed constexpr sized array - * - ``std::map`` - - ``map`` - - ``std`` sorted dictionary - * - ``std::unordered_map`` - - ``unordered_map`` - - ``std`` unsorted dictionary (hash-map) - -.. dropdown:: Example: Container usage - :icon: code - - .. code-block:: c++ - - // BAD - Direct std containers - std::vector data; - std::array fixedData; - std::map orderedMap; - - // GOOD - GEOS wrappers - stdVector data; - stdArray fixedData; - map orderedMap; - -The following standard library components are allowed: -- ``std::pair``, ``std::tuple`` for temporary structures, `but sometimes struct are to prefer `_. -- ``std::function`` for callbacks -- ``std::optional`` for optional return values -- ``std::move``, ``std::forward`` for move semantics -- Algorithms from ```` where appropriate -- ``std::unique_ptr``, ``std::shared_ptr`` for memory management - -Contribution -============ - -Documentation -------------- - -**All new features must include appropriate documentation, for developpers and user, at the following levels:** - -Wrapper Documentation (User-Oriented) -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -All data repository wrappers must be documented with ``setDescription()``. -As much as possible, valid values rules should be provided, and default values should be provided with ``setApplyDefaultValue()``. - -.. dropdown:: Example: Wrapper documentation - :icon: code - - .. code-block:: c++ - - registerWrapper( "timeStep", &m_timeStep ). - setInputFlag( InputFlags::REQUIRED ). - setApplyDefaultValue( 1.0e-6 ). - setDescription( "Initial time step value. This parameter controls " - "the starting time increment for the simulation. " - "Valid range: (0, 1e-3]. See User Guide Chapter 5." ); - -RST Documentation (User Guide) -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -User-facing features must be documented in the Sphinx RST documentation (``src/docs/sphinx/``): - -- Explain **what** the feature does -- Provide **usage examples** -- Document **input parameters** -- Include **expected output** or behaviour -- Add **references** to related features - -Doxygen & Naming (Developer-Oriented) -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -Keep Doxygen comments and naming as clear as possible. Variable names and inline comments should provide help for developers even if they are not experts in the domain. - -.. dropdown:: Example: Doxygen documentation - :icon: code - - .. code-block:: c++ - // GOOD - The documentation gives useful information about what the function *does*, without - // going into theorical knowledge nor implementation details. - /** - * @brief Computes the residual for the nonlinear solver. - * @param[in] domain The domain partition containing mesh and field data - * @param[in] time Current simulation time - * @param[in] dt Time step size - * @param[out] residual The computed residual vector - * @return Residual L2 norm - * - * This method assembles the discrete residual for the physics equations - * by looping over all elements in the target regions. The residual is - * stored in the provided vector and its L2 norm is returned. - */ - real64 computeResidual( DomainPartition & domain, - real64 const time, - real64 const dt, - arrayView1d const & residual ); - - // BAD - documentation does not provide useful information beyond the function name. - /** - * @brief Compute CFL numbers - * @param[in] time current time - * @param[in] dt the time step size - * @param[in] domain the domain partition - */ - void computeCFLNumbers( real64 const time, - real64 const dt, - DomainPartition & domain ) const; - - -Testing -------- - -Unit Tests -^^^^^^^^^^ - -**Always provide unit tests** to ensure every function behaves according to expectations. -Unit testing is not about testing every single line of code (an unrealistic goal), but about *verifying assumptions and preventing regressions*. Do not assume that your code will always work as intended — protect it with unit tests. - -- Focus on **critical logic** (core functionality and algorithms), **edge cases** (extreme, erroneous, empty values). -- **Use GoogleTest framework**. -- In GEOS, the goal of unit test is never to take position on the validity of physical results (this is the point of integrated-tests). -- Place tests in ``src/coreComponents//tests/``. -- Test GPU code paths if applicable (use ``#ifdef GEOS_USE_DEVICE``). - -.. dropdown:: Example: Unit test structure - :icon: code - - .. code-block:: c++ - - #include - - // Test core functionality - TEST( MyComponentTest, ComputesCorrectResult ) - { - MyComponent component; - EXPECT_DOUBLE_EQ( component.compute( 2.0 ), 4.0 ); - } - - // Test edge case - TEST( MyComponentTest, HandlesZeroInput ) - { - MyComponent component; - EXPECT_DOUBLE_EQ( component.compute( 0.0 ), 0.0 ); - } - - // Test error detection - TEST( MyComponentTest, RejectsInvalidInput ) - { - MyComponent component; - EXPECT_THROW( component.compute( -1.0 ), std::invalid_argument ); - } - - // Test the `testMyFunc()` function for GPU (geos::parallelDevicePolicy, if applicable) and - // CPU scenario (serialPolicy) - #ifdef GEOS_USE_DEVICE - TEST( MyFuncTests, testMyFuncDevice ) - { - testMyFunc< geos::parallelDevicePolicy< > >(); - } - #endif - TEST( MyFuncTests, testMyFuncHost ) - { - testMyFunc< serialPolicy >(); - } - - // entry point of the test binary - int main( int argc, char * * argv ) - { - ::testing::InitGoogleTest( &argc, argv ); - g_commandLineOptions = *geos::basicSetup( argc, argv ); - int const result = RUN_ALL_TESTS(); - geos::basicCleanup(); - return result; - } - -Code Coverage -^^^^^^^^^^^^^ - -**Code coverage should never decrease.** New code contributions must maintain or improve overall code coverage. -Use code-cov to report untested code paths. - -Integrated Tests & Examples -^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**Always provide an example or an integrated test for every XML-usable & serializable component:** - -- Each usage pattern should have its own dedicated example, -- Examples must run as part of the integrated tests to protect features against regressions, -- Place examples in ``examples/`` directory, and integrated tests in ``inputFiles/``, - -For more info, please refer to :ref:`IntegratedTests.rst`. - -Logging -======= - -Rules ------ - -GEOS Logging Macros (``GEOS_LOG*``) -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**Use of ``std::cout``, ``std::cerr``, or ``printf`` for logging must never appear** on the develop branch. -Always use GEOS logging macros defined in ``Logger.hpp``: - -.. dropdown:: Example: Correct logging usage - :icon: code - - .. code-block:: c++ - - // BAD - Do not use - std::cout << "Starting simulation" << std::endl; - printf("Value = %f\n", value); - - // GOOD - Use GEOS macros - GEOS_LOG("Starting simulation"); - GEOS_LOG_RANK_0("Master rank initialization complete"); - -**Available logging macros:** - -- ``GEOS_LOG(...)`` - Log on all ranks -- ``GEOS_LOG_VAR(...)`` - Log variable name and value -- ``GEOS_LOG_IF(condition, msg)`` - Conditional logging -- ``GEOS_LOG_RANK_0(msg)`` - Log only on rank 0 (MPI master) -- ``GEOS_LOG_RANK_0_IF(condition, msg)`` - Conditional logging on rank 0 -- ``GEOS_LOG_RANK(msg)`` - Log to rank-specific output stream -- ``GEOS_LOG_RANK_VAR(var)`` - Log variable on rank stream - -Using Log Levels -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**Logging should be meaningful and controlled**: - -- When possible, **use ``GEOS_LOG_LEVEL_RANK*`` macros with appropriate ``logInfo``** to allow output filtering, -- **Avoid unnecessary logging**, prefer rank 0 logging for global information to avoid redundant output, -- **Consider logs performance impact** (output flush frequency, formatting). - -See :doc:`LogLevels` for using / adding log level (e.g., ``logInfo::Convergence``, ``logInfo::TimeStep``). - -.. dropdown:: Example: Using log levels - :icon: code - - .. code-block:: c++ - - // physics solver cpp - effective logging - GEOS_LOG_LEVEL_RANK_0( logInfo::Convergence, - GEOS_FMT( "Newton iteration {}: residual = {:4.2e}", - iter, residual ) ); - - // src/coreComponents/physicsSolvers/LogLevelsInfo.hpp - convergence logging config - struct Convergence - { - static constexpr int getMinLogLevel() { return 1; } - static constexpr std::string_view getDescription() { return "Convergence information"; } - }; - -.. dropdown:: Example: Avoiding log spam - :icon: code - - .. code-block:: c++ - - // BAD - Will spam output for every element on every ranks - forAll< parallelDevicePolicy<> >( numElems, [=] GEOS_HOST_DEVICE ( localIndex const ei ) - { - GEOS_LOG("Processing element " << ei); // DO NOT DO THIS - }); - - // GOOD - Log summary information - GEOS_LOG_RANK_0(GEOS_FMT("Processed {} elements", numElems)); - -Logging in RAJA Kernels -^^^^^^^^^^^^^^^^^^^^^^^^ - -**Logging inside RAJA kernels should only be done for debugging purposes** and must never appear on the develop branch. - -.. dropdown:: Why logging on GPU is not allowed on production code? - :icon: info - - Logging on GPU can: - - - Reserve cache lines and degrade performance, - - Cause race conditions in output, - - Generate massive amounts of output, - - Produce unpredictable behaviour (e.g. runtime crashes on AMD). - -.. dropdown:: Example: Debug-only kernel logging - :icon: code - - .. code-block:: c++ - - forAll< parallelDevicePolicy<> >( numElems, [=] GEOS_HOST_DEVICE ( localIndex const ei ) - { - GEOS_LOG("processing..."); // Compiles, but remove before merge to develop! - GEOS_LOG(GEOS_FMT("processing elem {}/{}", ei, numElems)); // Will crash on CUDA builds (unsupported GEOS_FMT) - printf("DEBUG: ei=%d\n", ei); // Compiles, but remove before merge to develop! - }); - -Structured Data Output ----------------------- - -For tabular or statistical output, use structured logging facilities instead of plain text logging: - -- Use ``TableTextFormatter`` with a ``logLevel`` XML setting for formatted text table output, -- Use ``TableCSVFormatter`` with a ``writeCSV()`` XML setting for structured data output. - -.. dropdown:: Example: Log table output - :icon: code - - .. code-block:: c++ - // log output preparation, layout precomputation can be kept for periodic outputs - TableLayout statsLogLayout( "", { "Flux(es)", "Region", "Element Count", massColumn, rateColumn } ); - m_logLayout = statsLogLayout; - - // filling table data (log format is for current timestep) - m_logData.addRow( fluxName, - elementSetName, - GEOS_FMT( "{}", wrappedStats.stats().m_elementCount ), - GEOS_FMT( "{}", wrappedStats.stats().m_producedMass ), - GEOS_FMT( "{}", wrappedStats.stats().m_productionRate ) ); - - // statistics CSV output (only if enabled and from rank 0) - if( logLevelActive && MpiWrapper::commRank() == 0 ) - { - string const title = GEOS_FMT( "{}, flux statistics for: {}", - getName(), fluxesStr ); - string const multilineTitle = wrapTextToMaxLength( title, 80 ); // enabling automatic line wrapping - m_logLayout.setTitle( multilineTitle ); - TableTextFormatter const tableStatFormatter( m_logLayout ); - GEOS_LOG_RANK( tableStatFormatter.toString( statsData ) ); - } - -.. dropdown:: Example: CSV output - :icon: code - - .. code-block:: c++ - - // CSV output preparation, layout precomputation can be kept for periodic outputs - TableLayout const statsCSVLayout( "", {"Time [s]", "Flux(es)", "Region", "Element Count", massColumn, rateColumn} ); - m_csvLayout = statsCSVLayout; - - // filling table data (CSV format is all timestep series) - m_csvData.addRow( wrappedStats.getStatsPeriodStart(), - fluxName, - elementSetName, - wrappedStats.stats().m_elementCount, - wrappedStats.stats().m_producedMass, - wrappedStats.stats().m_productionRate ); - - // statistics CSV output (only if enabled and from rank 0) - if( m_writeCSV > 0 && MpiWrapper::commRank() == 0 ) - { - std::ofstream outputFile( m_csvFilename ); - TableCSVFormatter const tableStatFormatter( m_csvLayout ); - outputFile << tableStatFormatter.toString( csvData ); - outputFile.close(); - csvData.clear(); - } - -Error Management -================ - -- **Developpers must provide runtime validations** of their code behaviour, -- For runtime validation fails that can be triggered by users, **meaningful messages and context must be provided** to help troubleshooting. - -Errors (GEOS_ERROR*) --------------------- - -Why Use Errors? -^^^^^^^^^^^^^^^ - -Use ``GEOS_ERROR*`` macros when encountering a **blocking error** where: - -- The simulation result would be certifiably invalid, -- The application state is unrecoverable, -- Continuing execution would be unsafe or meaningless. - -.. dropdown:: Example: Error handling - :icon: code - - .. code-block:: c++ - - // Example: not implemented feature - GEOS_ERROR( "Rounding interpolation with derivatives not implemented" ); - - // Example: invalid mesh topology - GEOS_ERROR_IF( numFaces == 0, - GEOS_FMT("Element {} has zero faces - invalid mesh", elemId) ); - - // Example: comparison-based error - GEOS_ERROR_IF_LT( dt, 0.0, "Time step must be positive" ); - -Exceptions (GEOS_THROW*) ------------------------- - -Why Use Exceptions? -^^^^^^^^^^^^^^^^^^^ - -Use ``GEOS_THROW*`` macros for the same reasons as ``GEOS_ERROR*`` (**unrecoverable state**), and when: - -- You want to **add context** at higher call stack levels, -- You need to **propagate detailed error information** upward -- If you use custom exception class, you need to document them here. - -.. list-table:: Available exception types - :header-rows: 1 - :widths: 35 65 - - * - Error class - - Used for - * - ``std::runtime_error`` - - General runtime errors - * - ``geos::InputError`` - - Errors in user input (XML, data files) - * - ``geos::SimulationError`` - - Errors during simulation execution - * - ``geos::BadTypeError`` - - Type conversion errors - * - ``geos::NotAnError`` - - Non critical, volontary stopping state - -.. dropdown:: Example: Throwing exceptions - :icon: code - - .. code-block:: c++ - - // Example 1: throw with context - GEOS_THROW_IF( !file.is_open(), - GEOS_FMT("Cannot open file {}", filename), - InputError ); - - // Example 2: rethrow with context info (because of unclear exception) - try - { - tolerance = stod( inputParams[8] ); - } - catch( const std::invalid_argument & e ) - { - GEOS_THROW( GEOS_FMT( "{}: invalid model parameter value: {}", functionName, e.what() ), - InputError, getDataContext() ); - } - -Never Catch and Continue -^^^^^^^^^^^^^^^^^^^^^^^^^ - -**Do not catch exceptions and continue execution**. Exceptions in GEOS indicate serious problems, they -do not serve as an alternative code path to consult a behaviour state. - -Why not to use exceptions for anything else than passing information up to the call stack: - -- **Performance:** Exception handling is expensive, especially with stack unwinding, -- **Predictability:** Catching and continuing makes behaviour unpredictable and hard to debug, -- **Safety:** The error state may have corrupted simulation data, broken the code path at an unexpected place, and/or left the system in an invalid state. - -.. dropdown:: Example: Exception handling anti-pattern - :icon: code - - .. code-block:: c++ - - // BAD - Do not catch to continue - // Here, the "solver" should have a "isFailed" flag or a returned resulting state - try - { - solver.solve(); - } - catch( std::exception const & e ) // <- Note: this exception could have been throw for unexpected reason - { - GEOS_LOG("Solver failed, trying alternative"); - alternativeSolver.solve(); - } - - // GOOD - Manage your code paths manually - if(! solver.solve() ) - { - GEOS_LOG("Solver failed, trying alternative"); - alternativeSolver.solve(); - } - -Warnings (GEOS_WARNING*) ------------------------- - -Why Use Warnings? -^^^^^^^^^^^^^^^^^ - -Use ``GEOS_WARNING*`` macros when: - -- An issue is detected but **simulation can continue** -- Simulation may be **affected but not completly invalid** -- The user should be notified of a **potential problem** - -.. dropdown:: Example: Warning usage - :icon: code - - .. code-block:: c++ - - // Example: suboptimal but valid configuration - GEOS_WARNING_IF( dt > recommendedDt, - GEOS_FMT("Time step {} exceeds recommended value {}. This may affect solution accuracy.", - dt, recommendedDt), - getDataContext() ); - -Errors & Warning in RAJA Kernels --------------------------------- - -Similarly to Logging, **error and warning output inside RAJA kernels is only for debugging purposes** and must never appear on the develop branch. - -.. dropdown:: Why errors are forbidden on GPU in production code? - :icon: info - - - GPU kernel errors cause immediate kernel termination, adding potentially costly branches, - - Error handling on device has significant performance and cache impact, - - GPU error handling may be unsupported depending on the platform / device, - - Can cause deadlocks in parallel execution (`as discussed here`_, a rank can die, then the other ranks will wait for it at next MPI call). - -.. dropdown:: Example: Kernel error handling (debug only) - :icon: code - - .. code-block:: c++ - - // On GPU, errors will cause kernel abort - forAll< parallelDevicePolicy<> >( numElems, [=] GEOS_HOST_DEVICE ( localIndex const ei ) - { - #if !defined( NDEBUG ) - // This will trap on GPU - use sparingly for debugging only - GEOS_ERROR_IF( badCondition, "Bad condition detected" ); - #endif - }); - -Contextualization with DataContext ------------------------------------ - -GEOS provides a abstract contextualization mechanism using ``DataContext`` to provide detailed error -information on the source of the data, which can be: - -- the file & line of the simulation XML file, -- if no source exists, ``Group`` / ``Wrapper`` path in the data repository. - -**Add ``DataContext`` when** an error occurs because of: - -- a ``Wrapper`` when validation fails for a given user input setting (``getWrapperDataContext()``), -- a ``Group`` when the error is implied by the instance state (``getDataContext()``), -- multiple data-contexts when the error is due to multiple objects / settings (first = more important). - -**Use ``getDataContext()`` or ``getWrapperDataContext()`` as last parameters of the logging macros to pass the data context metadata.** -It will be provided properly in the error message depending on the format (log / YAML). - -.. dropdown:: Example: Using DataContext - :icon: code - - .. code-block:: c++ - - // Example 1, using the getWrapperDataContext() because the error is - // related to a precise input setting. - GEOS_ERROR_IF( !meshBodies.hasGroup( meshBodyName ), - GEOS_FMT( "MeshBody ({}) is specified in targetRegions, but does not exist.", - meshBodyName ), - getWrapperDataContext( viewKeyStruct::targetRegionsString() ) ); - - // Exemple 2, using getDataContext() because the error is related with the entire instance - // state (a physics solver for instance) - GEOS_ERROR_IF( convergenceFailed, - "Failed to converge!", - getDataContext() ); - - // Exemple 3, using multiple DataContext because the error is related with multiple objects - // The first ones are considered the highest priority. - GEOS_THROW_IF( m_isThermal && !isFluidModelThermal, - GEOS_FMT( "CompositionalMultiphaseBase {}: the thermal option is enabled in the solver, but the fluid model {} is incompatible with the thermal option", - getDataContext(), fluid.getDataContext() ), - InputError, getDataContext(), fluid.getDataContext() ); - -Parallelism -=========== -RAJA Kernels / CHAI Memory ---------------------------- +This section provides comprehensive coding rules and best practices for developing with GEOS. These rules complement the :doc:`CodeStyle` guidelines and focus on typing, error handling, parallelism, performance, and general architecture principles. -Do Not Launch Kernels from Within Kernels -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**Never launch RAJA kernels from within other RAJA kernels.** This causes portability issues and is not supported on GPU. - -.. dropdown:: Example: Nested kernel anti-pattern - :icon: code - - .. code-block:: c++ - - // BAD - Do not nest kernels - forAll>( n, [=] GEOS_HOST_DEVICE ( localIndex i ) - { - forAll>( m, [=] GEOS_HOST_DEVICE ( localIndex j ) - { - // ... - }); - }); - - // GOOD - Use nested loops or restructure algorithm - forAll>( n, [=] GEOS_HOST_DEVICE ( localIndex i ) - { - for( localIndex j = 0; j < m; ++j ) - { - // ... - } - }); - -CHAI Memory Management -^^^^^^^^^^^^^^^^^^^^^^ - -GEOS uses CHAI for **automatic kernel memory migration**. Follow these rules: - -- **Use CHAI-managed arrays:** ``array1d``, ``array2d``, etc. -- **Call** ``.move(space)`` **to explicitly control memory location** -- **Never use manual CUDA calls** - RAJA has wrappers for device operations - -.. dropdown:: Example: CHAI memory management - :icon: code - - .. code-block:: c++ - - array1d data( 1000 ); - - // Explicitly move to device before kernel - // (would have been done automatically when giving an LvArray view to a RAJA kernel) - data.move( parallelDeviceMemorySpace ); - - // Use in kernel - auto dataView = data.toView(); - forAll>( numElems, [=] GEOS_HOST_DEVICE ( localIndex i ) - { - dataView[i] = computeValue( i ); - }); - - // Move back to host if needed for CPU operations (costly) - data.move( hostMemorySpace ); - -Use RAJA Reductions -^^^^^^^^^^^^^^^^^^^ - -Use RAJA reductions for summation, min, max, and other parallel operations: - -.. dropdown:: Example: RAJA reductions - :icon: code - - .. code-block:: c++ - - RAJA::ReduceSum totalSum( 0.0 ); - RAJA::ReduceMin minValue( std::numeric_limits::max() ); - RAJA::ReduceMax maxValue( std::numeric_limits::lowest() ); - - forAll>( n, [=] GEOS_HOST_DEVICE ( localIndex i ) - { - totalSum += values[i]; - minValue.min( values[i] ); - maxValue.max( values[i] ); - }); - - real64 sum = totalSum.get(); - real64 min = minValue.get(); - real64 max = maxValue.get(); - -MPI Communication ------------------ - -Ensure All Ranks Participate in Collective Operations -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**All MPI ranks must participate in collective operations.** Beware of branch deadlocks. - -.. dropdown:: Example: MPI collective operations - :icon: code - - .. code-block:: c++ - - // BAD - only rank 0 calls MPI_AllReduce() - DEADLOCK! - if( MpiWrapper::commRank() == 0 ) - { - MpiWrapper::allReduce( ... ); - } - - // BAD - some ranks may have this condition to false - not every ranks call MPI_Allreduce - DEADLOCK! - if( localWellElementsCount > 0 ) - { - MpiWrapper::allReduce( ... ); // WRONG - other ranks waiting! - } - - // GOOD - All ranks participate, even if they may not have influence in the result - real64 localSum = computeLocalSum(); - real64 globalSum = MpiWrapper::allReduce( localSum ); - - // GOOD - Conditional execution AFTER collective operation - real64 globalSum = MpiWrapper::allReduce( localSum ); - if( MpiWrapper::commRank() == 0 ) - { - GEOS_LOG_RANK_0( GEOS_FMT("Global sum = {}", globalSum) ); - } - -Always Use MpiWrapper -^^^^^^^^^^^^^^^^^^^^^ - -**Always use ``MpiWrapper`` instead of raw MPI calls.** - -Refer to ``MpiWrapper.hpp`` for more. - -.. dropdown:: Example: Using MpiWrapper - :icon: code - - .. code-block:: c++ - - // BAD - Raw MPI calls - int rank; - MPI_Comm_rank( MPI_COMM_WORLD, &rank ); - - // GOOD - Use MpiWrapper - int rank = MpiWrapper::commRank(); - int size = MpiWrapper::commSize(); - MpiWrapper::barrier(); - - // Common patterns - real64 globalMax = MpiWrapper::max( localValue ); - real64 globalSum = MpiWrapper::sum( localValue ); - MpiWrapper::allGather( sendBuf, recvBuf, count ); - -Communication Batching -^^^^^^^^^^^^^^^^^^^^^^^ - -**Batch small communications as much as possible**. - -.. dropdown:: Why? - :icon: info - - As MPI has significant per-message overhead, combining multiple small messages into larger ones can reduce communication overhead. - -.. dropdown:: Example: Batching MPI communication - :icon: code - - .. code-block:: cpp - - // WRONG - many small communications - void sendDataInefficiently( arrayView1d< real64 const > const & data, - int const destination ) - { - for( localIndex i = 0; i < data.size(); ++i ) - { - MpiWrapper::send( &data[i], 1, destination, tag + i ); // N messages! - } - } - - // CORRECT - single batched communication - void sendDataEfficiently( arrayView1d< real64 const > const & data, - int const destination ) - { - MpiWrapper::send( data.data(), data.size(), destination, tag ); // 1 message - } - - -Validation / Precision -====================== - -Computation Validation ------------------------ - -Use Tolerance-Based Comparisons -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**Always consider proper tolerance** for floating-point numbers comparisons, taking into account rounding errors, even for extreme values. - -.. dropdown:: Example: Correct float comparison - :icon: code - - .. code-block:: c++ - - // BAD - Direct comparison - if( value == 1.0 ) { ... } - if( a == b ) { ... } - - // GOOD - Absolute tolerance - real64 const absTol = 1.0e-12; - if( LvArray::math::abs(value) < absTol ) { ... } - - // GOOD - Relative tolerance - real64 const relTol = 1.0e-8; - real64 const scale = std::max( std::abs(a), std::abs(b) ); - if( std::abs(a - b) < relTol * scale ) { ... } - -Division Safety, NaN/Inf Values -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -- **Always verify that denominators are not zero or near-zero before a division**. -- In General, **we should not make any computation that could result in NaN/Inf values**. - -.. dropdown:: Example: Division Safety - :icon: code - - .. code-block:: cpp - - // WRONG - unprotected division - real64 const normalizedResidual = residual / initialResidual; - real64 const strainRate = velocityGradient / thickness; - - // CORRECT - protected division - real64 computeNormalizedResidual( real64 const residual, - real64 const initialResidual ) - { - if( initialResidual > machinePrecision ) - return residual / initialResidual; - else - return residual; // or return a flag indicating special case - } - - // CORRECT - with error reporting - real64 safeDivide( real64 const numerator, - real64 const denominator, - string const & context ) - { - GEOS_ERROR_IF( LvArray::math::abs( denominator ) < machinePrecision, - GEOS_FMT( "Division by zero in {}: denominator = {:.2e}", - context, denominator ) ); - return numerator / denominator; - } - -Overflow Prevention -^^^^^^^^^^^^^^^^^^^^ - -Overflow leads to undefined behavior and memory corruption. **Validate that operations won't exceed type / container limits**, especially for index calculations. - -Use LvArray Math Utilities -^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -The ``LvArray::math`` namespace provides safe numerical utilities, avoid LvArray math utilities instead of platform specific math function. - -.. dropdown:: Example: LvArray math utilities - :icon: code - - .. code-block:: c++ - - #include "LvArray/src/math.hpp" - - real64 const absValue = LvArray::math::abs( value ); - real64 const maxValue = LvArray::math::max( a, b ); - real64 const minValue = LvArray::math::min( a, b ); - -Data repository validation -------------------------------- - -Validate in postInputInitialization -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**Perform cross-parameter validation as soon as possible**, and when possible **only once **. -This can be done in the ``postInputInitialization()`` method, or later if needed, (i.e. in the -``registerDataOnMesh()`` or during the simulation). - -At the ``postInputInitialization()`` stage: - -- All XML input has been parsed -- Default values have been applied -- Cross-references can be resolved -- User can receive clear error messages before simulation starts. - -.. dropdown:: Example: Post-input validation - :icon: code - - .. code-block:: c++ - - void MyClass::postInputInitialization() - { - // Call parent first - Base::postInputInitialization(); - - // Validate individual parameters - GEOS_THROW_IF_LT( m_timeStep, 0.0, - "timeStep must be positive", - InputError, getDataContext() ); - - // Validate cross-parameter constraints - GEOS_THROW_IF( m_endTime <= m_startTime, - "endTime must be greater than startTime", - InputError, getDataContext() ); - - // Validate against other components - GEOS_THROW_IF( !hasRequiredSolver(), - "Required solver not found in configuration", - InputError, getDataContext() ); - } - -Additional Validation Guidelines -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -- Use ``setInputFlag()`` to mark parameters as ``REQUIRED`` or ``OPTIONAL`` -- Use ``setApplyDefaultValue()`` to provide sensible defaults -- Use ``setRTTypeName()`` for runtime type validation (e.g., ``rtTypes::CustomTypes::positive``) -- Document valid ranges in ``setDescription()`` - -Performance -=========== - -Profiling ---------- - -Profile When Affecting Kernel Behaviour -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -When modifying performance-critical code (kernels, assembly loops, solvers): - -1. **Use Caliper profiling** on representative test cases before changes -2. **Measure performance** after changes -3. **Document performance impact** in pull request description - -.. dropdown:: Example: Performance profiling workflow - :icon: code - - .. code-block:: bash - - # Profile before changes - geos -i test_case.xml --caliper-output=baseline.cali - - # Make code changes - - # Profile after changes - geos -i test_case.xml --caliper-output=modified.cali - - # Compare results - cali-query -q "SELECT time.duration WHERE annotation='myKernel'" baseline.cali - cali-query -q "SELECT time.duration WHERE annotation='myKernel'" modified.cali - -Caliper Integration -^^^^^^^^^^^^^^^^^^^^ - -Use ``GEOS_MARK_FUNCTION``, ``GEOS_CALIPER_MARK_SCOPE`` and ``Timer`` for performance tracking for the main computation functions / scopes. - -.. dropdown:: Example: Performance instrumentation - :icon: code - - .. code-block:: c++ - - void PhysicsSolverBase::solverStep( ... ) - { - GEOS_MARK_FUNCTION; // Caliper tracking for entire function - - { - GEOS_CALIPER_MARK_SCOPE("assemble system"); - Timer timer( m_timers.get_inserted( "assembly" ) ); - assembleSystem(); - } - - { - GEOS_CALIPER_MARK_SCOPE("linear solve"); - Timer timer( m_timers.get_inserted( "linear solver" ) ); - linearSolve(); - } - } - -Speed Optimization Rules ------------------------------ - -- **Hoist Loop Invariants**, move computations that don't change during iterations outside the loop. -- When it does not critically affect the code architecture and clarity, **fuse multiple related kernels to reduce memory traffic and launch overhead** (i.e., statistics kernels process all physics field at once). -- **Optimize Memory Access for Cache and Coalescing**. Access memory sequentially and ensure coalesced access, especially on GPUs. -- **Minimize Host-Device Transfers**. Keep data on the appropriate memory space and minimize transfers. - -Memory Management Rules ---------------------------------- - -Generalities -^^^^^^^^^^^^^^^^^^^^ - -- **Minimize dynamic memory allocation** as much as possible, particularly in performance-critical code, -- For frequently allocated/deallocated objects, **consider memory pools**, -- For containers with known size, **reserve capacity upfront**. - -Be Const / Constexpr -^^^^^^^^^^^^^^^^^^^^ - -**Use ``const`` and ``constexpr`` when possible**, enabling: - -- **Compiler optimization,** enables better code generation by giving constraints to the code, -- **Code safety,** prevents accidental modification for constant contexts, -- **Show code intention,** make code more readable. - -Also, **mark methods const if the method is not designed to modify** the object state. - -Value / Const Function Arguments -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**Pass large and non-trivially copyable objects by const reference**: - -- A "large" size can be defined as "more that 16 bytes, ", which is 2 pointers, 2 integer / index, or 2 double, -- An object is non-trivially copyable objects when it needs to perform sub-allocation when being copied, - -.. dropdown:: Value / Const Reference function parameters practices - :icon: code - - .. code-block:: c++ - - // 16 bytes (string_view are equivalent to 2 constant pointers): PASS BY VALUE - static constexpr string_view str0 = "Hello"; - string const str1 = getTimeStr(); // constant string are passed as string_view, beware of lifetime! - void func( string_view param ); - - arrayView2d< int > myDataView; - - // 16 bytes: PASS BY VALUE - stdArray< double, 2 > vec2D; - void func( vec2D param ); - - // 12 to 16 bytes (globalIndex is 8, integer is 4 to 8): PASS BY VALUE - struct SmallStruct - { - globalIndex id; - integer count; - }; - void func( SmallStruct param ); - - // 16 to 20 bytes (globalIndex is 8, integer is 4 to 8): PASS BY REFERENCE - struct BigStruct - { - globalIndex id; - integer countA; - integer countB; - }; - void func( BigStruct const & param ); - - // Does dynamic allocation: PASS BY REFERENCE - map< int, long > myMap; - stdVector< int > myList { 123 }; - void func( map< int, long > const & myMap, - stdVector< int > const & myList ); - - // LvArray types parameter practices depends on the intent - array2d< int > myConstantData; - array2d< int > myMutableData;, - // passing as non-view means that we will potencially resize the arrays: PASS BY REFERENCE - void arrayResizer( array2d< const int > & myConstantData - array2d< int > & myMutableData ); - // passing as view means that we may just affect the data, only if non-const template: PASS BY VALUE - void arrayProcess( arrayView2d< const int > myConstantData - arrayView2d< int > myMutableData ); - - // 16 bytes superficially... but does dynamic allocation: PASS BY REFERENCE - struct DynAllocStruct - { - string const a; - }; - void func( DynAllocStruct const & param ); - - // Any Group subclass is large: PASS BY REFERENCE - void modify( DomainPartition & domain ); - -Pointer / Reference Function Arguments -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -When passing nullable objects as parameter, pointers can be passed. - -**Pointer parameters must always be null-checked when dereferenced**. -If applying this rule produces repetitive code, consider using a const reference. - -This rule imply the following: -- **A reference may never be null**, -- In a given method, **``this`` may never be null**. - -.. dropdown:: Why? - :icon: info - - Passing a pointer or reference parameter show the intention of modifying the underlying object. - The difference of the two practices is that passing a pointer means that we consider that the - underlying object can be null. - -Constexpr for Compile-Time Constants -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**Use ``constexpr`` for values known at compile time**. - -The rule is not absolute: when the impact is not significative, and the code needs is getting really unclear, the rule can be ignored. - -.. dropdown:: Example: Constexpr usage - :icon: code - - .. code-block:: c++ - - // Compile-time constants - // A rule of thumb is that oftenly any time a value is constexpr, it can also be static. - static constexpr localIndex numDimensions = 3; - static constexpr real64 tolerance = 1.0e-10; - - // Constexpr functions (evaluated at compile time when possible) - constexpr localIndex computeArraySize( localIndex n ) - { return n * n + 2 * n + 1; } - -Provide Views to Arrays -^^^^^^^^^^^^^^^^^^^^^^^^ - -- When possible, prefer provide views to arrays (to const data if possible). -- This **must** be done especially for RAJA kernels, and views must be **captured by value**. - -The rule is generalizable to ``string_view`` for strings, but not applicable in GPU context. - -.. dropdown:: Why Use Views? - :icon: info - - - **No memory allocation:** Views are lightweight references - - **Mutability correctness:** Can provide ``const`` read-only access to inner data - - **GPU compatibility:** LvArray views work seamlessly on device - -.. dropdown:: Example: Views for arrays - :icon: code - - .. code-block:: c++ - - // BAD - Creates a copy - array1d< real64 > copy = originalArray; - - // GOOD - Creates a view (lightweight) - arrayView1d< real64 > view = originalArray.toView(); - - // GOOD - Const view for read-only access - arrayView1d< real64 const > constView = originalArray.toViewConst(); - -View Lifetime Management ------------------------- - -**Never Outlive Parent Arrays** -Dangling views cause segmentation faults and undefined behavior, that can be particularly hard to diagnose. -The rule is applicable to ``arrayView*`` and ``string_view``. - -.. dropdown:: Example: Lifetime Management - :icon: code - - .. code-block:: cpp - - // WRONG - returning view of local array - arrayView1d< real64 > getDanglingView() - { - array1d< real64 > tempArray( 100 ); - return tempArray.toView(); // DANGER: tempArray destroyed, view dangles! - } - - // WRONG - storing view of temporary - class DataHolder - { - arrayView1d< real64 > m_view; - - void setData() - { - array1d< real64 > temp( 100 ); - m_view = temp.toView(); // DANGER: temp destroyed at end of scope! - } - }; - - // CORRECT - return the array - array1d< real64 > getSafeArray() - { - array1d< real64 > result( 100 ); - // populate... - return result; // Move semantics ensure safety - } - - // CORRECT - store array, create view when needed - class SafeDataHolder - { - public: - arrayView1d< real64 > getView() { return m_data.toView(); } - arrayView1d< real64 const > getViewConst() const { return m_data.toViewConst(); } - private: - array1d< real64 > m_data; - }; - -General Architecture -==================== - -Avoid Coupling --------------- - -Minimize coupling between components when possible **without compromising memory efficiency or execution speed**. - -Principles -^^^^^^^^^^ - -- **Loose coupling:** Components should depend on interfaces, not concrete implementations, -- **No circular dependencies:** Consider the existing GEOS dependencies to not make components co-dependent (headers inclusion, packages referencing in ``CMakeLists.txt``, avoid tightly coupled objects), -- **Dependency injection:** Public components should receive their dependencies from external sources. Pass required dependencies using intermediate types instead of direct implementation types, using lambda, templates, and minimal interfaces, relying on **lambda**, **templates** and **minimal interfaces** (loose coupling, testability), -- **Performance exceptions:** Tight coupling is acceptable when required for performance, -- **Minimize header inclusions and dependancies**. - -.. dropdown:: Example: Reducing coupling - :icon: code - - .. code-block:: c++ - - // BAD - Tight coupling to concrete class / implementation - class SolverA - { - void registerSomeMeshData( VTKMesh & mesh ); - void solveHypre( HypreInterface * solver, HypreMatrix * matrix ); - }; - - // GOOD - Depends on interface - class SolverB - { - void registerSomeMeshData( MeshLevel & mesh ); // More general interface - void solve( LAInterface * solver, MatrixBase * matrix ); - }; - - // Performance-critical tight coupling - template< typename FIELD_T > - class Kernel - { - // Direct access to specific data layout for performance - void compute( FIELD_T const & data ); - }; - - template<> Kernel::compute( arrayView1d< double > & field ) - { /* template specialization */ } - - template<> Kernel::compute( arrayView2d< double > & field ); - { /* template specialization */ } - -Avoid Globals and Mutable State --------------------------------- - -**Minimize mutable global states when possible.** -Prefer passing context explicitly. - -.. dropdown:: Why to avoid globals? - :icon: info - - - **Thread safety:** Globals can cause race conditions in parallel code - - **Testability:** Hard to test code that depends on global state - - **Predictability:** Global state makes code behaviour harder to understand - - **Encapsulation:** Violates modularity principles - -.. dropdown:: Example: Avoiding global state - :icon: code - - .. code-block:: c++ - - // BAD - Global mutable state - static real64 g_tolerance = 1.0e-6; - - void solve() - { - if( residual < g_tolerance ) { ... } // Depends on global - } - - // GOOD - Pass as parameter - void solve( real64 tolerance ) - { - if( residual < tolerance ) { ... } - } - - // GOOD - Member variable - class Solver - { - real64 const m_tolerance; - public: - Solver( real64 tolerance ) : m_tolerance( 1.0e-6 ) {} - - void solve() - { - if( residual < m_tolerance ) { ... } - } - }; - -Acceptable Global Usage: - -- **Library wrappers** (MPI wrapper, system-level resources, dependancies interface) -- **Read-only configuration** (immutable after initialization) -- **Performance counters** (for profiling) - -Magic values, ``Group`` paths -------------------------------- - -**Avoid magic values**: - -- **arbitrary values should not be written more than once**, define constants, consider using or extending ``PhysicsConstants.hpp`` / ``Units.hpp``, -- Because of architecture considerations, **Avoid using full raw ``Group`` path**, especially when the path has parts unrelated to the component consulting it -> prefer as relative paths as possible, -- ``Group`` / ``Wrapper`` name in the data repository is considered as a magic value, **use the ``getCatalogName()`` / ``getName()`` getters**, -- **Prefer to let appear the calculus of constants** rather than writing its value directly without explaination (constexpr has no runtime cost). - -====================== -Physics-Specific Rules -====================== - -Unit Consistency ----------------- - -**Verify physical units**, and **document any non-S.I. units** (variable name, comments). -**Non-S.I. units are only for internal computation use.** - -Physical Bounds Checking ------------------------- - -Unphysical values indicate errors and can cause solver failures. -**Validate that state variables remain within physically meaningful bounds**. - -If a value is not strictly disallowed but does not seem possible for the model, **show a warning to the user that he can disable**. +All contributors should familiarize themselves with these rules to maintain code quality, consistency, and reliability across the GEOS codebase. -.. dropdown:: Example: Avoiding global state - :icon: code +.. toctree:: + :maxdepth: 2 - .. code-block:: cpp + CodingPractices.rst - class ConstitutiveModel - { - void updateState( real64 const pressure, - real64 const temperature, - real64 const porosity ) - { - // Check physical bounds - GEOS_ERROR_IF( pressure < 1e-5, - GEOS_FMT( "Pressure {:.2e} Pa below cavitation limit", pressure ) ); - - GEOS_ERROR_IF( temperature < 0.0, - GEOS_FMT( "Absolute temperature {:.2f} K is negative", temperature ) ); - - GEOS_ERROR_IF( porosity < 0.0 || porosity > 1.0, - GEOS_FMT( "Porosity {:.3f} outside [0,1]", porosity ) ); - - // Warn about unusual but valid values, allow user to disable warning - GEOS_WARNING_IF( isHighTemperatureWarningEnabled && temperature > 1000.0, - GEOS_FMT( "High temperature {:.0f} K may be outside model validity", - temperature ) ); - } - }; + CodeGeosFeatures.rst diff --git a/src/docs/sphinx/developerGuide/Contributing/CodingPractices.rst b/src/docs/sphinx/developerGuide/Contributing/CodingPractices.rst new file mode 100644 index 00000000000..ae2c7682c17 --- /dev/null +++ b/src/docs/sphinx/developerGuide/Contributing/CodingPractices.rst @@ -0,0 +1,415 @@ +.. _CodingPractices: + + +############################################################################### +Coding Practice Rules +############################################################################### + +This section covers general C++ coding practices in place for the GEOS project. + +All contributors should familiarize themselves with these rules to maintain code quality, consistency, and reliability across the GEOS codebase. + +.. note:: + This document uses collapsible code blocks. Click on "Show/Hide Code" to expand examples. + +Memory Managment +================ + +Allocations +^^^^^^^^^^^^^^^^^^^^ + +**Minimize dynamic memory allocation** as much as possible, particularly in performance-critical code, + +- For frequently allocated/deallocated objects, **consider memory pools**, +- For containers with known size, **reserve capacity upfront**. + +As they are intended to low-level allocator codes, **avoid ``new``/``delete`` and ``malloc``/``free`` as much as possible**: + +- When applicable, **Prefer composition** over pointers, +- **Use RAII principles** as much as possible, +- **Prefer ``std::unique_ptr``** over ``std::shared_ptr``. + +Pointer / Reference Function Arguments +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +**Pass object by reference rather than passing pointer types.** + +When passing nullable objects as parameter, pointers can be passed. + +**Pointer parameters must always be null-checked when dereferenced**. +If applying this rule produces repetitive code, consider using a const reference. + +This rule imply the following: +- **A reference may never be null**, +- In a given method, **``this`` may never be null**. + +.. dropdown:: Why? + :icon: info + + Passing a pointer or reference parameter show the intention of modifying the underlying object. + The difference of the two practices is that passing a pointer means that we consider that the + underlying object can be null. + +Provide Views to Arrays +^^^^^^^^^^^^^^^^^^^^^^^^ + +- When possible, prefer provide views to arrays (to const data if possible). +- views must be **passed / captured by value** for function / lambdas. + +The rule is generalizable to ``string_view`` for strings, but not applicable in GPU context. + +.. dropdown:: Why Use Views? + :icon: info + + - **No memory allocation:** Views are lightweight references + - **Mutability correctness:** Depending on the view type, can provide ``const`` read-only access to inner data + - **GPU compatibility:** LvArray views work seamlessly on device + +.. dropdown:: Example: Views for arrays + :icon: code + + .. code-block:: c++ + + // BAD - Creates a copy + array1d< real64 > copy = originalArray; + + // GOOD - Creates a view (lightweight) + arrayView1d< real64 > view = originalArray.toView(); + + // GOOD - Const view for read-only access + arrayView1d< real64 const > constView = originalArray.toViewConst(); + +View Lifetime Management +------------------------ + +**Never Outlive Parent Arrays** +Dangling views cause segmentation faults and undefined behavior, that can be particularly hard to diagnose. +The rule is applicable to ``arrayView*`` and ``string_view``. + +.. dropdown:: Example: Lifetime Management + :icon: code + + .. code-block:: cpp + + // WRONG - returning view of local array + arrayView1d< real64 > getDanglingView() + { + array1d< real64 > tempArray( 100 ); + return tempArray.toView(); // DANGER: tempArray destroyed, view dangles! + } + + // WRONG - storing view of temporary + class DataHolder + { + arrayView1d< real64 > m_view; + + void setData() + { + array1d< real64 > temp( 100 ); + m_view = temp.toView(); // DANGER: temp destroyed at end of scope! + } + }; + + // CORRECT - return the array + array1d< real64 > getSafeArray() + { + array1d< real64 > result( 100 ); + // populate... + return result; // Move semantics ensure safety + } + + // CORRECT - store array, create view when needed + class SafeDataHolder + { + public: + arrayView1d< real64 > getView() { return m_data.toView(); } + arrayView1d< real64 const > getViewConst() const { return m_data.toViewConst(); } + private: + array1d< real64 > m_data; + }; + +Value / Const Function Arguments +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +**Pass large and non-trivially copyable objects by const reference**: + +- A "large" size can be defined as "more that 16 bytes, ", which is 2 pointers, 2 integer / index, or 2 double, +- An object is non-trivially copyable objects when it needs to perform sub-allocation when being copied, + +.. dropdown:: Value / Const Reference function parameters practices + :icon: code + + .. code-block:: c++ + + // 16 bytes (string_view are equivalent to 2 constant pointers): PASS BY VALUE + static constexpr string_view str0 = "Hello"; + string const str1 = getTimeStr(); // constant string are passed as string_view, beware of lifetime! + void func( string_view param ); + + arrayView2d< int > myDataView; + + // 16 bytes: PASS BY VALUE + stdArray< double, 2 > vec2D; + void func( vec2D param ); + + // 12 to 16 bytes (globalIndex is 8, integer is 4 to 8): PASS BY VALUE + struct SmallStruct + { + globalIndex id; + integer count; + }; + void func( SmallStruct param ); + + // 16 to 20 bytes (globalIndex is 8, integer is 4 to 8): PASS BY REFERENCE + struct BigStruct + { + globalIndex id; + integer countA; + integer countB; + }; + void func( BigStruct const & param ); + + // Does dynamic allocation: PASS BY REFERENCE + map< int, long > myMap; + stdVector< int > myList { 123 }; + void func( map< int, long > const & myMap, + stdVector< int > const & myList ); + + // LvArray types parameter practices depends on the intent + array2d< int > myConstantData; + array2d< int > myMutableData;, + // passing as non-view means that we will potencially resize the arrays: PASS BY REFERENCE + void arrayResizer( array2d< const int > & myConstantData + array2d< int > & myMutableData ); + // passing as view means that we may just affect the data, only if non-const template: PASS BY VALUE + void arrayProcess( arrayView2d< const int > myConstantData + arrayView2d< int > myMutableData ); + + // 16 bytes superficially... but does dynamic allocation: PASS BY REFERENCE + struct DynAllocStruct + { + string const a; + }; + void func( DynAllocStruct const & param ); + + // Any Group subclass is large: PASS BY REFERENCE + void modify( DomainPartition & domain ); + +Avoid undesired mutability +^^^^^^^^^^^^^^^^^^^^^^^^^^ + +**Use ``const`` and ``constexpr`` when possible**, enabling: + +- **Compiler optimization,** enables better code generation by giving constraints to the code, +- **Code safety,** prevents accidental modification for constant contexts, +- **Show code intention,** make code more readable. + +Also, **mark methods const if the method is not designed to modify** the object state. + +Constexpr for Compile-Time Constants +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +**Use ``constexpr`` for values known at compile time**. + +The rule is not absolute: when the impact is not significative, and the code needs is getting really unclear, the rule can be ignored. + +.. dropdown:: Example: Constexpr usage + :icon: code + + .. code-block:: c++ + + // Compile-time constants + // A rule of thumb is that oftenly any time a value is constexpr, it can also be static. + static constexpr localIndex numDimensions = 3; + static constexpr real64 tolerance = 1.0e-10; + + // Constexpr functions (evaluated at compile time when possible) + constexpr localIndex computeArraySize( localIndex n ) + { return n * n + 2 * n + 1; } + +Validation / Precision +====================== + +Use Tolerance-Based Comparisons +------------------------------- + +**Always consider proper tolerance** for floating-point numbers comparisons, taking into account rounding errors, even for extreme values. + +.. dropdown:: Example: Correct float comparison + :icon: code + + .. code-block:: c++ + + // BAD - Direct comparison + if( value == 1.0 ) { ... } + if( a == b ) { ... } + + // GOOD - Absolute tolerance + real64 const absTol = 1.0e-12; + if( LvArray::math::abs(value) < absTol ) { ... } + + // GOOD - Relative tolerance + real64 const relTol = 1.0e-8; + real64 const scale = LvArray::math::max( LvArray::math::abs(a), LvArray::math::abs(b) ); + if( LvArray::math::abs(a - b) < relTol * scale ) { ... } + +Division Safety, NaN/Inf Values +------------------------------- + +- **Always verify that denominators are not zero or near-zero before a division**. +- In General, **we should not make any computation that could result in NaN/Inf values**. + +.. dropdown:: Example: Division Safety + :icon: code + + .. code-block:: cpp + + // WRONG - unprotected division + real64 const normalizedResidual = residual / initialResidual; + real64 const strainRate = velocityGradient / thickness; + + // CORRECT - protected division + real64 computeNormalizedResidual( real64 const residual, + real64 const initialResidual ) + { + if( initialResidual > machinePrecision ) + return residual / initialResidual; + else + return residual; // or return a flag indicating special case + } + + // CORRECT - with error reporting + real64 safeDivide( real64 const numerator, + real64 const denominator, + string const & context ) + { + GEOS_ERROR_IF( LvArray::math::abs( denominator ) < machinePrecision, + GEOS_FMT( "Division by zero in {}: denominator = {:.2e}", + context, denominator ) ); + return numerator / denominator; + } + +Overflow Prevention +------------------- + +Overflow leads to undefined behavior and memory corruption. **Validate that operations won't exceed type / container limits**, especially for index calculations. + +Performance +=========== + +Speed Optimization Rules +----------------------------- + +- **Hoist Loop Invariants**, move computations that don't change during iterations outside the loop. +- When it does not critically affect the code architecture and clarity, **fuse multiple related kernels to reduce memory traffic and launch overhead** (i.e., statistics kernels process all physics field at once). +- **Optimize Memory Access for Cache and Coalescing**. Access memory sequentially and ensure coalesced access, especially on GPUs. +- **Minimize Host-Device Transfers**. Keep data on the appropriate memory space and minimize transfers. + +General Architecture +==================== + +Avoid Coupling +-------------- + +Minimize coupling between components when possible **without compromising memory efficiency or execution speed**. + +Principles: + +- **Loose coupling:** Components should depend on interfaces, not concrete implementations, +- **No circular dependencies:** Consider the existing GEOS dependencies to not make components co-dependent (headers inclusion, packages referencing in ``CMakeLists.txt``, avoid tightly coupled objects), +- **Dependency injection:** Public components should receive their dependencies from external sources. Pass required dependencies using intermediate types instead of direct implementation types, using lambda, templates, and minimal interfaces, relying on **lambda**, **templates** and **minimal interfaces** (loose coupling, testability), +- **Performance exceptions:** Tight coupling is acceptable when required for performance, +- **Minimize header inclusions and dependancies**. + +.. dropdown:: Example: Reducing coupling + :icon: code + + .. code-block:: c++ + + // BAD - Tight coupling to concrete class / implementation + class SolverA + { + void registerSomeMeshData( VTKMesh & mesh ); + void solveHypre( HypreInterface * solver, HypreMatrix * matrix ); + }; + + // GOOD - Depends on interface + class SolverB + { + void registerSomeMeshData( MeshLevel & mesh ); // More general interface + void solve( LAInterface * solver, MatrixBase * matrix ); + }; + + // Performance-critical tight coupling + template< typename FIELD_T > + class Kernel + { + // Direct access to specific data layout for performance + void compute( FIELD_T const & data ); + }; + + template<> Kernel::compute( arrayView1d< double > & field ) + { /* template specialization */ } + + template<> Kernel::compute( arrayView2d< double > & field ); + { /* template specialization */ } + +Avoid Globals and Mutable State +-------------------------------- + +**Minimize mutable global states when possible.** +Prefer passing context explicitly. + +.. dropdown:: Why to avoid globals? + :icon: info + + - **Thread safety:** Globals can cause race conditions in parallel code + - **Testability:** Hard to test code that depends on global state + - **Predictability:** Global state makes code behaviour harder to understand + - **Encapsulation:** Violates modularity principles + +.. dropdown:: Example: Avoiding global state + :icon: code + + .. code-block:: c++ + + // BAD - Global mutable state + static real64 g_tolerance = 1.0e-6; + + void solve() + { + if( residual < g_tolerance ) { ... } // Depends on global + } + + // GOOD - Pass as parameter + void solve( real64 tolerance ) + { + if( residual < tolerance ) { ... } + } + + // GOOD - Member variable + class Solver + { + real64 const m_tolerance; + public: + Solver( real64 tolerance ) : m_tolerance( 1.0e-6 ) {} + + void solve() + { + if( residual < m_tolerance ) { ... } + } + }; + +Acceptable Global Usage: + +- **Library wrappers** (MPI wrapper, system-level resources, dependancies interface) +- **Read-only configuration** (immutable after initialization) +- **Performance counters** (for profiling) + +Magic values +------------ + +**Avoid magic values**: + +- **arbitrary values should not be written more than once**, define constants, consider using or extending ``PhysicsConstants.hpp`` / ``Units.hpp``, +- **Prefer to let appear the calculus of constants** rather than writing its value directly without explaination (constexpr has no runtime cost). diff --git a/src/docs/sphinx/developerGuide/Contributing/index_contributing.rst b/src/docs/sphinx/developerGuide/Contributing/index_contributing.rst index dbc125fa326..e81d3a38a72 100644 --- a/src/docs/sphinx/developerGuide/Contributing/index_contributing.rst +++ b/src/docs/sphinx/developerGuide/Contributing/index_contributing.rst @@ -9,6 +9,8 @@ Contributing CodeStyle.rst + CodeRules.rst + GitWorkflow.rst Sphinx.rst