Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/Transport/Problem/ShearDiff.cpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.
7 : : All rights reserved. See the LICENSE file for details.
8 : : \brief Problem configuration for transport equations
9 : : \details This file defines a Problem policy class for the transport
10 : : equations, defined in PDE/Transport/CGTransport.h implementing
11 : : node-centered continuous Galerkin (CG) and PDE/Transport/DGTransport.h
12 : : implementing cell-centered discontinuous Galerkin (DG) discretizations.
13 : : See PDE/Transport/Problem.h for general requirements on Problem policy
14 : : classes for cg::Transport and dg::Transport.
15 : : */
16 : : // *****************************************************************************
17 : :
18 : : #include "ShearDiff.hpp"
19 : : #include "Inciter/InputDeck/InputDeck.hpp"
20 : :
21 : : namespace inciter {
22 : :
23 : : extern ctr::InputDeck g_inputdeck;
24 : :
25 : : } // ::inciter
26 : :
27 : : using inciter::TransportProblemShearDiff;
28 : :
29 : : std::vector< tk::real >
30 : 0 : TransportProblemShearDiff::initialize( ncomp_t ncomp,
31 : : const std::vector< EOS >&, tk::real x, tk::real y, tk::real z,
32 : : tk::real t )
33 : : // *****************************************************************************
34 : : // Evaluate analytical solution at (x,y,z,t) for all components
35 : : //! \param[in] ncomp Number of components in this transport equation system
36 : : //! \param[in] x X coordinate where to evaluate the solution
37 : : //! \param[in] y Y coordinate where to evaluate the solution
38 : : //! \param[in] z Z coordinate where to evaluate the solution
39 : : //! \param[in] t Time where to evaluate the solution
40 : : //! \return Values of all components evaluated at (x,y,t)
41 : : // *****************************************************************************
42 : : {
43 : 0 : const auto& u0 = g_inputdeck.get< eq, tag::u0 >();
44 : 0 : const auto& d = g_inputdeck.get< eq, tag::diffusivity >();
45 : 0 : const auto& l = g_inputdeck.get< eq, tag::lambda >();
46 : :
47 [ - - ]: 0 : std::vector< tk::real > r( ncomp );
48 [ - - ]: 0 : for (ncomp_t c=0; c<ncomp; ++c) {
49 : 0 : const auto li = 2*c;
50 : 0 : const auto di = 3*c;
51 : 0 : const auto phi3s = (l[li+0]*l[li+0]*d[di+1]/d[di+0] +
52 : 0 : l[li+1]*l[li+1]*d[di+2]/d[di+0]) / 12.0;
53 : 0 : r[c] =
54 : 0 : 1.0 / ( 8.0 * std::pow(M_PI,3.0/2.0) *
55 : 0 : std::sqrt(d[di+0]*d[di+1]*d[di+2]) *
56 : 0 : std::pow(t,3.0/2.0) * std::sqrt(1.0+phi3s*t*t) ) *
57 : 0 : exp( -std::pow( x - u0[c]*t -
58 : 0 : 0.5*(l[li+0]*y + l[li+1]*z)*t, 2.0 ) /
59 : 0 : ( 4.0 * d[di+0] * t * (1.0 + phi3s*t*t) )
60 : 0 : -y*y / ( 4.0 * d[di+1] * t )
61 : 0 : -z*z / ( 4.0 * d[di+2] * t ) );
62 : : }
63 : :
64 : 0 : return r;
65 : : }
66 : :
67 : : void
68 : 0 : TransportProblemShearDiff::errchk( ncomp_t ncomp ) const
69 : : // *****************************************************************************
70 : : // Do error checking on PDE parameters
71 : : //! \param[in] ncomp Number of components in this transport equation
72 : : // *****************************************************************************
73 : : {
74 [ - - ]: 0 : auto u0 = g_inputdeck.get< eq, tag::u0 >();
75 [ - - ][ - - ]: 0 : ErrChk( ncomp == u0.size(),
[ - - ][ - - ]
76 : : "Wrong number of advection-diffusion PDE parameters 'u0'" );
77 : :
78 [ - - ]: 0 : auto lambda = g_inputdeck.get< eq, tag::lambda >();
79 [ - - ][ - - ]: 0 : ErrChk( 2*ncomp == lambda.size(),
[ - - ][ - - ]
80 : : "Wrong number of advection-diffusion PDE parameters 'lambda'" );
81 : :
82 : 0 : auto& d = g_inputdeck.get< eq, tag::diffusivity >();
83 [ - - ][ - - ]: 0 : ErrChk( 3*ncomp == d.size(),
[ - - ][ - - ]
84 : : "Wrong number of advection-diffusion PDE parameters 'diffusivity'" );
85 : 0 : }
86 : :
87 : : std::vector< std::array< tk::real, 3 > >
88 : 0 : TransportProblemShearDiff::prescribedVelocity( ncomp_t ncomp,
89 : : tk::real, tk::real y, tk::real z,
90 : : tk::real )
91 : : // *****************************************************************************
92 : : // Assign prescribed shear velocity at a point
93 : : //! \param[in] ncomp Number of components in this transport equation
94 : : //! \param[in] y y coordinate at which to assign velocity
95 : : //! \param[in] z Z coordinate at which to assign velocity
96 : : //! \return Velocity assigned to all vertices of a tetrehedron, size:
97 : : //! ncomp * ndim = [ncomp][3]
98 : : // *****************************************************************************
99 : : {
100 [ - - ]: 0 : auto u0 = g_inputdeck.get< eq, tag::u0 >();
101 [ - - ]: 0 : auto l = g_inputdeck.get< eq, tag::lambda >();
102 : :
103 [ - - ]: 0 : std::vector< std::array< tk::real, 3 > > vel( ncomp );
104 [ - - ]: 0 : for (ncomp_t c=0; c<ncomp; ++c)
105 : 0 : vel[c] = {{ u0[c] + l[2*c+0]*y + l[2*c+1]*z, 0.0, 0.0 }};
106 : :
107 : 0 : return vel;
108 : : }
|
