Xyst test code coverage report
Current view: top level - Physics - Problems.cpp (source / functions) Hit Total Coverage
Commit: d790e211db0f2cf155e72db869cf1a5a372c10f5 Lines: 449 552 81.3 %
Date: 2025-02-17 15:57:47 Functions: 37 45 82.2 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 172 360 47.8 %

           Branch data     Line data    Source code
       1                 :            : // *****************************************************************************
       2                 :            : /*!
       3                 :            :   \file      src/Physics/Problems.cpp
       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     Problem-specific functions. Initial conditions, source terms.
      10                 :            : */
      11                 :            : // *****************************************************************************
      12                 :            : 
      13                 :            : #include "Problems.hpp"
      14                 :            : #include "EOS.hpp"
      15                 :            : #include "InciterConfig.hpp"
      16                 :            : #include "Box.hpp"
      17                 :            : 
      18                 :            : namespace inciter {
      19                 :            : 
      20                 :            : extern ctr::Config g_cfg;
      21                 :            : 
      22                 :            : } // ::inciter
      23                 :            : 
      24                 :            : namespace problems {
      25                 :            : 
      26                 :            : using inciter::g_cfg;
      27                 :            : 
      28                 :            : namespace userdef {
      29                 :            : 
      30                 :            : static std::vector< tk::real >
      31                 :    1276019 : ic( tk::real, tk::real, tk::real, tk::real )
      32                 :            : // *****************************************************************************
      33                 :            : //! Set homogeneous initial conditions for a generic user-defined problem
      34                 :            : //! \return Values of conserved variables
      35                 :            : // *****************************************************************************
      36                 :            : {
      37                 :            :   // pressure-based solvers
      38                 :            : 
      39                 :    1276019 :   const auto& solver = g_cfg.get< tag::solver >();
      40         [ +  + ]:    1276019 :   if (solver == "chocg") {
      41                 :     659594 :     const auto& ncomp = g_cfg.get< tag::problem_ncomp >();
      42         [ +  - ]:     659594 :     std::vector< tk::real > u( ncomp, 0.0 );
      43         [ +  - ]:     659594 :     auto ic_velocity = g_cfg.get< tag::ic_velocity >();
      44                 :     659594 :     auto large = std::numeric_limits< double >::max();
      45         [ +  - ]:     659594 :     if (std::abs(ic_velocity[0] - large) > 1.0e-12) u[0] = ic_velocity[0];
      46         [ +  - ]:     659594 :     if (std::abs(ic_velocity[1] - large) > 1.0e-12) u[1] = ic_velocity[1];
      47         [ +  - ]:     659594 :     if (std::abs(ic_velocity[2] - large) > 1.0e-12) u[2] = ic_velocity[2];
      48                 :     659594 :     return u;
      49                 :     659594 :   }
      50         [ +  + ]:     616425 :   else if (solver == "lohcg") {
      51                 :     475985 :     const auto& ncomp = g_cfg.get< tag::problem_ncomp >();
      52         [ +  - ]:     475985 :     std::vector< tk::real > u( ncomp, 0.0 );
      53         [ +  - ]:     475985 :     auto ic_velocity = g_cfg.get< tag::ic_velocity >();
      54                 :     475985 :     auto large = std::numeric_limits< double >::max();
      55         [ +  - ]:     475985 :     if (std::abs(ic_velocity[0] - large) > 1.0e-12) u[1] = ic_velocity[0];
      56         [ +  - ]:     475985 :     if (std::abs(ic_velocity[1] - large) > 1.0e-12) u[2] = ic_velocity[1];
      57         [ +  - ]:     475985 :     if (std::abs(ic_velocity[2] - large) > 1.0e-12) u[3] = ic_velocity[2];
      58                 :     475985 :     return u;
      59                 :     475985 :   }
      60                 :            : 
      61                 :            :   // density-based solvers
      62                 :            : 
      63                 :     140440 :   auto ic_density = g_cfg.get< tag::ic_density >();
      64                 :     140440 :   const auto& ic_velocity = g_cfg.get< tag::ic_velocity >();
      65 [ -  + ][ -  - ]:     140440 :   ErrChk( ic_velocity.size() == 3, "ic_velocity must have 3 components" );
         [ -  - ][ -  - ]
      66                 :            : 
      67         [ +  - ]:     140440 :   std::vector< tk::real > u( 5, 0.0 );
      68                 :            : 
      69                 :     140440 :   u[0] = ic_density;
      70                 :     140440 :   u[1] = u[0] * ic_velocity[0];
      71                 :     140440 :   u[2] = u[0] * ic_velocity[1];
      72                 :     140440 :   u[3] = u[0] * ic_velocity[2];
      73                 :            : 
      74                 :     140440 :   auto ic_pressure = g_cfg.get< tag::ic_pressure >();
      75                 :     140440 :   auto ic_energy = g_cfg.get< tag::ic_energy >();
      76                 :     140440 :   auto ic_temperature = g_cfg.get< tag::ic_temperature >();
      77                 :            : 
      78                 :     140440 :   auto largereal = std::numeric_limits< double >::max();
      79                 :            : 
      80         [ +  - ]:     140440 :   if (std::abs(ic_pressure - largereal) > 1.0e-12) {
      81                 :            : 
      82                 :     140440 :     u[4] = eos::totalenergy( u[0], u[1]/u[0], u[2]/u[0], u[3]/u[0],
      83                 :            :                              ic_pressure );
      84                 :            : 
      85         [ -  - ]:          0 :   } else if (std::abs(ic_energy - largereal) > 1.0e-12) {
      86                 :            : 
      87                 :          0 :     u[4] = u[0] * ic_energy;
      88                 :            : 
      89         [ -  - ]:          0 :   } else if (std::abs(ic_temperature - largereal) > 1.0e-12) {
      90                 :            : 
      91                 :          0 :     auto cv = g_cfg.get< tag::mat_spec_heat_const_vol >();
      92         [ -  - ]:          0 :     if (std::abs(cv - largereal) > 1.0e-12) {
      93                 :          0 :       u[4] = u[0] * ic_temperature * cv;
      94                 :            :     }
      95                 :            : 
      96                 :            :   } else {
      97                 :            : 
      98 [ -  - ][ -  - ]:          0 :     Throw( "IC background energy cannot be computed. Must specify "
                 [ -  - ]
      99                 :            :            "one of background pressure, energy, or velocity." );
     100                 :            : 
     101                 :            :   }
     102                 :            : 
     103                 :     140440 :   return u;
     104                 :     140440 : }
     105                 :            : 
     106                 :            : static tk::real
     107                 :       8084 : pic( tk::real, tk::real, tk::real )
     108                 :            : // *****************************************************************************
     109                 :            : //! Set homogeneous initial conditions for a generic user-defined problem
     110                 :            : //! \return Value of pressure
     111                 :            : // *****************************************************************************
     112                 :            : {
     113                 :       8084 :   return 0.0;
     114                 :            : }
     115                 :            : 
     116                 :            : } // userdef::
     117                 :            : 
     118                 :            : namespace nonlinear_energy_growth {
     119                 :            : 
     120                 :            : static std::vector< tk::real >
     121                 :      73908 : ic( tk::real x, tk::real y, tk::real z, tk::real t )
     122                 :            : // *****************************************************************************
     123                 :            : //! Set initial conditions prescribing nonlinear energy growth
     124                 :            : //! \param[in] x X coordinate where to evaluate the solution
     125                 :            : //! \param[in] y Y coordinate where to evaluate the solution
     126                 :            : //! \param[in] z Z coordinate where to evaluate the solution
     127                 :            : //! \param[in] t Time where to evaluate the solution
     128                 :            : //! \return Values of conserved variables
     129                 :            : // *****************************************************************************
     130                 :            : {
     131                 :            :   using std::cos;
     132                 :            : 
     133                 :            :   // manufactured solution parameters
     134                 :      73908 :   auto ce = g_cfg.get< tag::problem_ce >();
     135                 :      73908 :   auto r0 = g_cfg.get< tag::problem_r0 >();
     136                 :      73908 :   auto a = g_cfg.get< tag::problem_alpha >();
     137                 :      73908 :   auto k = g_cfg.get< tag::problem_kappa >();
     138                 :      73908 :   const auto& b = g_cfg.get< tag::problem_beta >();
     139                 :            : 
     140                 :      73908 :   auto ec = [ ce, t ]( tk::real kappa, tk::real h, tk::real p ) {
     141                 :      73908 :     return std::pow( -3.0*(ce + kappa*h*h*t), p );
     142                 :      73908 :   };
     143                 :            : 
     144                 :     295632 :   auto hx = [ x, y, z, b ]() {
     145                 :      73908 :     return cos(b[0]*M_PI*x) * cos(b[1]*M_PI*y) * cos(b[2]*M_PI*z);
     146         [ +  - ]:      73908 :   };
     147                 :            : 
     148                 :            :   // density
     149                 :      73908 :   auto r = r0 + std::exp(-a*t) * (1.0 - x*x - y*y - z*z);
     150                 :            :   // energy
     151                 :      73908 :   auto re = r * ec(k,hx(),-1.0/3.0);
     152                 :            : 
     153         [ +  - ]:     147816 :   return { r, 0.0, 0.0, 0.0, re };
     154                 :      73908 : }
     155                 :            : 
     156                 :            : static std::vector< tk::real >
     157                 :     178308 : src( tk::real x, tk::real y, tk::real z, tk::real t )
     158                 :            : // *****************************************************************************
     159                 :            : //! Compute and return source term for nonlinear energy growth
     160                 :            : //! \param[in] x X coordinate where to evaluate the source
     161                 :            : //! \param[in] y Y coordinate where to evaluate the source
     162                 :            : //! \param[in] z Z coordinate where to evaluate the source
     163                 :            : //! \param[in] t Time where to evaluate the source
     164                 :            : //! \return Source for flow variables + transported scalars
     165                 :            : // *****************************************************************************
     166                 :            : {
     167                 :            :   using std::sin; using std::cos; using std::pow;
     168                 :            : 
     169                 :            :   // manufactured solution parameters
     170                 :     178308 :   auto a = g_cfg.get< tag::problem_alpha >();
     171                 :     178308 :   const auto& b = g_cfg.get< tag::problem_beta >();
     172                 :     178308 :   auto ce = g_cfg.get< tag::problem_ce >();
     173                 :     178308 :   auto kappa = g_cfg.get< tag::problem_kappa >();
     174                 :     178308 :   auto r0 = g_cfg.get< tag::problem_r0 >();
     175                 :            :   // ratio of specific heats
     176                 :     178308 :   auto g = g_cfg.get< tag::mat_spec_heat_ratio >();
     177                 :            :   // spatial component of density field
     178                 :     178308 :   auto gx = 1.0 - x*x - y*y - z*z;
     179                 :            :   // derivative of spatial component of density field
     180                 :     178308 :   std::array< tk::real, 3 > dg{ -2.0*x, -2.0*y, -2.0*z };
     181                 :            :   // spatial component of energy field
     182                 :     178308 :   auto h = cos(b[0]*M_PI*x) * cos(b[1]*M_PI*y) * cos(b[2]*M_PI*z);
     183                 :            :   // derivative of spatial component of energy field
     184                 :            :   std::array< tk::real, 3 >
     185                 :     178308 :     dh{ -b[0]*M_PI*sin(b[0]*M_PI*x)*cos(b[1]*M_PI*y)*cos(b[2]*M_PI*z),
     186                 :     178308 :         -b[1]*M_PI*cos(b[0]*M_PI*x)*sin(b[1]*M_PI*y)*cos(b[2]*M_PI*z),
     187                 :     356616 :         -b[2]*M_PI*cos(b[0]*M_PI*x)*cos(b[1]*M_PI*y)*sin(b[2]*M_PI*z) };
     188                 :            :   // temporal function f and its derivative
     189                 :     178308 :   auto ft = std::exp(-a*t);
     190                 :     178308 :   auto dfdt = -a*ft;
     191                 :            :   // density and its derivatives
     192                 :     178308 :   auto rho = r0 + ft*gx;
     193                 :     178308 :   std::array< tk::real, 3 > drdx{ ft*dg[0], ft*dg[1], ft*dg[2] };
     194                 :     178308 :   auto drdt = gx*dfdt;
     195                 :            :   // internal energy and its derivatives
     196                 :     178308 :   auto ie = pow( -3.0*(ce + kappa*h*h*t), -1.0/3.0 );
     197                 :     178308 :   std::array< tk::real, 3 > dedx{ 2.0 * pow(ie,4.0) * kappa * h * dh[0] * t,
     198                 :     356616 :                                   2.0 * pow(ie,4.0) * kappa * h * dh[1] * t,
     199                 :     178308 :                                   2.0 * pow(ie,4.0) * kappa * h * dh[2] * t };
     200                 :     178308 :   const auto dedt = kappa * h * h * pow(ie,4.0);
     201                 :            : 
     202         [ +  - ]:     178308 :   std::vector< tk::real > s( 5, 0.0 );
     203                 :            :   // density source
     204                 :     178308 :   s[0] = drdt;
     205                 :            :   // momentum source
     206                 :     178308 :   s[1] = (g-1.0)*(rho*dedx[0] + ie*drdx[0]);
     207                 :     178308 :   s[2] = (g-1.0)*(rho*dedx[1] + ie*drdx[1]);
     208                 :     178308 :   s[3] = (g-1.0)*(rho*dedx[2] + ie*drdx[2]);
     209                 :            :   // energy source
     210                 :     178308 :   s[4] = rho*dedt + ie*drdt;
     211                 :            : 
     212                 :     356616 :   return s;
     213                 :            : }
     214                 :            : 
     215                 :            : } // nonlinear_energy_growth::
     216                 :            : 
     217                 :            : namespace rayleigh_taylor {
     218                 :            : 
     219                 :            : static std::vector< tk::real >
     220                 :    1058942 : ic( tk::real x, tk::real y, tk::real z, tk::real t )
     221                 :            : // *****************************************************************************
     222                 :            : //! Set initial conditions prescribing a Rayleigh-Taylor flow
     223                 :            : //! \param[in] x X coordinate where to evaluate the solution
     224                 :            : //! \param[in] y Y coordinate where to evaluate the solution
     225                 :            : //! \param[in] z Z coordinate where to evaluate the solution
     226                 :            : //! \param[in] t Time where to evaluate the solution
     227                 :            : //! \return Values of conserved variables
     228                 :            : // *****************************************************************************
     229                 :            : {
     230                 :            :   using std::sin; using std::cos;
     231                 :            : 
     232                 :            :   // manufactured solution parameters
     233                 :    1058942 :   auto a = g_cfg.get< tag::problem_alpha >();
     234                 :    1058942 :   const auto& b = g_cfg.get< tag::problem_beta >();
     235                 :    1058942 :   auto p0 = g_cfg.get< tag::problem_p0 >();
     236                 :    1058942 :   auto r0 = g_cfg.get< tag::problem_r0 >();
     237                 :    1058942 :   auto k = g_cfg.get< tag::problem_kappa >();
     238                 :            : 
     239                 :            :   // spatial component of density and pressure fields
     240                 :    1058942 :   tk::real gx = b[0]*x*x + b[1]*y*y + b[2]*z*z;
     241                 :            :   // density
     242                 :    1058942 :   tk::real r = r0 - gx;
     243                 :            :   // velocity
     244                 :    1058942 :   tk::real ft = cos(k*M_PI*t);
     245                 :    1058942 :   tk::real u = ft * z * sin(M_PI*x);
     246                 :    1058942 :   tk::real v = ft * z * cos(M_PI*y);
     247                 :    1058942 :   tk::real w = ft * ( -0.5*M_PI*z*z*(cos(M_PI*x) - sin(M_PI*y)) );
     248                 :            :   // total specific energy
     249                 :    1058942 :   tk::real rE = eos::totalenergy( r, u, v, w, p0 + a*gx );
     250                 :            : 
     251         [ +  - ]:    1058942 :   return { r, r*u, r*v, r*w, rE };
     252                 :            : }
     253                 :            : 
     254                 :            : static std::vector< tk::real >
     255                 :     768460 : src( tk::real x, tk::real y, tk::real z, tk::real t )
     256                 :            : // *****************************************************************************
     257                 :            : //! Compute and return source term for a Rayleigh-Taylor flow
     258                 :            : //! \param[in] x X coordinate where to evaluate the source
     259                 :            : //! \param[in] y Y coordinate where to evaluate the source
     260                 :            : //! \param[in] z Z coordinate where to evaluate the source
     261                 :            : //! \param[in] t Time where to evaluate the source
     262                 :            : //! \return Source for flow variables + transported scalars
     263                 :            : // *****************************************************************************
     264                 :            : {
     265                 :            :   using std::sin; using std::cos;
     266                 :            : 
     267                 :            :   // manufactured solution parameters
     268                 :     768460 :   auto a = g_cfg.get< tag::problem_alpha >();
     269                 :     768460 :   const auto& b = g_cfg.get< tag::problem_beta >();
     270                 :     768460 :   auto k = g_cfg.get< tag::problem_kappa >();
     271                 :     768460 :   auto p0 = g_cfg.get< tag::problem_p0 >();
     272                 :     768460 :   auto g = g_cfg.get< tag::mat_spec_heat_ratio >();
     273                 :            : 
     274                 :            :   // evaluate solution at x,y,z,t
     275         [ +  - ]:     768460 :   auto U = ic( x, y, z, t );
     276                 :            : 
     277                 :            :   // density, velocity, energy, pressure
     278                 :     768460 :   auto rho = U[0];
     279                 :     768460 :   auto u = U[1]/U[0];
     280                 :     768460 :   auto v = U[2]/U[0];
     281                 :     768460 :   auto w = U[3]/U[0];
     282                 :     768460 :   auto E = U[4]/U[0];
     283                 :     768460 :   auto p = p0 + a*(b[0]*x*x + b[1]*y*y + b[2]*z*z);
     284                 :            : 
     285                 :            :   // spatial gradients
     286                 :     768460 :   std::array< tk::real, 3 > drdx{{ -2.0*b[0]*x, -2.0*b[1]*y, -2.0*b[2]*z }};
     287                 :     768460 :   std::array< tk::real, 3 > dpdx{{ 2.0*a*b[0]*x, 2.0*a*b[1]*y, 2.0*a*b[2]*z }};
     288                 :     768460 :   tk::real ft = cos(k*M_PI*t);
     289                 :     768460 :   std::array< tk::real, 3 > dudx{{ ft*M_PI*z*cos(M_PI*x),
     290                 :            :                                    0.0,
     291                 :     768460 :                                    ft*sin(M_PI*x) }};
     292                 :     768460 :   std::array< tk::real, 3 > dvdx{{ 0.0,
     293                 :     768460 :                                    -ft*M_PI*z*sin(M_PI*y),
     294                 :     768460 :                                    ft*cos(M_PI*y) }};
     295                 :     768460 :   std::array< tk::real, 3 > dwdx{{ ft*M_PI*0.5*M_PI*z*z*sin(M_PI*x),
     296                 :     768460 :                                    ft*M_PI*0.5*M_PI*z*z*cos(M_PI*y),
     297                 :     768460 :                                   -ft*M_PI*z*(cos(M_PI*x) - sin(M_PI*y)) }};
     298                 :            :   std::array< tk::real, 3 > dedx{{
     299                 :     768460 :     dpdx[0]/rho/(g-1.0) - p/(g-1.0)/rho/rho*drdx[0]
     300                 :     768460 :     + u*dudx[0] + v*dvdx[0] + w*dwdx[0],
     301                 :     768460 :     dpdx[1]/rho/(g-1.0) - p/(g-1.0)/rho/rho*drdx[1]
     302                 :     768460 :     + u*dudx[1] + v*dvdx[1] + w*dwdx[1],
     303                 :     768460 :     dpdx[2]/rho/(g-1.0) - p/(g-1.0)/rho/rho*drdx[2]
     304                 :    2305380 :     + u*dudx[2] + v*dvdx[2] + w*dwdx[2] }};
     305                 :            : 
     306                 :            :   // time derivatives
     307                 :     768460 :   auto dudt = -k*M_PI*sin(k*M_PI*t)*z*sin(M_PI*x);
     308                 :     768460 :   auto dvdt = -k*M_PI*sin(k*M_PI*t)*z*cos(M_PI*y);
     309                 :     768460 :   auto dwdt =  k*M_PI*sin(k*M_PI*t)/2*M_PI*z*z*(cos(M_PI*x) - sin(M_PI*y));
     310                 :     768460 :   auto dedt = u*dudt + v*dvdt + w*dwdt;
     311                 :            : 
     312         [ +  - ]:     768460 :   std::vector< tk::real > s( 5, 0.0 );
     313                 :            :   // density source
     314                 :     768460 :   s[0] = u*drdx[0] + v*drdx[1] + w*drdx[2];
     315                 :            :   // momentum source
     316                 :     768460 :   s[1] = rho*dudt+u*s[0]+dpdx[0] + U[1]*dudx[0]+U[2]*dudx[1]+U[3]*dudx[2];
     317                 :     768460 :   s[2] = rho*dvdt+v*s[0]+dpdx[1] + U[1]*dvdx[0]+U[2]*dvdx[1]+U[3]*dvdx[2];
     318                 :     768460 :   s[3] = rho*dwdt+w*s[0]+dpdx[2] + U[1]*dwdx[0]+U[2]*dwdx[1]+U[3]*dwdx[2];
     319                 :            :   // energy source
     320                 :     768460 :   s[4] = rho*dedt + E*s[0] + U[1]*dedx[0]+U[2]*dedx[1]+U[3]*dedx[2]
     321                 :     768460 :        + u*dpdx[0]+v*dpdx[1]+w*dpdx[2];
     322                 :            : 
     323                 :    1536920 :   return s;
     324                 :     768460 : }
     325                 :            : 
     326                 :            : } // rayleigh_taylor::
     327                 :            : 
     328                 :            : namespace sedov {
     329                 :            : static std::vector< tk::real >
     330                 :      66650 : ic( tk::real x, tk::real y, tk::real z, tk::real )
     331                 :            : // *****************************************************************************
     332                 :            : //! Set initial conditions prescribing the Sedov blast wave
     333                 :            : //! \param[in] x X coordinate where to evaluate the solution
     334                 :            : //! \param[in] y Y coordinate where to evaluate the solution
     335                 :            : //! \param[in] z Z coordinate where to evaluate the solution
     336                 :            : //! \return Values of conserved variables
     337                 :            : // *****************************************************************************
     338                 :            : {
     339                 :            :   using std::abs;
     340                 :            : 
     341                 :            :   // pressure
     342                 :      66650 :   auto eps = std::numeric_limits< tk::real >::epsilon();
     343                 :            :   tk::real p;
     344 [ +  + ][ +  + ]:      66650 :   if (abs(x) < eps && abs(y) < eps && abs(z) < eps) {
         [ +  + ][ +  + ]
     345                 :          5 :     p = g_cfg.get< tag::problem_p0 >();
     346                 :            :   } else {
     347                 :      66645 :     p = 0.67e-4;
     348                 :            :   }
     349                 :            : 
     350                 :            :   // density
     351                 :      66650 :   tk::real r = 1.0;
     352                 :            :   // velocity
     353                 :      66650 :   tk::real u = 0.0;
     354                 :      66650 :   tk::real v = 0.0;
     355                 :      66650 :   tk::real w = 0.0;
     356                 :            :   // total specific energy
     357                 :      66650 :   tk::real rE = eos::totalenergy( r, u, v, w, p );
     358                 :            : 
     359         [ +  - ]:      66650 :   return { r, r*u, r*v, r*w, rE };
     360                 :            : 
     361                 :            : }
     362                 :            : } // sedov::
     363                 :            : 
     364                 :            : namespace sod {
     365                 :            : static std::vector< tk::real >
     366                 :      13179 : ic( tk::real x, tk::real, tk::real, tk::real )
     367                 :            : // *****************************************************************************
     368                 :            : //! Set initial conditions prescribing the Sod shocktube
     369                 :            : //! \param[in] x X coordinate where to evaluate the solution
     370                 :            : //! \return Values of conserved variables
     371                 :            : // *****************************************************************************
     372                 :            : {
     373                 :            :   tk::real r, p, u, v, w, rE;
     374                 :            : 
     375         [ +  + ]:      13179 :   if (x<0.5) {
     376                 :            :     // density
     377                 :       6662 :     r = 1.0;
     378                 :            :     // pressure
     379                 :       6662 :     p = 1.0;
     380                 :            :   }
     381                 :            :   else {
     382                 :            :     // density
     383                 :       6517 :     r = 0.125;
     384                 :            :     // pressure
     385                 :       6517 :     p = 0.1;
     386                 :            :   }
     387                 :            : 
     388                 :            :   // velocity
     389                 :      13179 :   u = 0.0;
     390                 :      13179 :   v = 0.0;
     391                 :      13179 :   w = 0.0;
     392                 :            : 
     393                 :            :   // total specific energy
     394                 :      13179 :   rE = eos::totalenergy( r, u, v, w, p );
     395                 :            : 
     396         [ +  - ]:      13179 :   return { r, r*u, r*v, r*w, rE };
     397                 :            : }
     398                 :            : } // sod::
     399                 :            : 
     400                 :            : namespace taylor_green {
     401                 :            : 
     402                 :            : static std::vector< tk::real >
     403                 :     280356 : ic( tk::real x, tk::real y, tk::real, tk::real )
     404                 :            : // *****************************************************************************
     405                 :            : //! Set initial conditions prescribing the Taylor-Green vortex
     406                 :            : //! \param[in] x X coordinate where to evaluate the solution
     407                 :            : //! \param[in] y Y coordinate where to evaluate the solution
     408                 :            : //! \return Values of conserved variables
     409                 :            : // *****************************************************************************
     410                 :            : {
     411                 :            :   // density
     412                 :     280356 :   tk::real r = 1.0;
     413                 :            :   // pressure
     414                 :     280356 :   tk::real p = 10.0 + r/4.0*(cos(2.0*M_PI*x) + cos(2.0*M_PI*y));
     415                 :            :   // velocity
     416                 :     280356 :   tk::real u =  sin(M_PI*x) * cos(M_PI*y);
     417                 :     280356 :   tk::real v = -cos(M_PI*x) * sin(M_PI*y);
     418                 :     280356 :   tk::real w = 0.0;
     419                 :            :   // total specific energy
     420                 :     280356 :   auto rE = eos::totalenergy( r, u, v, w, p );
     421                 :            : 
     422         [ +  - ]:     280356 :   return { r, r*u, r*v, r*w, rE };
     423                 :            : }
     424                 :            : 
     425                 :            : static std::vector< tk::real >
     426                 :     932688 : src( tk::real x, tk::real y, tk::real, tk::real )
     427                 :            : // *****************************************************************************
     428                 :            : //! Compute and return source term for a the Taylor-Green vortex
     429                 :            : //! \param[in] x X coordinate where to evaluate the source
     430                 :            : //! \param[in] y Y coordinate where to evaluate the source
     431                 :            : //! \return Source for flow variables + transported scalars
     432                 :            : // *****************************************************************************
     433                 :            : {
     434                 :            :   using std::cos;
     435                 :            : 
     436         [ +  - ]:     932688 :   std::vector< tk::real > s( 5, 0.0 );
     437                 :    1865376 :   s[4] = 3.0*M_PI/8.0*( cos(3.0*M_PI*x)*cos(M_PI*y)
     438                 :     932688 :                       - cos(3.0*M_PI*y)*cos(M_PI*x) );
     439                 :            : 
     440                 :     932688 :   return s;
     441                 :            : }
     442                 :            : 
     443                 :            : } // taylor_green::
     444                 :            : 
     445                 :            : namespace vortical_flow {
     446                 :            : 
     447                 :            : static std::vector< tk::real >
     448                 :    2350800 : ic( tk::real x, tk::real y, tk::real z, tk::real )
     449                 :            : // *****************************************************************************
     450                 :            : //! Set initial conditions prescribing vortical flow
     451                 :            : //! \param[in] x X coordinate where to evaluate the solution
     452                 :            : //! \param[in] y Y coordinate where to evaluate the solution
     453                 :            : //! \param[in] z Z coordinate where to evaluate the solution
     454                 :            : //! \return Values of conserved variables
     455                 :            : // *****************************************************************************
     456                 :            : {
     457                 :            :   // manufactured solution parameters
     458                 :    2350800 :   tk::real a = g_cfg.get< tag::problem_alpha >();
     459                 :    2350800 :   tk::real k = g_cfg.get< tag::problem_kappa >();
     460                 :    2350800 :   tk::real p0 = g_cfg.get< tag::problem_p0 >();
     461                 :            :   // ratio of specific heats
     462                 :    2350800 :   auto g = g_cfg.get< tag::mat_spec_heat_ratio >();
     463                 :            :   // velocity
     464                 :    2350800 :   tk::real ru = a*x - k*y;
     465                 :    2350800 :   tk::real rv = k*x + a*y;
     466                 :    2350800 :   tk::real rw = -2.0*a*z;
     467                 :            :   // total specific energy
     468                 :    2350800 :   tk::real rE = (ru*ru + rv*rv + rw*rw)/2.0 + (p0 - 2.0*a*a*z*z) / (g - 1.0);
     469                 :            : 
     470         [ +  - ]:    2350800 :   return { 1.0, ru, rv, rw, rE };
     471                 :            : }
     472                 :            : 
     473                 :            : static std::vector< tk::real >
     474                 :    1474909 : src( tk::real x, tk::real y, tk::real z, tk::real )
     475                 :            : // *****************************************************************************
     476                 :            : //! Compute and return source term for vortical flow
     477                 :            : //! \param[in] x X coordinate where to evaluate the source
     478                 :            : //! \param[in] y Y coordinate where to evaluate the source
     479                 :            : //! \param[in] z Z coordinate where to evaluate the source
     480                 :            : //! \return Source for flow variables + transported scalars
     481                 :            : // *****************************************************************************
     482                 :            : {
     483                 :            :   // manufactured solution parameters
     484                 :    1474909 :   auto a = g_cfg.get< tag::problem_alpha >();
     485                 :    1474909 :   auto k = g_cfg.get< tag::problem_kappa >();
     486                 :            :   // ratio of specific heats
     487                 :    1474909 :   auto g = g_cfg.get< tag::mat_spec_heat_ratio >();
     488                 :            :   // evaluate solution at x,y,z
     489         [ +  - ]:    1474909 :   auto u = ic( x, y, z, 0.0 );
     490                 :            : 
     491         [ +  - ]:    1474909 :   std::vector< tk::real > s( 5, 0.0 );
     492                 :            :   // momentum source
     493                 :    1474909 :   s[1] = a*u[1]/u[0] - k*u[2]/u[0];
     494                 :    1474909 :   s[2] = k*u[1]/u[0] + a*u[2]/u[0];
     495                 :            :   // energy source
     496                 :    1474909 :   s[4] = (s[1]*u[1] + s[2]*u[2])/u[0] + 8.0*a*a*a*z*z/(g-1.0);
     497                 :            : 
     498                 :    2949818 :   return s;
     499                 :    1474909 : }
     500                 :            : 
     501                 :            : } // vortical_flow::
     502                 :            : 
     503                 :            : namespace slot_cyl {
     504                 :            : 
     505                 :            : static std::vector< tk::real >
     506                 :    6179899 : ic( tk::real x, tk::real y, tk::real, tk::real t )
     507                 :            : // *****************************************************************************
     508                 :            : //! Set initial conditions prescribing slotted cylinder, cone, Gauss hump
     509                 :            : //! \param[in] x X coordinate where to evaluate the solution
     510                 :            : //! \param[in] y Y coordinate where to evaluate the solution
     511                 :            : //! \param[in] t Time where to evaluate the solution
     512                 :            : //! \return Values of conserved variables
     513                 :            : // *****************************************************************************
     514                 :            : {
     515                 :            :   using std::sin; using std::cos; using std::sqrt;
     516                 :            : 
     517                 :            :   // manufactured solution parameters
     518                 :    6179899 :   tk::real p0 = 1.0;
     519                 :            : 
     520                 :            :   // configure number of scalar components
     521                 :    6179899 :   std::size_t ncomp = 6;
     522                 :    6179899 :   const auto& solver = g_cfg.get< tag::solver >();
     523         [ +  + ]:    6179899 :   if (solver == "chocg") {
     524                 :     997120 :     ncomp = 4;
     525                 :            :   }
     526         [ +  + ]:    5182779 :   else if (solver == "lohcg") {
     527                 :     882080 :     ncomp = 5;
     528                 :            :   }
     529                 :            : 
     530         [ +  - ]:    6179899 :   std::vector< tk::real > u( ncomp, 0.0 );
     531                 :            : 
     532                 :            :   // prescribed velocity: rotate in x-y plane
     533                 :    6179899 :   std::size_t sc = 5;
     534         [ +  + ]:    6179899 :   if (solver == "chocg") {
     535                 :     997120 :     u[0] = 0.5 - y;
     536                 :     997120 :     u[1] = x - 0.5;
     537                 :     997120 :     u[2] = 0.0;
     538                 :     997120 :     sc = 3;
     539                 :            :   }
     540         [ +  + ]:    5182779 :   else if (solver == "lohcg") {
     541                 :     882080 :     u[0] = 0.0;
     542                 :     882080 :     u[1] = 0.5 - y;
     543                 :     882080 :     u[2] = x - 0.5;
     544                 :     882080 :     u[3] = 0.0;
     545                 :     882080 :     sc = 4;
     546                 :            :   }
     547                 :            :   else {
     548                 :    4300699 :     u[0] = 1.0;
     549                 :    4300699 :     u[1] = u[0] * (0.5 - y);
     550                 :    4300699 :     u[2] = u[0] * (x - 0.5);
     551                 :    4300699 :     u[3] = 0.0;
     552                 :    4300699 :     u[4] = eos::totalenergy( u[0], u[1]/u[0], u[2]/u[0], u[3]/u[0], p0 );
     553                 :            :   }
     554                 :            : 
     555                 :    6179899 :   const tk::real R0 = 0.15;
     556                 :            : 
     557                 :            :   // center of the cone
     558                 :    6179899 :   tk::real x0 = 0.5;
     559                 :    6179899 :   tk::real y0 = 0.25;
     560                 :    6179899 :   tk::real r = sqrt((x0-0.5)*(x0-0.5) + (y0-0.5)*(y0-0.5));
     561                 :    6179899 :   tk::real kx = 0.5 + r*sin( t );
     562                 :    6179899 :   tk::real ky = 0.5 - r*cos( t );
     563                 :            : 
     564                 :            :   // center of the hump
     565                 :    6179899 :   x0 = 0.25;
     566                 :    6179899 :   y0 = 0.5;
     567                 :    6179899 :   r = sqrt((x0-0.5)*(x0-0.5) + (y0-0.5)*(y0-0.5));
     568                 :    6179899 :   tk::real hx = 0.5 + r*sin( t-M_PI/2.0 ),
     569                 :    6179899 :            hy = 0.5 - r*cos( t-M_PI/2.0 );
     570                 :            : 
     571                 :            :   // center of the slotted cylinder
     572                 :    6179899 :   x0 = 0.5;
     573                 :    6179899 :   y0 = 0.75;
     574                 :    6179899 :   r = sqrt((x0-0.5)*(x0-0.5) + (y0-0.5)*(y0-0.5));
     575                 :    6179899 :   tk::real cx = 0.5 + r*sin( t+M_PI ),
     576                 :    6179899 :            cy = 0.5 - r*cos( t+M_PI );
     577                 :            : 
     578                 :            :   // end points of the cylinder slot
     579                 :    6179899 :   tk::real i1x = 0.525, i1y = cy - r*cos( std::asin(0.025/r) ),
     580                 :    6179899 :            i2x = 0.525, i2y = 0.8,
     581                 :    6179899 :            i3x = 0.475, i3y = 0.8;
     582                 :            : 
     583                 :            :   // rotate end points of cylinder slot
     584                 :    6179899 :   tk::real ri1x = 0.5 + cos(t)*(i1x-0.5) - sin(t)*(i1y-0.5),
     585                 :    6179899 :            ri1y = 0.5 + sin(t)*(i1x-0.5) + cos(t)*(i1y-0.5),
     586                 :    6179899 :            ri2x = 0.5 + cos(t)*(i2x-0.5) - sin(t)*(i2y-0.5),
     587                 :    6179899 :            ri2y = 0.5 + sin(t)*(i2x-0.5) + cos(t)*(i2y-0.5),
     588                 :    6179899 :            ri3x = 0.5 + cos(t)*(i3x-0.5) - sin(t)*(i3y-0.5),
     589                 :    6179899 :            ri3y = 0.5 + sin(t)*(i3x-0.5) + cos(t)*(i3y-0.5);
     590                 :            : 
     591                 :            :   // direction of slot sides
     592                 :    6179899 :   tk::real v1x = ri2x-ri1x, v1y = ri2y-ri1y,
     593                 :    6179899 :            v2x = ri3x-ri2x, v2y = ri3y-ri2y;
     594                 :            : 
     595                 :            :   // lengths of direction of slot sides vectors
     596                 :    6179899 :   tk::real v1 = sqrt(v1x*v1x + v1y*v1y),
     597                 :    6179899 :            v2 = sqrt(v2x*v2x + v2y*v2y);
     598                 :            : 
     599                 :            :   // cone
     600                 :    6179899 :   r = sqrt((x-kx)*(x-kx) + (y-ky)*(y-ky)) / R0;
     601         [ +  + ]:    6179899 :   if (r<1.0) u[sc] = 0.6*(1.0-r);
     602                 :            : 
     603                 :            :   // hump
     604                 :    6179899 :   r = sqrt((x-hx)*(x-hx) + (y-hy)*(y-hy)) / R0;
     605         [ +  + ]:    6179899 :   if (r<1.0) u[sc] = 0.2*(1.0+cos(M_PI*std::min(r,1.0)));
     606                 :            : 
     607                 :            :   // cylinder
     608                 :    6179899 :   r = sqrt((x-cx)*(x-cx) + (y-cy)*(y-cy)) / R0;
     609                 :    6179899 :   const std::array< tk::real, 2 > r1{{ v1x, v1y }},
     610                 :    6179899 :                                   r2{{ x-ri1x, y-ri1y }};
     611                 :    6179899 :   const auto d1 = (r1[0]*r2[1] - r2[0]*r1[1]) / v1;
     612                 :    6179899 :   const std::array< tk::real, 2 > r3{{ v2x, v2y }},
     613                 :    6179899 :                                   r4{{ x-ri2x, y-ri2y }};
     614                 :    6179899 :   const auto d2 = (r3[0]*r4[1] - r4[0]*r3[1]) / v2;
     615 [ +  + ][ +  + ]:    6179899 :   if (r<1.0 && (d1>0.05 || d1<0.0 || d2<0.0)) u[sc] = 0.6;
         [ +  + ][ +  + ]
     616                 :            : 
     617                 :   12359798 :   return u;
     618                 :            : }
     619                 :            : 
     620                 :            : static std::vector< tk::real >
     621                 :    3681250 : src( tk::real x, tk::real y, tk::real z, tk::real t )
     622                 :            : // *****************************************************************************
     623                 :            : //! Compute and return source term for slotted cylinder, cone, Gauss hump
     624                 :            : //! \param[in] x X coordinate where to evaluate the source
     625                 :            : //! \param[in] y Y coordinate where to evaluate the source
     626                 :            : //! \param[in] z Z coordinate where to evaluate the source
     627                 :            : //! \param[in] t Time where to evaluate the source
     628                 :            : //! \return Source for flow variables + transported scalars
     629                 :            : // *****************************************************************************
     630                 :            : {
     631                 :            :   // evaluate solution at x,y,z,t
     632         [ +  - ]:    3681250 :   auto u = ic( x, y, z, t );
     633                 :            : 
     634                 :            :   // configure number of scalar components
     635                 :    3681250 :   std::size_t ncomp = 6;
     636                 :    3681250 :   const auto& solver = g_cfg.get< tag::solver >();
     637         [ +  + ]:    3681250 :   if (solver == "chocg") {
     638                 :     369140 :     ncomp = 4;
     639                 :            :   }
     640         [ +  + ]:    3312110 :   else if (solver == "lohcg") {
     641                 :     369140 :     ncomp = 5;
     642                 :            :   }
     643                 :            : 
     644         [ +  - ]:    3681250 :   std::vector< tk::real > s( ncomp, 0.0 );
     645                 :            : 
     646                 :            :   // momentum source
     647         [ +  + ]:    3681250 :   if (solver == "chocg") {
     648                 :     369140 :     s[0] = -u[1];
     649                 :     369140 :     s[1] =  u[0];
     650                 :            :   }
     651         [ +  + ]:    3312110 :   else if (solver == "lohcg") {
     652                 :     369140 :     s[1] = -u[2];
     653                 :     369140 :     s[2] =  u[1];
     654                 :            :   }
     655                 :            :   else {
     656                 :    2942970 :     s[1] = -u[2];
     657                 :    2942970 :     s[2] =  u[1];
     658                 :            :   }
     659                 :            : 
     660                 :    7362500 :   return s;
     661                 :    3681250 : }
     662                 :            : 
     663                 :            : } // slot_cyl::
     664                 :            : 
     665                 :            : namespace sheardiff {
     666                 :            : 
     667                 :            : static std::vector< tk::real >
     668                 :          0 : ic( tk::real x, tk::real y, tk::real, tk::real t )
     669                 :            : // *****************************************************************************
     670                 :            : //! Set initial conditions prescribing shear-diffusion in 2D
     671                 :            : //! \param[in] x X coordinate where to evaluate the solution
     672                 :            : //! \param[in] y Y coordinate where to evaluate the solution
     673                 :            : //! \param[in] t Time where to evaluate the solution
     674                 :            : //! \return Values of conserved variables
     675                 :            : //! \see A. Okubo, M.J. Karweit, Diffusion from a continuous source in a uniform
     676                 :            : //!      shear flow, Limnology and Oceanography, 14, 1969,
     677                 :            : //!      https://doi.org/10.4319/lo.1969.14.4.0514.
     678                 :            : // *****************************************************************************
     679                 :            : {
     680                 :            :   using std::exp;
     681                 :            :   using std::pow;
     682                 :            :   using std::sqrt;
     683                 :            : 
     684                 :            :   // manufactured solution parameters
     685                 :          0 :   auto V0 = g_cfg.get< tag::problem_p0 >();  // translation velocity in x
     686                 :          0 :   auto L = g_cfg.get< tag::problem_alpha >();// shear velocity
     687                 :          0 :   auto t0 = g_cfg.get< tag::t0 >();          // initial time
     688                 :          0 :   auto dif = g_cfg.get< tag::mat_dyn_diffusivity >();  // diffusivity
     689                 :          0 :   auto eps = std::numeric_limits< tk::real >::epsilon();
     690 [ -  - ][ -  - ]:          0 :   if (dif < eps) Throw( "Diffusivity must be positive" );
         [ -  - ][ -  - ]
     691                 :            : 
     692                 :            :   // configure number of scalar components
     693                 :          0 :   std::size_t ncomp = 4;
     694                 :          0 :   const auto& solver = g_cfg.get< tag::solver >();
     695         [ -  - ]:          0 :   if (solver == "chocg") {
     696                 :            :   } else
     697         [ -  - ]:          0 :   if (solver == "lohcg") {
     698                 :          0 :     ncomp = 5;
     699                 :            :   }
     700 [ -  - ][ -  - ]:          0 :   else Throw( "Shear-diff IC not setup for this solver" );
                 [ -  - ]
     701                 :            : 
     702                 :            :   // prescribed velocity
     703         [ -  - ]:          0 :   std::vector< tk::real > u( ncomp, 0.0 );
     704                 :          0 :   std::size_t sc = 3;
     705         [ -  - ]:          0 :   if (solver == "chocg") {
     706                 :          0 :     u[0] = V0 + L*y;
     707                 :          0 :     u[1] = u[2] = 0.0;
     708                 :            :   }
     709         [ -  - ]:          0 :   else if (solver == "lohcg") {
     710                 :          0 :     u[1] = V0 + L*y;
     711                 :          0 :     u[2] = u[3] = 0.0;
     712                 :          0 :     sc = 4;
     713                 :            :   }
     714                 :            : 
     715                 :          0 :   auto M = [=](tk::real T){ return 4.0*M_PI*T*sqrt(1.0 + L*L*T*T/12.0); };
     716                 :            : 
     717                 :          0 :   u[sc] = M(t0) / M(t) *
     718                 :          0 :        exp( -pow( x - V0*t - 0.5*L*y*t, 2.0 ) /
     719                 :          0 :              (4.0*dif*t*(1.0 + L*L*t*t/12.0)) - y*y / (4.0*dif*t) );
     720                 :            : 
     721                 :          0 :   return u;
     722                 :            : }
     723                 :            : 
     724                 :            : static std::vector< tk::real >
     725                 :          0 : src( tk::real, tk::real, tk::real, tk::real )
     726                 :            : // *****************************************************************************
     727                 :            : //! Compute and return source term for shear-diffusion in 2D
     728                 :            : //! \return Source for flow variables + transported scalars
     729                 :            : //! \see A. Okubo, M.J. Karweit, Diffusion from a continuous source in a uniform
     730                 :            : //!      shear flow, Limnology and Oceanography, 14, 1969,
     731                 :            : //!      https://doi.org/10.4319/lo.1969.14.4.0514.
     732                 :            : // *****************************************************************************
     733                 :            : {
     734                 :            :   // configure number of scalar components
     735                 :          0 :   std::size_t ncomp = 4;
     736                 :          0 :   const auto& solver = g_cfg.get< tag::solver >();
     737         [ -  - ]:          0 :   if (solver == "lohcg") {
     738                 :          0 :     ncomp = 5;
     739                 :            :   }
     740                 :            : 
     741                 :            :   // manufactured solution parameters
     742                 :          0 :   auto L = g_cfg.get< tag::problem_alpha >();// shear velocity
     743                 :            : 
     744                 :            :   // source
     745         [ -  - ]:          0 :   std::vector< tk::real > s( ncomp, 0.0 );
     746         [ -  - ]:          0 :   if (solver == "lohcg") {
     747                 :          0 :     s[0] = -L;
     748                 :            :   }
     749                 :            : 
     750                 :          0 :   return s;
     751                 :            : }
     752                 :            : 
     753                 :            : } // sheardiff::
     754                 :            : 
     755                 :            : namespace point_src {
     756                 :            : 
     757                 :            : static std::vector< tk::real >
     758                 :       1858 : ic( tk::real x, tk::real y, tk::real z, tk::real t )
     759                 :            : // *****************************************************************************
     760                 :            : //! Set initial conditions for point source problem
     761                 :            : //! \param[in] x X coordinate where to evaluate initial conditions
     762                 :            : //! \param[in] y Y coordinate where to evaluate initial conditions
     763                 :            : //! \param[in] z Z coordinate where to evaluate initial conditions
     764                 :            : //! \param[in] t Time where to evaluate initial conditions
     765                 :            : //! \return Values of conserved variables
     766                 :            : // *****************************************************************************
     767                 :            : {
     768                 :       1858 :   auto u = userdef::ic( x, y, z, t );
     769         [ +  - ]:       1858 :   u.push_back( 0.0 );
     770                 :       1858 :   return u;
     771                 :          0 : }
     772                 :            : 
     773                 :            : static void
     774                 :        600 : src( const std::array< std::vector< tk::real >, 3 >& coord,
     775                 :            :      tk::real t,
     776                 :            :      tk::Fields& U )
     777                 :            : // *****************************************************************************
     778                 :            : //! Apply point-source directly to numerical solution
     779                 :            : //! \param[in] coord Mesh node coordinates
     780                 :            : //! \param[in] t Physical time
     781                 :            : //! \param[in,out] U Solution vector at recent time step
     782                 :            : //! \note This is different from other source terms, because this directly
     783                 :            : //!   modifies the solution instead of applied as a source term mathematically.
     784                 :            : //!   Hence the function signature is also different.
     785                 :            : // *****************************************************************************
     786                 :            : {
     787         [ -  + ]:        600 :   if (U.nprop() == 5) return;
     788                 :            : 
     789                 :        600 :   const auto& source = g_cfg.get< tag::problem_src >();
     790                 :        600 :   const auto& location = source.get< tag::location >();
     791                 :        600 :   auto radius = source.get< tag::radius >();
     792                 :        600 :   auto release_time = source.get< tag::release_time >();
     793                 :        600 :   auto largereal = std::numeric_limits< double >::max();
     794                 :            : 
     795         [ +  - ]:       1200 :   if (location.size() != 3 ||
     796 [ +  - ][ -  + ]:       1200 :       std::abs(radius - largereal) < 1.0e-12 ||
                 [ -  + ]
     797                 :        600 :       std::abs(release_time - largereal) < 1.0e-12)
     798                 :            :   {
     799                 :          0 :     return;
     800                 :            :   }
     801                 :            : 
     802                 :        600 :   auto sx = location[0];
     803                 :        600 :   auto sy = location[1];
     804                 :        600 :   auto sz = location[2];
     805                 :        600 :   auto sr = radius;
     806                 :        600 :   auto st = release_time;
     807                 :            : 
     808         [ -  + ]:        600 :   if (t < st) return;
     809                 :            : 
     810                 :        600 :   const auto& x = coord[0];
     811                 :        600 :   const auto& y = coord[1];
     812                 :        600 :   const auto& z = coord[2];
     813                 :            : 
     814         [ +  + ]:     279300 :   for (std::size_t i=0; i<U.nunk(); ++i) {
     815                 :     278700 :     auto rx = sx - x[i];
     816                 :     278700 :     auto ry = sy - y[i];
     817                 :     278700 :     auto rz = sz - z[i];
     818         [ +  + ]:     278700 :     if (rx*rx + ry*ry + rz*rz < sr*sr) U(i,5) = 1.0;
     819                 :            :   }
     820                 :            : 
     821                 :        600 :   return;
     822                 :            : }
     823                 :            : 
     824                 :            : } // point_src::
     825                 :            : 
     826                 :            : namespace poisson {
     827                 :            : 
     828                 :            : static std::vector< tk::real >
     829                 :       8568 : ic( tk::real, tk::real, tk::real, tk::real )
     830                 :            : // *****************************************************************************
     831                 :            : //! Set velocity initial conditions for testing a Poisson solve only
     832                 :            : //! \return Values for initial conditions
     833                 :            : // *****************************************************************************
     834                 :            : {
     835         [ +  - ]:       8568 :   return { 0, 0, 0 };
     836                 :            : }
     837                 :            : 
     838                 :            : } // poisson::
     839                 :            : 
     840                 :            : namespace poisson_const {
     841                 :            : 
     842                 :            : static tk::real
     843                 :        990 : pr( tk::real, tk::real, tk::real )
     844                 :            : // *****************************************************************************
     845                 :            : //! Set pressure rhs for testing a Poisson solve
     846                 :            : //! \return Value for pressure rhs
     847                 :            : // *****************************************************************************
     848                 :            : {
     849                 :        990 :   return 6.0;
     850                 :            : }
     851                 :            : 
     852                 :            : static tk::real
     853                 :       3781 : ic( tk::real x, tk::real y, tk::real z )
     854                 :            : // *****************************************************************************
     855                 :            : //! Evaluate pressure boundary condition
     856                 :            : //! \param[in] x X coordinate where to evaluate the BC
     857                 :            : //! \param[in] y Y coordinate where to evaluate the BC
     858                 :            : //! \param[in] z Z coordinate where to evaluate the BC
     859                 :            : //! \return Value for pressure BC
     860                 :            : // *****************************************************************************
     861                 :            : {
     862                 :       3781 :   return x*x + y*y + z*z;
     863                 :            : }
     864                 :            : 
     865                 :            : } // poisson_const::
     866                 :            : 
     867                 :            : namespace poisson_harmonic {
     868                 :            : 
     869                 :            : static tk::real
     870                 :          0 : pr( tk::real, tk::real, tk::real )
     871                 :            : // *****************************************************************************
     872                 :            : //! Set pressure rhs for testing a Laplace solve
     873                 :            : //! \return Value for pressure rhs
     874                 :            : // *****************************************************************************
     875                 :            : {
     876                 :          0 :   return 0.0;
     877                 :            : }
     878                 :            : 
     879                 :            : static tk::real
     880                 :          0 : ic( tk::real x, tk::real y, tk::real z )
     881                 :            : // *****************************************************************************
     882                 :            : //! Evaluate pressure boundary condition
     883                 :            : //! \param[in] x X coordinate where to evaluate the BC
     884                 :            : //! \param[in] y Y coordinate where to evaluate the BC
     885                 :            : //! \param[in] z Z coordinate where to evaluate the BC
     886                 :            : //! \return Value for pressure BC
     887                 :            : // *****************************************************************************
     888                 :            : {
     889                 :          0 :   const auto& b = g_cfg.get< tag::problem_beta >();
     890                 :          0 :   auto x0 = b[0];
     891                 :          0 :   auto y0 = b[1];
     892                 :          0 :   auto z0 = b[2];
     893                 :            : 
     894                 :          0 :   return 1.0 / std::sqrt( (x-x0)*(x-x0) + (y-y0)*(y-y0) + (z-z0)*(z-z0) );
     895                 :            : }
     896                 :            : 
     897                 :            : } // poisson_harmonic::
     898                 :            : 
     899                 :            : namespace poisson_sine {
     900                 :            : 
     901                 :            : static tk::real
     902                 :        278 : pr( tk::real x, tk::real y, tk::real z )
     903                 :            : // *****************************************************************************
     904                 :            : //! Set pressure rhs for testing a Poisson solve
     905                 :            : //! \return Value for pressure rhs
     906                 :            : // *****************************************************************************
     907                 :            : {
     908                 :        278 :   return -M_PI * M_PI * x * y * std::sin( M_PI * z );
     909                 :            : }
     910                 :            : 
     911                 :            : static tk::real
     912                 :       1067 : ic( tk::real x, tk::real y, tk::real z )
     913                 :            : // *****************************************************************************
     914                 :            : //! Evaluate pressure boundary condition
     915                 :            : //! \param[in] x X coordinate where to evaluate the BC
     916                 :            : //! \param[in] y Y coordinate where to evaluate the BC
     917                 :            : //! \param[in] z Z coordinate where to evaluate the BC
     918                 :            : //! \return Value for pressure BC
     919                 :            : // *****************************************************************************
     920                 :            : {
     921                 :       1067 :   return x * y * std::sin( M_PI * z );
     922                 :            : }
     923                 :            : 
     924                 :            : } // poisson_sine::
     925                 :            : 
     926                 :            : namespace poisson_sine3 {
     927                 :            : 
     928                 :            : static tk::real
     929                 :        596 : pr( tk::real x, tk::real y, tk::real z )
     930                 :            : // *****************************************************************************
     931                 :            : //! Set pressure rhs for testing a Poisson solve
     932                 :            : //! \return Value for pressure rhs
     933                 :            : // *****************************************************************************
     934                 :            : {
     935                 :            :   using std::sin;
     936                 :            : 
     937                 :        596 :   return -3.0 * M_PI * M_PI * sin(M_PI*x) * sin(M_PI*y) * sin(M_PI*z);
     938                 :            : }
     939                 :            : 
     940                 :            : static tk::real
     941                 :       2281 : ic( tk::real x, tk::real y, tk::real z )
     942                 :            : // *****************************************************************************
     943                 :            : //! Evaluate pressure boundary condition
     944                 :            : //! \param[in] x X coordinate where to evaluate the BC
     945                 :            : //! \param[in] y Y coordinate where to evaluate the BC
     946                 :            : //! \param[in] z Z coordinate where to evaluate the BC
     947                 :            : //! \return Value for pressure BC
     948                 :            : // *****************************************************************************
     949                 :            : {
     950                 :       2281 :   return sin(M_PI*x) * sin(M_PI*y) * sin(M_PI*z);
     951                 :            : }
     952                 :            : 
     953                 :            : } // poisson_sine3::
     954                 :            : 
     955                 :            : namespace poisson_neumann {
     956                 :            : 
     957                 :            : static tk::real
     958                 :        278 : pr( tk::real x, tk::real y, tk::real )
     959                 :            : // *****************************************************************************
     960                 :            : //! Set pressure rhs for testing a Poisson solve
     961                 :            : //! \param[in] x X coordinate where to evaluate the rhs
     962                 :            : //! \param[in] y Y coordinate where to evaluate the rhs
     963                 :            : //! \return Value for pressure rhs
     964                 :            : // *****************************************************************************
     965                 :            : {
     966                 :        278 :   return -3.0 * std::cos(2.0*x) * std::exp(y);
     967                 :            : }
     968                 :            : 
     969                 :            : static std::array< tk::real, 3 >
     970                 :        612 : pg( tk::real x, tk::real y, tk::real )
     971                 :            : // *****************************************************************************
     972                 :            : //! Set pressure gradient for testing a Poisson solve
     973                 :            : //! \param[in] x X coordinate where to evaluate the pressure gradient
     974                 :            : //! \param[in] y Y coordinate where to evaluate the pressure gradient
     975                 :            : //! \return Value for pressure gradient at a point
     976                 :            : // *****************************************************************************
     977                 :            : {
     978                 :        612 :   return { -2.0 * std::sin( 2.0 * x ) * std::exp( y ),
     979                 :        612 :            std::cos(2.0*x) * std::exp(y),
     980                 :        612 :            0.0 };
     981                 :            : }
     982                 :            : 
     983                 :            : static tk::real
     984                 :        972 : ic( tk::real x, tk::real y, tk::real )
     985                 :            : // *****************************************************************************
     986                 :            : //! Evaluate pressure boundary condition
     987                 :            : //! \param[in] x X coordinate where to evaluate the IC / analytic solution
     988                 :            : //! \param[in] y Y coordinate where to evaluate the IC / analytic solution
     989                 :            : //! \return Value for pressure
     990                 :            : // *****************************************************************************
     991                 :            : {
     992                 :        972 :   return std::cos(2.0*x) * std::exp(y);
     993                 :            : }
     994                 :            : 
     995                 :            : } // poisson_neumann::
     996                 :            : 
     997                 :            : 
     998                 :            : namespace poiseuille {
     999                 :            : 
    1000                 :            : static std::vector< tk::real >
    1001                 :          0 : ic( tk::real, tk::real y, tk::real, tk::real )
    1002                 :            : // *****************************************************************************
    1003                 :            : //! Set initial conditions prescribing the Poisuille problem
    1004                 :            : //! \param[in] y Y coordinate where to evaluate the solution
    1005                 :            : //! \return Values of physics variables
    1006                 :            : // *****************************************************************************
    1007                 :            : {
    1008                 :          0 :   auto eps = std::numeric_limits< tk::real >::epsilon();
    1009                 :          0 :   auto nu = g_cfg.get< tag::mat_dyn_viscosity >();
    1010 [ -  - ][ -  - ]:          0 :   if (nu < eps) Throw( "Poiseuille flow needs nonzero viscosity" );
         [ -  - ][ -  - ]
    1011                 :            : 
    1012                 :          0 :   auto dpdx = -0.12;
    1013                 :          0 :   auto u = -dpdx * y * (1.0 - y) / 2.0 / nu;
    1014                 :            : 
    1015                 :          0 :   const auto& solver = g_cfg.get< tag::solver >();
    1016         [ -  - ]:          0 :   if (solver == "chocg") {
    1017         [ -  - ]:          0 :     return { u, 0.0, 0.0 };
    1018                 :            :   } else
    1019         [ -  - ]:          0 :   if (solver == "lohcg") {
    1020         [ -  - ]:          0 :     return { 0.0, 0.0, 0.0, 0.0 };
    1021                 :            :   }
    1022 [ -  - ][ -  - ]:          0 :   else Throw( "Poiseuille IC not setup for this solver" );
                 [ -  - ]
    1023                 :            : }
    1024                 :            : 
    1025                 :            : static std::vector< tk::real >
    1026                 :          0 : sol( tk::real, tk::real y, tk::real, tk::real )
    1027                 :            : // *****************************************************************************
    1028                 :            : //! Set analytic solution of the Poisuille problem
    1029                 :            : //! \param[in] y Y coordinate where to evaluate the solution
    1030                 :            : //! \return Values of physics variables
    1031                 :            : // *****************************************************************************
    1032                 :            : {
    1033                 :          0 :   auto eps = std::numeric_limits< tk::real >::epsilon();
    1034                 :          0 :   auto nu = g_cfg.get< tag::mat_dyn_viscosity >();
    1035 [ -  - ][ -  - ]:          0 :   if (nu < eps) Throw( "Poiseuille flow needs nonzero viscosity" );
         [ -  - ][ -  - ]
    1036                 :            : 
    1037                 :          0 :   auto dpdx = -0.12;
    1038                 :          0 :   auto u = -dpdx * y * (1.0 - y) / 2.0 / nu;
    1039                 :            : 
    1040                 :          0 :   const auto& solver = g_cfg.get< tag::solver >();
    1041         [ -  - ]:          0 :   if (solver == "chocg") {
    1042         [ -  - ]:          0 :     return { u, 0.0, 0.0 };
    1043                 :            :   } else
    1044         [ -  - ]:          0 :   if (solver == "lohcg") {
    1045         [ -  - ]:          0 :     return { 0.0, u, 0.0, 0.0 };
    1046                 :            :   }
    1047 [ -  - ][ -  - ]:          0 :   else Throw( "Poiseuille IC not setup for this solver" );
                 [ -  - ]
    1048                 :            : }
    1049                 :            : 
    1050                 :            : } // poiseuille::
    1051                 :            : 
    1052                 :            : std::function< std::vector< tk::real >
    1053                 :            :              ( tk::real, tk::real, tk::real, tk::real ) >
    1054                 :      65648 : IC()
    1055                 :            : // *****************************************************************************
    1056                 :            : //  Query user config and assign function to set initial conditions
    1057                 :            : //! \return The function to call to set initial conditions
    1058                 :            : // *****************************************************************************
    1059                 :            : {
    1060                 :      65648 :   const auto& problem = inciter::g_cfg.get< tag::problem >();
    1061                 :            : 
    1062                 :            :   std::function< std::vector< tk::real >
    1063                 :      65648 :                ( tk::real, tk::real, tk::real, tk::real ) > ic;
    1064                 :            : 
    1065         [ +  + ]:      65648 :   if (problem == "userdef")
    1066                 :       6356 :     ic = userdef::ic;
    1067         [ +  + ]:      59292 :   else if (problem == "nonlinear_energy_growth")
    1068                 :       1254 :     ic = nonlinear_energy_growth::ic;
    1069         [ +  + ]:      58038 :   else if (problem == "rayleigh_taylor")
    1070                 :       3624 :     ic = rayleigh_taylor::ic;
    1071         [ +  + ]:      54414 :   else if (problem == "sedov")
    1072                 :         59 :     ic = sedov::ic;
    1073         [ +  + ]:      54355 :   else if (problem == "sod")
    1074                 :        629 :     ic = sod::ic;
    1075         [ +  + ]:      53726 :   else if (problem == "taylor_green")
    1076                 :       4812 :     ic = taylor_green::ic;
    1077         [ +  + ]:      48914 :   else if (problem == "vortical_flow")
    1078                 :      29151 :     ic = vortical_flow::ic;
    1079         [ +  + ]:      19763 :   else if (problem == "slot_cyl")
    1080                 :      19624 :     ic = slot_cyl::ic;
    1081         [ -  + ]:        139 :   else if (problem == "sheardiff")
    1082                 :          0 :     ic = sheardiff::ic;
    1083         [ +  + ]:        139 :   else if (problem == "point_src")
    1084                 :          4 :     ic = point_src::ic;
    1085         [ +  - ]:        135 :   else if (problem.find("poisson") != std::string::npos)
    1086                 :        135 :     ic = poisson::ic;
    1087         [ -  - ]:          0 :   else if (problem == "poiseuille")
    1088                 :          0 :     ic = poiseuille::ic;
    1089                 :            :   else
    1090 [ -  - ][ -  - ]:          0 :     Throw( "problem type ic not hooked up" );
                 [ -  - ]
    1091                 :            : 
    1092                 :      65648 :   return ic;
    1093                 :          0 : }
    1094                 :            : 
    1095                 :            : std::function< std::vector< tk::real >
    1096                 :            :              ( tk::real, tk::real, tk::real, tk::real ) >
    1097                 :      47235 : SOL()
    1098                 :            : // *****************************************************************************
    1099                 :            : //  Query user config and assign function to query analytic solutions
    1100                 :            : //! \return The function to call to query analytic solutions
    1101                 :            : // *****************************************************************************
    1102                 :            : {
    1103                 :      47235 :   const auto& problem = inciter::g_cfg.get< tag::problem >();
    1104                 :            : 
    1105         [ +  + ]:      71692 :   if (problem == "userdef" or
    1106         [ +  + ]:      46499 :       problem == "sod" or
    1107 [ +  + ][ +  + ]:      93734 :       problem == "sedov" or
                 [ +  + ]
    1108                 :      20028 :       problem == "point_src")
    1109                 :      27263 :     return {};
    1110         [ -  + ]:      19972 :   else if (problem == "poiseuille")
    1111                 :          0 :     return poiseuille::sol;
    1112                 :            :   else
    1113                 :      19972 :     return IC();
    1114                 :            : }
    1115                 :            : 
    1116                 :            : void
    1117                 :       2786 : initialize( const std::array< std::vector< tk::real >, 3 >& coord,
    1118                 :            :             tk::Fields& U,
    1119                 :            :             tk::real t,
    1120                 :            :             const std::vector< std::unordered_set< std::size_t > >& boxnodes )
    1121                 :            : // *****************************************************************************
    1122                 :            : //  Set inital conditions
    1123                 :            : //! \param[in] coord Mesh node coordinates
    1124                 :            : //! \param[in,out] U Array of unknowns
    1125                 :            : //! \param[in] t Physical time
    1126                 :            : //! \param[in] boxnodes Nodes at which box user ICs are set (for each box IC)
    1127                 :            : // *****************************************************************************
    1128                 :            : {
    1129 [ -  + ][ -  - ]:       2786 :   Assert( coord[0].size() == U.nunk(), "Size mismatch" );
         [ -  - ][ -  - ]
    1130                 :            : 
    1131         [ +  - ]:       2786 :   auto ic = IC();
    1132                 :       2786 :   const auto& x = coord[0];
    1133                 :       2786 :   const auto& y = coord[1];
    1134                 :       2786 :   const auto& z = coord[2];
    1135                 :            : 
    1136                 :            :   // Set initial conditions dependeing on problem configured
    1137         [ +  + ]:     385551 :   for (std::size_t i=0; i<x.size(); ++i) {
    1138                 :            : 
    1139                 :            :     // Set background ICs
    1140         [ +  - ]:     382765 :     auto s = ic( x[i], y[i], z[i], t );
    1141 [ -  + ][ -  - ]:     382765 :     Assert( s.size() == U.nprop(), "Size mismatch" );
         [ -  - ][ -  - ]
    1142                 :            : 
    1143                 :            :     // Initialize user-defined ICs in boxes
    1144         [ +  - ]:     382765 :     box( i, s, boxnodes );
    1145                 :            : 
    1146                 :            :     // Set values for ICs
    1147 [ +  - ][ +  + ]:    2324142 :     for (std::size_t c=0; c<s.size(); ++c) U(i,c) = s[c];
    1148                 :            : 
    1149                 :     382765 :   }
    1150                 :       2786 : }
    1151                 :            : 
    1152                 :            : std::function< tk::real( tk::real, tk::real, tk::real ) >
    1153                 :       6198 : PRESSURE_RHS()
    1154                 :            : // *****************************************************************************
    1155                 :            : //  Query user config and assign function to set pressure rhs
    1156                 :            : //! \return The function to call to set pressure rhs
    1157                 :            : // *****************************************************************************
    1158                 :            : {
    1159                 :       6198 :   const auto& problem = inciter::g_cfg.get< tag::problem >();
    1160                 :            : 
    1161                 :       6198 :   std::function< tk::real( tk::real, tk::real, tk::real ) > pr;
    1162                 :            : 
    1163         [ +  + ]:       6198 :   if (problem == "poisson_const")
    1164                 :         22 :     pr = poisson_const::pr;
    1165         [ -  + ]:       6176 :   else if (problem == "poisson_harmonic")
    1166                 :          0 :     pr = poisson_harmonic::pr;
    1167         [ +  + ]:       6176 :   else if (problem == "poisson_sine")
    1168                 :          2 :     pr = poisson_sine::pr;
    1169         [ +  + ]:       6174 :   else if (problem == "poisson_sine3")
    1170                 :          6 :     pr = poisson_sine3::pr;
    1171         [ +  + ]:       6168 :   else if (problem == "poisson_neumann")
    1172                 :          2 :     pr = poisson_neumann::pr;
    1173                 :            : 
    1174                 :       6198 :   return pr;
    1175                 :            : }
    1176                 :            : 
    1177                 :            : std::function< tk::real( tk::real, tk::real, tk::real ) >
    1178                 :       3412 : PRESSURE_IC()
    1179                 :            : // *****************************************************************************
    1180                 :            : //  Query user config and assign function to set pressure initial conditions
    1181                 :            : //! \return The function to call to set pressure initial conditions
    1182                 :            : // *****************************************************************************
    1183                 :            : {
    1184                 :       3412 :   const auto& problem = inciter::g_cfg.get< tag::problem >();
    1185                 :            : 
    1186                 :       3412 :   std::function< tk::real( tk::real, tk::real, tk::real ) > ic;
    1187                 :            : 
    1188         [ +  + ]:       3650 :   if ( problem == "userdef" or
    1189         [ +  - ]:        380 :        problem == "slot_cyl" or
    1190 [ +  + ][ -  + ]:       3792 :        problem == "sheardiff" or
                 [ +  + ]
    1191                 :        142 :        problem == "poiseuille" )
    1192                 :       3270 :     ic = userdef::pic;
    1193         [ +  + ]:        142 :   else if (problem == "poisson_const")
    1194                 :         94 :     ic = poisson_const::ic;
    1195         [ -  + ]:         48 :   else if (problem == "poisson_harmonic")
    1196                 :          0 :     ic = poisson_harmonic::ic;
    1197         [ +  + ]:         48 :   else if (problem == "poisson_sine")
    1198                 :         10 :     ic = poisson_sine::ic;
    1199         [ +  + ]:         38 :   else if (problem == "poisson_sine3")
    1200                 :         28 :     ic = poisson_sine3::ic;
    1201         [ +  - ]:         10 :   else if (problem == "poisson_neumann")
    1202                 :         10 :     ic = poisson_neumann::ic;
    1203                 :            :   else
    1204 [ -  - ][ -  - ]:          0 :     Throw( "problem type not hooked up" );
                 [ -  - ]
    1205                 :            : 
    1206                 :       3412 :   return ic;
    1207                 :          0 : }
    1208                 :            : 
    1209                 :            : std::function< tk::real( tk::real, tk::real, tk::real ) >
    1210                 :       5612 : PRESSURE_SOL()
    1211                 :            : // *****************************************************************************
    1212                 :            : //  Query user config and assign function to query analytic pressure solutions
    1213                 :            : //! \return The function to call to query analytic solutions
    1214                 :            : // *****************************************************************************
    1215                 :            : {
    1216                 :       5612 :   const auto& problem = inciter::g_cfg.get< tag::problem >();
    1217                 :            : 
    1218         [ +  + ]:       6906 :   if ( problem == "userdef" or
    1219         [ +  - ]:       1404 :        problem == "slot_cyl" or
    1220 [ +  + ][ -  + ]:       7016 :        problem == "poiseuille" or
                 [ +  + ]
    1221                 :        110 :        problem == "sheardiff" )
    1222                 :            :   {
    1223                 :       5502 :     return {};
    1224                 :            :   }
    1225                 :            :   else {
    1226                 :        110 :     return PRESSURE_IC();
    1227                 :            :   }
    1228                 :            : }
    1229                 :            : 
    1230                 :            : std::function< std::array< tk::real, 3 >( tk::real, tk::real, tk::real ) >
    1231                 :       6198 : PRESSURE_GRAD()
    1232                 :            : // *****************************************************************************
    1233                 :            : //  Assign function to query pressure gradient at a point
    1234                 :            : //! \return The function to call to query the pressure gradient
    1235                 :            : // *****************************************************************************
    1236                 :            : {
    1237                 :       6198 :   const auto& problem = inciter::g_cfg.get< tag::problem >();
    1238                 :            : 
    1239         [ +  + ]:       6198 :   if (problem == "poisson_neumann")
    1240                 :          2 :     return poisson_neumann::pg;
    1241                 :            : 
    1242                 :       6196 :   return {};
    1243                 :            : }
    1244                 :            : 
    1245                 :            : tk::real
    1246                 :          0 : initialize( tk::real x, tk::real y, tk::real z )
    1247                 :            : // *****************************************************************************
    1248                 :            : //  Evaluate initial condition for pressure
    1249                 :            : //! \param[in] x X coordinate where to evaluate the pressure initial condition
    1250                 :            : //! \param[in] y Y coordinate where to evaluate the pressure initial condition
    1251                 :            : //! \param[in] z Z coordinate where to evaluate the pressure initial condition
    1252                 :            : //! \return Pressure initial condition
    1253                 :            : // *****************************************************************************
    1254                 :            : {
    1255                 :          0 :   const auto& problem = inciter::g_cfg.get< tag::problem >();
    1256                 :            : 
    1257                 :          0 :   std::function< tk::real( tk::real, tk::real, tk::real ) > ic;
    1258                 :            : 
    1259         [ -  - ]:          0 :   if (problem == "poisson_const")
    1260                 :          0 :     ic = poisson_const::ic;
    1261         [ -  - ]:          0 :   else if (problem == "poisson_harmonic")
    1262                 :          0 :     ic = poisson_harmonic::ic;
    1263         [ -  - ]:          0 :   else if (problem == "poisson_sine")
    1264                 :          0 :     ic = poisson_sine::ic;
    1265         [ -  - ]:          0 :   else if (problem == "poisson_sine3")
    1266                 :          0 :     ic = poisson_sine3::ic;
    1267                 :            :   else
    1268 [ -  - ][ -  - ]:          0 :     Throw( "problem type not hooked up" );
                 [ -  - ]
    1269                 :            : 
    1270         [ -  - ]:          0 :   return ic( x, y, z );
    1271                 :          0 : }
    1272                 :            : 
    1273                 :            : std::function< std::vector< tk::real >
    1274                 :            :                  ( tk::real, tk::real, tk::real, tk::real ) >
    1275                 :      78156 : SRC()
    1276                 :            : // *****************************************************************************
    1277                 :            : //  Query user config and assign function to add a source term
    1278                 :            : //! \return The function to call to evaluate a problem-sepcific source term
    1279                 :            : // *****************************************************************************
    1280                 :            : {
    1281                 :      78156 :   const auto& problem = inciter::g_cfg.get< tag::problem >();
    1282                 :            : 
    1283                 :            :   std::function<
    1284                 :      78156 :     std::vector< tk::real >( tk::real, tk::real, tk::real, tk::real ) > src;
    1285                 :            : 
    1286         [ +  + ]:      78156 :   if (problem == "nonlinear_energy_growth")
    1287                 :        676 :     src = nonlinear_energy_growth::src;
    1288         [ +  + ]:      77480 :   else if (problem == "rayleigh_taylor")
    1289                 :       2000 :     src = rayleigh_taylor::src;
    1290         [ +  + ]:      75480 :   else if (problem == "taylor_green")
    1291                 :       3536 :     src = taylor_green::src;
    1292         [ +  + ]:      71944 :   else if (problem == "vortical_flow")
    1293                 :      18547 :     src = vortical_flow::src;
    1294         [ +  + ]:      53397 :   else if (problem == "slot_cyl")
    1295                 :      11537 :     src = slot_cyl::src;
    1296         [ -  + ]:      41860 :   else if (problem == "sheardiff")
    1297                 :          0 :     src = sheardiff::src;
    1298                 :            : 
    1299                 :      78156 :   return src;
    1300                 :            : }
    1301                 :            : 
    1302                 :            : std::function< void( const std::array< std::vector< tk::real >, 3 >&,
    1303                 :            :                      tk::real,
    1304                 :            :                      tk::Fields& ) >
    1305                 :      78156 : PHYS_SRC()
    1306                 :            : // *****************************************************************************
    1307                 :            : //  Query user config and assign function to apply source to numerical solution
    1308                 :            : //! \return The function to call to evaluate a problem-sepcific source term
    1309                 :            : // *****************************************************************************
    1310                 :            : {
    1311                 :      78156 :   const auto& problem = inciter::g_cfg.get< tag::problem >();
    1312                 :            : 
    1313                 :            :   std::function< void( const std::array< std::vector< tk::real >, 3 >&,
    1314                 :            :                        tk::real,
    1315                 :      78156 :                        tk::Fields& ) > src;
    1316                 :            : 
    1317         [ +  + ]:      78156 :   if (problem == "point_src") {
    1318                 :        600 :     src = point_src::src;
    1319                 :            :   }
    1320                 :            : 
    1321                 :      78156 :   return src;
    1322                 :            : }
    1323                 :            : 
    1324                 :            : } // problems::

Generated by: LCOV version 1.16