Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/Transport/Problem/SlotCyl.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 "SlotCyl.hpp"
19 : :
20 : : using inciter::TransportProblemSlotCyl;
21 : :
22 : : std::vector< tk::real >
23 : 1613744 : TransportProblemSlotCyl::initialize( ncomp_t ncomp,
24 : : const std::vector< EOS >&, tk::real x, tk::real y, tk::real,
25 : : tk::real t )
26 : : // *****************************************************************************
27 : : // Evaluate analytical solution at (x,y,t) for all components
28 : : //! \param[in] ncomp Number of components in this transport equation system
29 : : //! \param[in] x X coordinate where to evaluate the solution
30 : : //! \param[in] y Y coordinate where to evaluate the solution
31 : : //! \param[in] t Time where to evaluate the solution
32 : : //! \return Values of all components evaluated at (x,y,t)
33 : : // *****************************************************************************
34 : : {
35 : : using std::sin; using std::cos;
36 : :
37 [ + - ]: 1613744 : std::vector< tk::real > s( ncomp, 0.0 );
38 [ + + ]: 3227488 : for (ncomp_t c=0; c<ncomp; ++c) {
39 : 1613744 : auto T = t + 2.0*M_PI/static_cast<tk::real>(ncomp)*static_cast<tk::real>(c);
40 : 1613744 : const tk::real R0 = 0.15;
41 : :
42 : : // center of the cone
43 : 1613744 : tk::real x0 = 0.5;
44 : 1613744 : tk::real y0 = 0.25;
45 : 1613744 : tk::real r = std::sqrt((x0-0.5)*(x0-0.5) + (y0-0.5)*(y0-0.5));
46 : 1613744 : tk::real kx = 0.5 + r*sin( T );
47 : 1613744 : tk::real ky = 0.5 - r*cos( T );
48 : :
49 : : // center of the hump
50 : 1613744 : x0 = 0.25;
51 : 1613744 : y0 = 0.5;
52 : 1613744 : r = std::sqrt((x0-0.5)*(x0-0.5) + (y0-0.5)*(y0-0.5));
53 : 1613744 : tk::real hx = 0.5 + r*sin( T-M_PI/2.0 ),
54 : 1613744 : hy = 0.5 - r*cos( T-M_PI/2.0 );
55 : :
56 : : // center of the slotted cylinder
57 : 1613744 : x0 = 0.5;
58 : 1613744 : y0 = 0.75;
59 : 1613744 : r = std::sqrt((x0-0.5)*(x0-0.5) + (y0-0.5)*(y0-0.5));
60 : 1613744 : tk::real cx = 0.5 + r*sin( T+M_PI ),
61 : 1613744 : cy = 0.5 - r*cos( T+M_PI );
62 : :
63 : : // end points of the cylinder slot
64 : 1613744 : tk::real i1x = 0.525, i1y = cy - r*cos( std::asin(0.025/r) ),
65 : 1613744 : i2x = 0.525, i2y = 0.8,
66 : 1613744 : i3x = 0.475, i3y = 0.8;
67 : :
68 : : // rotate end points of cylinder slot
69 : 1613744 : tk::real ri1x = 0.5 + cos(T)*(i1x-0.5) - sin(T)*(i1y-0.5),
70 : 1613744 : ri1y = 0.5 + sin(T)*(i1x-0.5) + cos(T)*(i1y-0.5),
71 : 1613744 : ri2x = 0.5 + cos(T)*(i2x-0.5) - sin(T)*(i2y-0.5),
72 : 1613744 : ri2y = 0.5 + sin(T)*(i2x-0.5) + cos(T)*(i2y-0.5),
73 : 1613744 : ri3x = 0.5 + cos(T)*(i3x-0.5) - sin(T)*(i3y-0.5),
74 : 1613744 : ri3y = 0.5 + sin(T)*(i3x-0.5) + cos(T)*(i3y-0.5);
75 : :
76 : : // direction of slot sides
77 : 1613744 : tk::real v1x = ri2x-ri1x, v1y = ri2y-ri1y,
78 : 1613744 : v2x = ri3x-ri2x, v2y = ri3y-ri2y;
79 : :
80 : : // lengths of direction of slot sides vectors
81 : 1613744 : tk::real v1 = std::sqrt(v1x*v1x + v1y*v1y),
82 : 1613744 : v2 = std::sqrt(v2x*v2x + v2y*v2y);
83 : :
84 : : // cone
85 : 1613744 : r = std::sqrt((x-kx)*(x-kx) + (y-ky)*(y-ky)) / R0;
86 [ + + ]: 1613744 : if (r<1.0) s[c] = 0.6*(1.0-r);
87 : :
88 : : // hump
89 : 1613744 : r = std::sqrt((x-hx)*(x-hx) + (y-hy)*(y-hy)) / R0;
90 [ + + ]: 1613744 : if (r<1.0) s[c] = 0.2*(1.0+cos(M_PI*std::min(r,1.0)));
91 : :
92 : : // cylinder
93 : 1613744 : r = std::sqrt((x-cx)*(x-cx) + (y-cy)*(y-cy)) / R0;
94 : 1613744 : const std::array< tk::real, 2 > r1{{ v1x, v1y }},
95 : 1613744 : r2{{ x-ri1x, y-ri1y }};
96 : 1613744 : const auto d1 = (r1[0]*r2[1] - r2[0]*r1[1]) / v1;
97 : 1613744 : const std::array< tk::real, 2 > r3{{ v2x, v2y }},
98 : 1613744 : r4{{ x-ri2x, y-ri2y }};
99 : 1613744 : const auto d2 = (r3[0]*r4[1] - r4[0]*r3[1]) / v2;
100 [ + + ][ + + ]: 1613744 : if (r<1.0 && (d1>0.05 || d1<0.0 || d2<0.0)) s[c] = 0.6;
[ + + ][ + + ]
101 : : }
102 : :
103 : 1613744 : return s;
104 : : }
105 : :
106 : : std::vector< std::array< tk::real, 3 > >
107 : 8671479 : TransportProblemSlotCyl::prescribedVelocity( ncomp_t ncomp,
108 : : tk::real x, tk::real y, tk::real,
109 : : tk::real )
110 : : // *****************************************************************************
111 : : // Assign prescribed shear velocity at a point
112 : : //! \param[in] ncomp Number of components in this transport equation
113 : : //! \param[in] x X coordinate at which to assign velocity
114 : : //! \param[in] y y coordinate at which to assign velocity
115 : : //! \return Velocity assigned to all vertices of a tetrehedron, size:
116 : : //! ncomp * ndim = [ncomp][3]
117 : : // *****************************************************************************
118 : : {
119 [ + - ]: 8671479 : std::vector< std::array< tk::real, 3 > > vel( ncomp );
120 : :
121 [ + + ]: 17342958 : for (ncomp_t c=0; c<ncomp; ++c)
122 : 8671479 : vel[c] = {{ 0.5-y, x-0.5, 0.0 }};
123 : :
124 : 8671479 : return vel;
125 : : }
|
