Xyst test code coverage report
Current view: top level - Inciter - ChoCG.hpp (source / functions) Hit Total Coverage
Commit: d790e211db0f2cf155e72db869cf1a5a372c10f5 Lines: 51 51 100.0 %
Date: 2025-02-17 15:39:21 Functions: 4 4 100.0 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 0 0 -

           Branch data     Line data    Source code
       1                 :            : // *****************************************************************************
       2                 :            : /*!
       3                 :            :   \file      src/Inciter/ChoCG.hpp
       4                 :            :   \copyright 2012-2015 J. Bakosi,
       5                 :            :              2016-2018 Los Alamos National Security, LLC.,
       6                 :            :              2019-2021 Triad National Security, LLC.,
       7                 :            :              2022-2025 J. Bakosi
       8                 :            :              All rights reserved. See the LICENSE file for details.
       9                 :            :   \brief     ChoCG: Projection-based solver for incompressible flow
      10                 :            : */
      11                 :            : // *****************************************************************************
      12                 :            : 
      13                 :            : #pragma once
      14                 :            : 
      15                 :            : #include <vector>
      16                 :            : #include <map>
      17                 :            : 
      18                 :            : #include "Types.hpp"
      19                 :            : #include "Fields.hpp"
      20                 :            : #include "Table.hpp"
      21                 :            : #include "DerivedData.hpp"
      22                 :            : #include "NodeDiagnostics.hpp"
      23                 :            : #include "PUPUtil.hpp"
      24                 :            : #include "ConjugateGradients.hpp"
      25                 :            : 
      26                 :            : #include "NoWarning/chocg.decl.h"
      27                 :            : 
      28                 :            : namespace inciter {
      29                 :            : 
      30                 :            : //! ChoCG Charm++ chare array used to advance PDEs in time with ChoCG
      31                 :            : class ChoCG : public CBase_ChoCG {
      32                 :            : 
      33                 :            :   public:
      34                 :            :     #if defined(__clang__)
      35                 :            :       #pragma clang diagnostic push
      36                 :            :       #pragma clang diagnostic ignored "-Wunused-parameter"
      37                 :            :       #pragma clang diagnostic ignored "-Wdeprecated-declarations"
      38                 :            :     #elif defined(STRICT_GNUC)
      39                 :            :       #pragma GCC diagnostic push
      40                 :            :       #pragma GCC diagnostic ignored "-Wunused-parameter"
      41                 :            :       #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
      42                 :            :     #elif defined(__INTEL_COMPILER)
      43                 :            :       #pragma warning( push )
      44                 :            :       #pragma warning( disable: 1478 )
      45                 :            :     #endif
      46                 :            :     // Include Charm++ SDAG code. See http://charm.cs.illinois.edu/manuals/html/
      47                 :            :     // charm++/manual.html, Sec. "Structured Control Flow: Structured Dagger".
      48                 :            :     ChoCG_SDAG_CODE
      49                 :            :     #if defined(__clang__)
      50                 :            :       #pragma clang diagnostic pop
      51                 :            :     #elif defined(STRICT_GNUC)
      52                 :            :       #pragma GCC diagnostic pop
      53                 :            :     #elif defined(__INTEL_COMPILER)
      54                 :            :       #pragma warning( pop )
      55                 :            :     #endif
      56                 :            : 
      57                 :            :     //! Constructor
      58                 :            :     explicit ChoCG( const CProxy_Discretization& disc,
      59                 :            :                     const tk::CProxy_ConjugateGradients& cgpre,
      60                 :            :                     const tk::CProxy_ConjugateGradients& cgmom,
      61                 :            :                     const std::map< int, std::vector< std::size_t > >& bface,
      62                 :            :                     const std::map< int, std::vector< std::size_t > >& bnode,
      63                 :            :                     const std::vector< std::size_t >& triinpoel );
      64                 :            : 
      65                 :            :     #if defined(__clang__)
      66                 :            :       #pragma clang diagnostic push
      67                 :            :       #pragma clang diagnostic ignored "-Wundefined-func-template"
      68                 :            :     #endif
      69                 :            :     //! Migrate constructor
      70                 :            :     // cppcheck-suppress uninitMemberVar
      71                 :        827 :     explicit ChoCG( CkMigrateMessage* m ) : CBase_ChoCG( m ) {}
      72                 :            :     #if defined(__clang__)
      73                 :            :       #pragma clang diagnostic pop
      74                 :            :     #endif
      75                 :            : 
      76                 :            :     //! Configure Charm++ custom reduction types initiated from this chare array
      77                 :            :     static void registerReducers();
      78                 :            : 
      79                 :            :     //! Return from migration
      80                 :            :     void ResumeFromSync() override;
      81                 :            : 
      82                 :            :     //! Start setup for solution
      83                 :            :     void setup( tk::real v );
      84                 :            : 
      85                 :            :     //! Initialize Poisson solve
      86                 :            :     void pinit();
      87                 :            : 
      88                 :            :     //! Solve Poisson equation
      89                 :            :     void psolve();
      90                 :            : 
      91                 :            :     //! Continue after Poisson solve
      92                 :            :     void psolved();
      93                 :            : 
      94                 :            :     //! Solve momentum/transport equations
      95                 :            :     void msolve();
      96                 :            : 
      97                 :            :     //! Continue after momentum/transport ssolve
      98                 :            :     void msolved();
      99                 :            : 
     100                 :            :     // Start time stepping
     101                 :            :     void start();
     102                 :            : 
     103                 :            :     //! Advance equations to next time step
     104                 :            :     void advance( tk::real newdt );
     105                 :            : 
     106                 :            :     //! Compute righ-hand side vector of transport equations
     107                 :            :     void rhs();
     108                 :            : 
     109                 :            :     //! Evaluate diagnostics
     110                 :            :     void diag();
     111                 :            : 
     112                 :            :     //! Start (re-)computing domain and boundary integrals
     113                 :            :     void feop();
     114                 :            : 
     115                 :            :     //! Receive contributions to boundary point normals on chare-boundaries
     116                 :            :     void comnorm( const std::unordered_map< int,
     117                 :            :       std::unordered_map< std::size_t, std::array<tk::real,4> > >& inbnd );
     118                 :            : 
     119                 :            :     //! Receive contributions to velocity gradients
     120                 :            :     void comvgrad( const std::unordered_map< std::size_t,
     121                 :            :                            std::vector< tk::real > >& ingrad );
     122                 :            : 
     123                 :            :     //! Receive contributions to momentum flux on chare-boundaries
     124                 :            :     void comflux( const std::unordered_map< std::size_t,
     125                 :            :                           std::vector< tk::real > >& influx );
     126                 :            : 
     127                 :            :     //! Receive contributions to conjugate gradients solution gradient
     128                 :            :     void comsgrad( const std::unordered_map< std::size_t,
     129                 :            :                             std::vector< tk::real > >& ingrad );
     130                 :            : 
     131                 :            :     //! Receive contributions to pressure gradient
     132                 :            :     void compgrad( const std::unordered_map< std::size_t,
     133                 :            :                            std::vector< tk::real > >& ingrad );
     134                 :            : 
     135                 :            :     //! Receive contributions to right-hand side vector on chare-boundaries
     136                 :            :     void comrhs( const std::unordered_map< std::size_t,
     137                 :            :                          std::vector< tk::real > >& inrhs );
     138                 :            : 
     139                 :            :     //! Receive contributions to velocity divergence on chare-boundaries
     140                 :            :     void comdiv( const std::unordered_map< std::size_t, tk::real >& indiv );
     141                 :            : 
     142                 :            :     //! Evaluate residuals
     143                 :            :     void evalres( const std::vector< tk::real >& l2res );
     144                 :            : 
     145                 :            :     //! Receive new mesh from Refiner
     146                 :            :     void resizePostAMR(
     147                 :            :       const std::vector< std::size_t >& ginpoel,
     148                 :            :       const tk::UnsMesh::Chunk& chunk,
     149                 :            :       const tk::UnsMesh::Coords& coord,
     150                 :            :       const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& addedNodes,
     151                 :            :       const std::unordered_map< std::size_t, std::size_t >& addedTets,
     152                 :            :       const std::set< std::size_t >& removedNodes,
     153                 :            :       const std::unordered_map< int, std::unordered_set< std::size_t > >&
     154                 :            :         nodeCommMap,
     155                 :            :       const std::map< int, std::vector< std::size_t > >& bface,
     156                 :            :       const std::map< int, std::vector< std::size_t > >& bnode,
     157                 :            :       const std::vector< std::size_t >& triinpoel );
     158                 :            : 
     159                 :            :     //! Const-ref access to current solution
     160                 :            :     //! \return Const-ref to current solution
     161                 :            :     const tk::Fields& solution() const { return m_u; }
     162                 :            : 
     163                 :            :     //! Compute integral quantities for output
     164                 :            :     void integrals();
     165                 :            : 
     166                 :            :     //! Compute recent conjugate gradients solution gradient
     167                 :            :     void sgrad();
     168                 :            : 
     169                 :            :     //! Evaluate whether to continue with next time step
     170                 :            :     void step();
     171                 :            : 
     172                 :            :     // Evaluate whether to do load balancing
     173                 :            :     void evalLB( int nrestart );
     174                 :            : 
     175                 :            :     /** @name Charm++ pack/unpack serializer member functions */
     176                 :            :     ///@{
     177                 :            :     //! \brief Pack/Unpack serialize member function
     178                 :            :     //! \param[in,out] p Charm++'s PUP::er serializer object reference
     179                 :       2599 :     void pup( PUP::er &p ) override {
     180                 :            :       p | m_disc;
     181                 :            :       p | m_cgpre;
     182                 :            :       p | m_cgmom;
     183                 :       2599 :       p | m_nrhs;
     184                 :       2599 :       p | m_nnorm;
     185                 :       2599 :       p | m_nsgrad;
     186                 :       2599 :       p | m_npgrad;
     187                 :       2599 :       p | m_nvgrad;
     188                 :       2599 :       p | m_nflux;
     189                 :       2599 :       p | m_ndiv;
     190                 :       2599 :       p | m_nbpint;
     191                 :       2599 :       p | m_np;
     192                 :       2599 :       p | m_bnode;
     193                 :       2599 :       p | m_bface;
     194                 :       2599 :       p | m_triinpoel;
     195                 :       2599 :       p | m_u;
     196                 :       2599 :       p | m_un;
     197                 :       2599 :       p | m_pr;
     198                 :       2599 :       p | m_rhs;
     199                 :       2599 :       p | m_rhsc;
     200                 :       2599 :       p | m_sgrad;
     201                 :       2599 :       p | m_sgradc;
     202                 :       2599 :       p | m_pgrad;
     203                 :       2599 :       p | m_pgradc;
     204                 :       2599 :       p | m_vgrad;
     205                 :       2599 :       p | m_vgradc;
     206                 :       2599 :       p | m_flux;
     207                 :       2599 :       p | m_fluxc;
     208                 :       2599 :       p | m_div;
     209                 :       2599 :       p | m_divc;
     210                 :            :       p | m_diag;
     211                 :       2599 :       p | m_bnorm;
     212                 :       2599 :       p | m_bnormc;
     213                 :       2599 :       p | m_bndpoinint;
     214                 :       2599 :       p | m_domedgeint;
     215                 :       2599 :       p | m_bpint;
     216                 :            :       p | m_dsupedge;
     217                 :            :       p | m_dsupint;
     218                 :       2599 :       p | m_dirbcmask;
     219                 :       2599 :       p | m_dirbcval;
     220                 :       2599 :       p | m_dirbcmaskp;
     221                 :       2599 :       p | m_dirbcvalp;
     222                 :       2599 :       p | m_symbcnodes;
     223                 :       2599 :       p | m_symbcnorms;
     224                 :       2599 :       p | m_noslipbcnodes;
     225                 :       2599 :       p | m_surfint;
     226                 :       2599 :       p | m_stage;
     227                 :       2599 :       p | m_finished;
     228                 :       2599 :       p | m_rkcoef;
     229                 :       2599 :       p | m_timer;
     230                 :       2599 :     }
     231                 :            :     //! \brief Pack/Unpack serialize operator|
     232                 :            :     //! \param[in,out] p Charm++'s PUP::er serializer object reference
     233                 :            :     //! \param[in,out] i ChoCG object reference
     234                 :            :     friend void operator|( PUP::er& p, ChoCG& i ) { i.pup(p); }
     235                 :            :     //@}
     236                 :            : 
     237                 :            :   private:
     238                 :            :     //! Discretization proxy
     239                 :            :     CProxy_Discretization m_disc;
     240                 :            :     //! Conjugate Gradients Charm++ proxy for pressure solve
     241                 :            :     tk::CProxy_ConjugateGradients m_cgpre;
     242                 :            :     //! Conjugate Gradients Charm++ proxy for momentum solve
     243                 :            :     tk::CProxy_ConjugateGradients m_cgmom;
     244                 :            :     //! Counter for right-hand side vector nodes updated
     245                 :            :     std::size_t m_nrhs;
     246                 :            :     //! Counter for receiving boundary point normals
     247                 :            :     std::size_t m_nnorm;
     248                 :            :     //! Counter for receiving conjugrate gradient solution gradient
     249                 :            :     std::size_t m_nsgrad;
     250                 :            :     //! Counter for receiving pressure gradient
     251                 :            :     std::size_t m_npgrad;
     252                 :            :     //! Counter for receiving velocity gradient
     253                 :            :     std::size_t m_nvgrad;
     254                 :            :     //! Counter for receiving momentum flux
     255                 :            :     std::size_t m_nflux;
     256                 :            :     //! Counter for receiving boundary velocity divergences
     257                 :            :     std::size_t m_ndiv;
     258                 :            :     //! Counter for receiving boundary point integrals
     259                 :            :     std::size_t m_nbpint;
     260                 :            :     //! Count number of Poisson solves during setup
     261                 :            :     std::size_t m_np;
     262                 :            :     //! Boundary node lists mapped to side set ids used in the input file
     263                 :            :     std::map< int, std::vector< std::size_t > > m_bnode;
     264                 :            :     //! Boundary face lists mapped to side set ids used in the input file
     265                 :            :     std::map< int, std::vector< std::size_t > > m_bface;
     266                 :            :     //! Boundary triangle face connecitivity where BCs are set by user
     267                 :            :     std::vector< std::size_t > m_triinpoel;
     268                 :            :     //! Unknown/solution vector at mesh nodes
     269                 :            :     tk::Fields m_u;
     270                 :            :     //! Unknown/solution vector at mesh nodes at previous time step
     271                 :            :     tk::Fields m_un;
     272                 :            :     //! Pressure
     273                 :            :     std::vector< tk::real > m_pr;
     274                 :            :     //! Right-hand side vector (for the high order system)
     275                 :            :     tk::Fields m_rhs;
     276                 :            :     //! Receive buffer for communication of the right hand side
     277                 :            :     //! \details Key: global node id, value: rhs for all scalar components per
     278                 :            :     //!   node.
     279                 :            :     std::unordered_map< std::size_t, std::vector< tk::real > > m_rhsc;
     280                 :            :     //! Conjugate gradient solution gradient in mesh nodes
     281                 :            :     tk::Fields m_sgrad;
     282                 :            :     //! Conjugate gradient solution gradient receive buffer
     283                 :            :     std::unordered_map< std::size_t, std::vector< tk::real > > m_sgradc;
     284                 :            :     //! Pressure gradient in mesh nodes
     285                 :            :     tk::Fields m_pgrad;
     286                 :            :     //! Pressure gradient receive buffer
     287                 :            :     std::unordered_map< std::size_t, std::vector< tk::real > > m_pgradc;
     288                 :            :     //! Velocity gradient in mesh nodes
     289                 :            :     tk::Fields m_vgrad;
     290                 :            :     //! Velocity gradient receive buffer
     291                 :            :     std::unordered_map< std::size_t, std::vector< tk::real > > m_vgradc;
     292                 :            :     //! Momentum flux in mesh nodes
     293                 :            :     tk::Fields m_flux;
     294                 :            :     //! Momentum flux receive buffer
     295                 :            :     std::unordered_map< std::size_t, std::vector< tk::real > > m_fluxc;
     296                 :            :     //! Velocity divergence
     297                 :            :     std::vector< tk::real > m_div;
     298                 :            :     //! Receive buffer for communication of the velocity divergence
     299                 :            :     //! \details Key: global node id, value: velocity divergence
     300                 :            :     std::unordered_map< std::size_t, tk::real > m_divc;
     301                 :            :     //! Diagnostics object
     302                 :            :     NodeDiagnostics m_diag;
     303                 :            :     //! Boundary point normals
     304                 :            :     //! \details Outer key: side set id. Inner key: global node id of boundary
     305                 :            :     //!   point, value: weighted normal vector, inverse distance square.
     306                 :            :     std::unordered_map< int,
     307                 :            :       std::unordered_map< std::size_t, std::array< tk::real, 4 > > > m_bnorm;
     308                 :            :     //! Boundary point normals receive buffer
     309                 :            :     //! \details Outer key: side set id. Inner key: global node id of boundary
     310                 :            :     //!   point, value: weighted normals and inverse distance square.
     311                 :            :     decltype(m_bnorm) m_bnormc;
     312                 :            :     //! Boundary point integrals
     313                 :            :     //! \details Key: global node id of boundary point, value: boundary point
     314                 :            :     //!   integral contributions.
     315                 :            :     std::unordered_map< std::size_t, std::array<tk::real,3> > m_bndpoinint;
     316                 :            :     //! Domain edge integrals
     317                 :            :     std::unordered_map< tk::UnsMesh::Edge, std::array< tk::real, 5 >,
     318                 :            :       tk::UnsMesh::Hash<2>, tk::UnsMesh::Eq<2> > m_domedgeint;
     319                 :            :     //! Streamable boundary point integrals
     320                 :            :     std::vector< tk::real > m_bpint;
     321                 :            :     //! Superedge (tet, face, edge) end points with local ids for domain edges
     322                 :            :     std::array< std::vector< std::size_t >, 3 > m_dsupedge;
     323                 :            :     //! Superedge (tet, face, edge) domain edge integrals
     324                 :            :     std::array< std::vector< tk::real >, 3 > m_dsupint;
     325                 :            :     //! Nodes and their Dirichlet BC masks
     326                 :            :     std::vector< std::size_t > m_dirbcmask;
     327                 :            :     //! Nodes and their Dirichlet BC values
     328                 :            :     std::vector< double > m_dirbcval;
     329                 :            :     //! Nodes and their pressure Dirichlet BC masks
     330                 :            :     std::vector< std::size_t > m_dirbcmaskp;
     331                 :            :     //! Nodes and their pressure Dirichlet BC values
     332                 :            :     std::vector< double > m_dirbcvalp;
     333                 :            :     //! Streamable nodes at which symmetry BCs are set
     334                 :            :     std::vector< std::size_t > m_symbcnodes;
     335                 :            :     //! Streamable normals at nodes at which symmetry BCs are set
     336                 :            :     std::vector< tk::real > m_symbcnorms;
     337                 :            :     //! Streamable nodes at which noslip BCs are set
     338                 :            :     std::vector< std::size_t > m_noslipbcnodes;
     339                 :            :     //! Streamable surface integral nodes and normals * dA on surfaces
     340                 :            :     std::map< int, std::pair< std::vector< std::size_t >,
     341                 :            :                               std::vector< tk::real > > > m_surfint;
     342                 :            :     //! Runge-Kutta stage counter
     343                 :            :     std::size_t m_stage;
     344                 :            :     //! True in the last time step
     345                 :            :     int m_finished;
     346                 :            :     //! Runge-Kutta coefficients
     347                 :            :     std::vector< tk::real > m_rkcoef;
     348                 :            :     //! Timer
     349                 :            :     std::vector< tk::Timer > m_timer;
     350                 :            : 
     351                 :            :     //! Access bound Discretization class pointer
     352                 :     374318 :     Discretization* Disc() const {
     353                 :            :       Assert( m_disc[ thisIndex ].ckLocal() != nullptr, "ckLocal() null" );
     354                 :     748636 :       return m_disc[ thisIndex ].ckLocal();
     355                 :            :     }
     356                 :            : 
     357                 :            :     //! Access bound momentum solve matrix
     358                 :         40 :     tk::CSR& Lhs() {
     359                 :            :       Assert( m_cgmom[ thisIndex ].ckLocal() != nullptr, "ckLocal() null" );
     360                 :         80 :       return m_cgmom[ thisIndex ].ckLocal()->lhs();
     361                 :            :     }
     362                 :            : 
     363                 :            :    //! Prepare Dirichlet boundary condition data structures
     364                 :            :    void setupDirBC( const std::vector< std::vector< int > >& cfgmask,
     365                 :            :                     const std::vector< std::vector< double > >& cfgval,
     366                 :            :                     std::size_t ncomp,
     367                 :            :                     std::vector< std::size_t >& mask,
     368                 :            :                     std::vector< double >& val );
     369                 :            : 
     370                 :            :     //! Start computing velocity divergence
     371                 :            :     void div( const tk::Fields& u );
     372                 :            : 
     373                 :            :     //! Start computing velocity gradient
     374                 :            :     void velgrad();
     375                 :            : 
     376                 :            :     //! Start computing momentum flux
     377                 :            :     void flux();
     378                 :            : 
     379                 :            :     //! Finalize computing gradient
     380                 :            :     void fingrad( tk::Fields& grad,
     381                 :            :       std::unordered_map< std::size_t, std::vector< tk::real > >& gradc );
     382                 :            : 
     383                 :            :     //! Finalize computing pressure gradient
     384                 :            :     void finpgrad();
     385                 :            : 
     386                 :            :     //! Compute pressure gradient
     387                 :            :     void pgrad();
     388                 :            : 
     389                 :            :     //! Compute local contributions to domain edge integrals
     390                 :            :     void domint();
     391                 :            : 
     392                 :            :     //! Setup lhs matrix for pressure solve
     393                 :            :     std::tuple< tk::CSR, std::vector< tk::real >, std::vector< tk::real > >
     394                 :            :     prelhs( const std::pair< std::vector< std::size_t >,
     395                 :            :                              std::vector< std::size_t > >& psup );
     396                 :            : 
     397                 :            :     //! Setup empty lhs matrix for momentum solve
     398                 :            :     std::tuple< tk::CSR, std::vector< tk::real >, std::vector< tk::real > >
     399                 :            :     momlhs( const std::pair< std::vector< std::size_t >,
     400                 :            :                              std::vector< std::size_t > >& psup );
     401                 :            : 
     402                 :            :     //! Compute chare-boundary edges
     403                 :            :     void bndEdges();
     404                 :            : 
     405                 :            :     //! Compute local contributions to boundary normals and integrals
     406                 :            :     void bndint();
     407                 :            : 
     408                 :            :     //! Combine own and communicated portions of the boundary point normals
     409                 :            :     void bnorm();
     410                 :            : 
     411                 :            :     //! Convert integrals into streamable data structures
     412                 :            :     void streamable();
     413                 :            : 
     414                 :            :     //! Generate superedge-groups for domain-edge loops
     415                 :            :     void domsuped();
     416                 :            : 
     417                 :            :     //! Output mesh and particle fields to files
     418                 :            :     void out();
     419                 :            : 
     420                 :            :     //! Output mesh-based fields to file
     421                 :            :     void writeFields( CkCallback cb );
     422                 :            : 
     423                 :            :     //! Combine own and communicated portions of the integrals
     424                 :            :     void merge();
     425                 :            : 
     426                 :            :     //! Fill lhs matrix of transport equations
     427                 :            :     void lhs();
     428                 :            : 
     429                 :            :     //! Advance systems of equations
     430                 :            :     void solve();
     431                 :            : 
     432                 :            :     //! Compute advective-diffusive prediction of momentum/transport
     433                 :            :     void pred();
     434                 :            : 
     435                 :            :     //! Compute pressure correction
     436                 :            :     void corr();
     437                 :            : 
     438                 :            :     //! Optionally refine/derefine mesh
     439                 :            :     void refine();
     440                 :            : 
     441                 :            :     //! Compute time step size
     442                 :            :     void dt();
     443                 :            : 
     444                 :            :     //! Evaluate whether to save checkpoint/restart
     445                 :            :     void evalRestart();
     446                 :            : 
     447                 :            :     //! Apply boundary conditions
     448                 :            :     void BC( tk::Fields& u, tk::real t );
     449                 :            : 
     450                 :            :     //! Apply scalar source to solution
     451                 :            :     void src();
     452                 :            : };
     453                 :            : 
     454                 :            : } // inciter::

Generated by: LCOV version 1.16