Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/Integrate/SolidTerms.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 Functions for computing integrals for the RHS terms
9 : : of the inverse deformation equations in the multi-material
10 : : solid dynamics solver
11 : : \details This file contains routines to integrate right hand side terms
12 : : to the inverse deformation equations for the multi-material
13 : : solid dynamics solver.
14 : : */
15 : : // *****************************************************************************
16 : :
17 : : #include <array>
18 : : #include <vector>
19 : :
20 : : #include "SolidTerms.hpp"
21 : : #include "Vector.hpp"
22 : : #include "Quadrature.hpp"
23 : : #include "Reconstruction.hpp"
24 : : #include "MultiMatTerms.hpp"
25 : : #include "MultiMat/MultiMatIndexing.hpp"
26 : : #include "MultiMat/MiscMultiMatFns.hpp"
27 : : #include "Inciter/InputDeck/InputDeck.hpp"
28 : :
29 : : namespace inciter {
30 : : extern ctr::InputDeck g_inputdeck;
31 : : }
32 : :
33 : : namespace tk {
34 : :
35 : : void
36 : 0 : solidTermsVolInt(
37 : : std::size_t nmat,
38 : : const std::vector< inciter::EOS >& mat_blk,
39 : : const std::size_t ndof,
40 : : const std::size_t rdof,
41 : : const std::size_t nelem,
42 : : const std::vector< std::size_t >& inpoel,
43 : : const UnsMesh::Coords& coord,
44 : : const Fields& geoElem,
45 : : const Fields& U,
46 : : const Fields& P,
47 : : const std::vector< std::size_t >& ndofel,
48 : : const tk::real dt,
49 : : Fields& R,
50 : : int intsharp )
51 : : // *****************************************************************************
52 : : // Compute all RHS volume terms in the inverse deformation equations to
53 : : // satisfy the condition curl(g) = 0.
54 : : //! \param[in] nmat Number of materials in this PDE system
55 : : //! \param[in] mat_blk EOS material block
56 : : //! \param[in] ndof Maximum number of degrees of freedom
57 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
58 : : //! \param[in] nelem Maximum number of elements
59 : : //! \param[in] inpoel Element-node connectivity
60 : : //! \param[in] coord Array of nodal coordinates
61 : : //! \param[in] geoElem Element geometry array
62 : : //! \param[in] U Solution vector at recent time step
63 : : //! \param[in] P Vector of primitives at recent time step
64 : : //! \param[in] ndofel Vector of local number of degrees of freedom
65 : : //! \param[in] dt Delta time
66 : : //! \param[in,out] R Right-hand side vector computed
67 : : //! \param[in] intsharp Interface compression tag, an optional argument, with
68 : : //! default 0, so that it is unused for single-material and transport.
69 : : // *****************************************************************************
70 : : {
71 : :
72 : : using inciter::velocityDofIdx;
73 : : using inciter::volfracDofIdx;
74 : : using inciter::deformDofIdx;
75 : :
76 : : const auto& solidx =
77 : 0 : inciter::g_inputdeck.get< tag::matidxmap, tag::solidx >();
78 : :
79 : 0 : auto ncomp = R.nprop()/ndof;
80 : 0 : auto nprim = P.nprop()/rdof;
81 : 0 : const auto& cx = coord[0];
82 : 0 : const auto& cy = coord[1];
83 : 0 : const auto& cz = coord[2];
84 : :
85 [ - - ]: 0 : for (std::size_t e=0; e<nelem; ++e)
86 : : {
87 : :
88 : : // Relaxation and diffusion coefficients
89 [ - - ]: 0 : tk::real dx = geoElem(0,4);
90 : 0 : tk::real D = dx*dx/(12.0*dt) * 1.0e-00;
91 : 0 : tk::real eta = 1.0/(6.0*dt) * 1.0e+00;
92 : :
93 [ - - ]: 0 : auto ng = tk::NGvol(ndofel[e]);
94 : :
95 : : // arrays for quadrature points
96 : 0 : std::array< std::vector< real >, 3 > coordgp;
97 : 0 : std::vector< real > wgp;
98 : :
99 [ - - ]: 0 : coordgp[0].resize( ng );
100 [ - - ]: 0 : coordgp[1].resize( ng );
101 [ - - ]: 0 : coordgp[2].resize( ng );
102 [ - - ]: 0 : wgp.resize( ng );
103 : :
104 [ - - ]: 0 : GaussQuadratureTet( ng, coordgp, wgp );
105 : :
106 : : // Extract the element coordinates
107 : : std::array< std::array< real, 3>, 4 > coordel {{
108 : 0 : {{ cx[ inpoel[4*e ] ], cy[ inpoel[4*e ] ], cz[ inpoel[4*e ] ] }},
109 : 0 : {{ cx[ inpoel[4*e+1] ], cy[ inpoel[4*e+1] ], cz[ inpoel[4*e+1] ] }},
110 : 0 : {{ cx[ inpoel[4*e+2] ], cy[ inpoel[4*e+2] ], cz[ inpoel[4*e+2] ] }},
111 : 0 : {{ cx[ inpoel[4*e+3] ], cy[ inpoel[4*e+3] ], cz[ inpoel[4*e+3] ] }}}};
112 : :
113 : : // Compute the inverse of the Jacobian
114 : : auto jacInv =
115 : 0 : tk::inverseJacobian( coordel[0], coordel[1], coordel[2], coordel[3] );
116 : :
117 [ - - ]: 0 : for (std::size_t igp=0; igp<ng; ++igp)
118 : : {
119 : : // Compute the coordinates of quadrature point at physical domain
120 [ - - ]: 0 : auto gp = eval_gp( igp, coordel, coordgp );
121 : :
122 : : // Compute the basis function
123 : : auto B =
124 [ - - ]: 0 : eval_basis(rdof,coordgp[0][igp],coordgp[1][igp],coordgp[2][igp]);
125 : :
126 : : // Compute derivatives
127 : 0 : std::array< std::vector<tk::real>, 3 > dBdx;
128 [ - - ]: 0 : dBdx[0].resize( ndof, 0 );
129 [ - - ]: 0 : dBdx[1].resize( ndof, 0 );
130 [ - - ]: 0 : dBdx[2].resize( ndof, 0 );
131 [ - - ]: 0 : if (rdof > 1)
132 : : {
133 [ - - ]: 0 : dBdx = tk::eval_dBdx_p1( rdof, jacInv );
134 [ - - ]: 0 : if(ndof > 4) {
135 [ - - ]: 0 : tk::eval_dBdx_p2(igp, coordgp, jacInv, dBdx);
136 : : }
137 : : }
138 : :
139 : : // Get state
140 : 0 : std::vector< real > state;
141 [ - - ]: 0 : state = evalPolynomialSol(mat_blk, intsharp, ncomp, nprim, rdof,
142 : 0 : nmat, e, rdof, inpoel, coord, geoElem, gp, B, U, P);
143 : :
144 : : // Get velocity
145 : 0 : std::vector< real > v = {{state[ncomp+inciter::velocityIdx(nmat, 0)],
146 : 0 : state[ncomp+inciter::velocityIdx(nmat, 1)],
147 [ - - ]: 0 : state[ncomp+inciter::velocityIdx(nmat, 2)]}};
148 : :
149 : : // Loop through materials
150 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
151 : : {
152 [ - - ]: 0 : if (solidx[k] > 0)
153 : : {
154 : 0 : tk::real alpha = state[inciter::volfracIdx(nmat, k)];
155 : : std::array< std::array< tk::real, 3 >, 3 > g;
156 : : // Compute the source terms
157 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
158 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
159 : 0 : g[i][j] = state[inciter::deformIdx(nmat,solidx[k],i,j)]/alpha;
160 : :
161 : : // ** HARDCODED **
162 : : // Should be able to store rho_0 for every cell at every Gauss point.
163 : 0 : tk::real rho = state[inciter::densityIdx(nmat, k)];
164 : : tk::real rho0;
165 [ - - ]: 0 : if (k==0) rho0 = 8900.0;
166 : 0 : else rho0 = 2700.0;
167 : 0 : tk::real rfact = eta*(rho/(rho0*tk::determinant(g))-1.0);
168 : :
169 : : // Compute the source terms
170 [ - - ]: 0 : std::vector< real > s(9*ndof, 0.0);
171 [ - - ]: 0 : std::vector< real > deriv(6, 0.0);
172 [ - - ]: 0 : std::vector< std::size_t > defIdx(6, 0);
173 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
174 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
175 [ - - ]: 0 : for (std::size_t idof=0; idof<ndof; ++idof)
176 : : {
177 [ - - ]: 0 : for (std::size_t l=0; l<deriv.size(); ++l)
178 : 0 : deriv[l] = 0.0;
179 [ - - ]: 0 : for (std::size_t jdof=0; jdof<rdof; ++jdof)
180 : : {
181 : : // Find indeces for all unknowns used
182 : 0 : defIdx[0] = deformDofIdx(nmat,solidx[k],i,(j+1)%3,rdof,jdof);
183 : 0 : defIdx[1] = deformDofIdx(nmat,solidx[k],i, j,rdof,jdof);
184 : 0 : defIdx[2] = deformDofIdx(nmat,solidx[k],i, j,rdof,jdof);
185 : 0 : defIdx[3] = deformDofIdx(nmat,solidx[k],i,(j+2)%3,rdof,jdof);
186 : 0 : defIdx[4] = volfracDofIdx(nmat,k,rdof,jdof);
187 : 0 : defIdx[5] = volfracDofIdx(nmat,k,rdof,jdof);
188 : : // Compute derivatives
189 [ - - ]: 0 : deriv[0] += U(e,defIdx[0]) *dBdx[ j][jdof];
190 [ - - ]: 0 : deriv[1] += U(e,defIdx[1]) *dBdx[(j+1)%3][jdof];
191 [ - - ]: 0 : deriv[2] += U(e,defIdx[2]) *dBdx[(j+2)%3][jdof];
192 [ - - ]: 0 : deriv[3] += U(e,defIdx[3]) *dBdx[ j][jdof];
193 [ - - ]: 0 : deriv[4] += U(e,defIdx[4]) *dBdx[(j+1)%3][jdof];
194 [ - - ]: 0 : deriv[5] += U(e,defIdx[5]) *dBdx[(j+2)%3][jdof];
195 : : }
196 : 0 : s[(i*3+j)*ndof+idof] = B[idof]*rfact*alpha*g[i][j]
197 : 0 : + alpha*B[idof]*(v[(j+1)%3]*(deriv[0]-deriv[1])
198 : 0 : -v[(j+2)%3]*(deriv[2]-deriv[3]))
199 : 0 : + D*((alpha*dBdx[(j+1)%3][idof]+B[idof]*deriv[4])
200 : 0 : *(deriv[0]-deriv[1])
201 : 0 : -(alpha*dBdx[(j+2)%3][idof]+B[idof]*deriv[5])
202 : 0 : *(deriv[2]-deriv[3]));
203 : : }
204 : :
205 [ - - ]: 0 : auto wt = wgp[igp] * geoElem(e, 0);
206 : :
207 [ - - ]: 0 : update_rhs( nmat, ndof, wt, e, solidx[k], s, R );
208 : :
209 : : }
210 : : }
211 : : }
212 : : }
213 : :
214 : 0 : }
215 : :
216 : : void
217 : 0 : solidTermsSurfInt( std::size_t nmat,
218 : : const std::size_t ndof,
219 : : const std::size_t rdof,
220 : : const std::array< tk::real, 3 >& fn,
221 : : const std::size_t el,
222 : : const std::size_t er,
223 : : const std::vector< std::size_t >& solidx,
224 : : const Fields& geoElem,
225 : : const Fields& U,
226 : : const std::array< std::array< tk::real, 3>, 4 > coordel_l,
227 : : const std::array< std::array< tk::real, 3>, 4 > coordel_r,
228 : : const std::size_t igp,
229 : : const std::array< std::vector< tk::real >, 2 >& coordgp,
230 : : const tk::real dt,
231 : : std::vector< tk::real >& fl )
232 : : // *****************************************************************************
233 : : // Compute all RHS surface terms in the inverse deformation equations to
234 : : // satisfy the condition curl(g) = 0.
235 : : //! \param[in] nmat Number of materials in this PDE system
236 : : //! \param[in] ndof Maximum number of degrees of freedom
237 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
238 : : //! \param[in] fn Face/Surface normal
239 : : //! \param[in] el Left element index
240 : : //! \param[in] er Right element index
241 : : //! \param[in] B_l Basis function for the left element
242 : : //! \param[in] B_r Basis function for the right element
243 : : //! \param[in] solidx Solid material indicator
244 : : //! \param[in] geoElem Element geometry array
245 : : //! \param[in] U Solution vector at recent time step
246 : : //! \param[in] coordel_l Coordinates of left elements' nodes
247 : : //! \param[in] coordel_r Coordinates of right elements' nodes
248 : : //! \param[in] igp Index of quadrature points
249 : : //! \param[in] coordgp Gauss point coordinates for tetrahedron element
250 : : //! \param[in] dt Delta time
251 : : //! \param[in,out] fl Surface flux
252 : : // *****************************************************************************
253 : : {
254 : :
255 : : using inciter::deformIdx;
256 : : using inciter::volfracDofIdx;
257 : : using inciter::deformDofIdx;
258 : :
259 : : // Diffusion coefficient
260 : : //tk::real dx = std::min(geoElem(el,4),geoElem(er,4));
261 [ - - ]: 0 : tk::real dx = geoElem(0,4);
262 : 0 : tk::real D = dx*dx/(12.0*dt) * 1.0e-00;
263 : : // Compute the inverse of the Jacobian
264 : : auto jacInv_l =
265 : 0 : tk::inverseJacobian(coordel_l[0],coordel_l[1],coordel_l[2],coordel_l[3]);
266 : : auto jacInv_r =
267 : 0 : tk::inverseJacobian(coordel_r[0],coordel_r[1],coordel_r[2],coordel_r[3]);
268 : : // Compute derivatives
269 : 0 : std::array< std::vector<tk::real>, 3 > dBdx_l, dBdx_r;
270 [ - - ]: 0 : dBdx_l[0].resize( ndof, 0 );
271 [ - - ]: 0 : dBdx_l[1].resize( ndof, 0 );
272 [ - - ]: 0 : dBdx_l[2].resize( ndof, 0 );
273 [ - - ]: 0 : dBdx_r[0].resize( ndof, 0 );
274 [ - - ]: 0 : dBdx_r[1].resize( ndof, 0 );
275 [ - - ]: 0 : dBdx_r[2].resize( ndof, 0 );
276 [ - - ]: 0 : if (rdof > 1)
277 : : {
278 [ - - ]: 0 : dBdx_l = tk::eval_dBdx_p1( rdof, jacInv_l );
279 [ - - ]: 0 : dBdx_r = tk::eval_dBdx_p1( rdof, jacInv_r );
280 [ - - ]: 0 : if(ndof > 4) {
281 : : // Since eval_dBdx_p2 takes a coordgp of dimension ng by 3
282 : : // I need to add a row of zeros to my coordgp (UNTESTED)
283 : 0 : std::array< std::vector< tk::real >, 3 > coordgp_aux;
284 : 0 : coordgp_aux[0][igp] = coordgp[0][igp];
285 : 0 : coordgp_aux[1][igp] = coordgp[1][igp];
286 : 0 : coordgp_aux[2][igp] = 0.0;
287 [ - - ]: 0 : tk::eval_dBdx_p2(igp, coordgp_aux, jacInv_l, dBdx_l);
288 [ - - ]: 0 : tk::eval_dBdx_p2(igp, coordgp_aux, jacInv_r, dBdx_r);
289 : : }
290 : : }
291 : : // Loop through materials
292 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
293 : : {
294 [ - - ]: 0 : if (solidx[k] > 0)
295 : : {
296 : : // Compute all derivatives
297 : : std::array< std::array< std::array< tk::real, 3 >, 3 >, 3 >
298 : : deriv_l, deriv_r;
299 : : std::size_t defIdx;
300 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
301 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
302 : : {
303 : 0 : deriv_l[0][i][j] = 0.0;
304 : 0 : deriv_l[1][i][j] = 0.0;
305 : 0 : deriv_l[2][i][j] = 0.0;
306 : 0 : deriv_r[0][i][j] = 0.0;
307 : 0 : deriv_r[1][i][j] = 0.0;
308 : 0 : deriv_r[2][i][j] = 0.0;
309 [ - - ]: 0 : for (std::size_t jdof=0; jdof<rdof; ++jdof)
310 : : {
311 : 0 : defIdx = deformDofIdx(nmat,solidx[k],i,j,rdof,jdof);
312 [ - - ]: 0 : deriv_l[0][i][j] += U(el,defIdx)*dBdx_l[0][jdof];
313 [ - - ]: 0 : deriv_l[1][i][j] += U(el,defIdx)*dBdx_l[1][jdof];
314 [ - - ]: 0 : deriv_l[2][i][j] += U(el,defIdx)*dBdx_l[2][jdof];
315 [ - - ]: 0 : deriv_r[0][i][j] += U(er,defIdx)*dBdx_r[0][jdof];
316 [ - - ]: 0 : deriv_r[1][i][j] += U(er,defIdx)*dBdx_r[1][jdof];
317 [ - - ]: 0 : deriv_r[2][i][j] += U(er,defIdx)*dBdx_r[2][jdof];
318 : : }
319 : : }
320 : :
321 : : // Compute the source terms
322 [ - - ]: 0 : tk::real alpha_l = U(el,volfracDofIdx(nmat,k,rdof,0));
323 [ - - ]: 0 : tk::real alpha_r = U(er,volfracDofIdx(nmat,k,rdof,0));
324 : : tk::real sl, sr;
325 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
326 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
327 : : {
328 : : // Compute source
329 : 0 : sl = alpha_l*D*(fn[(j+1)%3]*(deriv_l[j][i][(j+1)%3]
330 : 0 : -deriv_l[(j+1)%3][i][j])
331 : 0 : -fn[(j+2)%3]*(deriv_l[(j+2)%3][i][j]
332 : 0 : -deriv_l[j][i][(j+2)%3]));
333 : 0 : sr = alpha_r*D*(fn[(j+1)%3]*(deriv_r[j][i][(j+1)%3]
334 : 0 : -deriv_r[(j+1)%3][i][j])
335 : 0 : -fn[(j+2)%3]*(deriv_r[(j+2)%3][i][j]
336 : 0 : -deriv_r[j][i][(j+2)%3]));
337 : : // Add to flux
338 : 0 : fl[deformIdx(nmat,solidx[k],i,j)] += 0.5*(sl+sr);
339 : : }
340 : : }
341 : : }
342 : 0 : }
343 : :
344 : : void
345 : 0 : update_rhs( std::size_t nmat,
346 : : const std::size_t ndof,
347 : : const tk::real wt,
348 : : const std::size_t e,
349 : : const std::size_t solidx,
350 : : const std::vector< real >& s,
351 : : Fields& R )
352 : : // *****************************************************************************
353 : : // Update the rhs by adding the source term integrals
354 : : //! \param[in] nmat Number of materials in this PDE system
355 : : //! \param[in] rdof Maximum number of degrees of freedom
356 : : //! \param[in] wt Weight of gauss quadrature point
357 : : //! \param[in] e Element index
358 : : //! \param[in] solidx Material index indicator
359 : : //! \param[in] s Source term vector
360 : : //! \param[in,out] R Right-hand side vector computed
361 : : // *****************************************************************************
362 : : {
363 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
364 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
365 [ - - ]: 0 : for (std::size_t idof=0; idof<ndof; ++idof)
366 : : {
367 : 0 : std::size_t dofId = inciter::deformDofIdx(nmat,solidx,i,j,ndof,idof);
368 : 0 : std::size_t srcId = (i*3+j)*ndof+idof;
369 : 0 : R(e, dofId) += wt * s[srcId];
370 : : }
371 : 0 : }
372 : :
373 : : } // tk::
|
