David Galiffi
0403aaa97f
* Updating clang-format to v18
- Updates the pre-commit-config
- Formats source files according to the utility
Signed-off-by: David Galiffi <David.Galiffi@amd.com>
* Update format source workflow
Signed-off-by: David Galiffi <David.Galiffi@amd.com>
* Update CONTRIBUTING
* Update comment in .clang-format
* Update CONTRIBUTING.md
* Update helper script
---------
Signed-off-by: David Galiffi <David.Galiffi@amd.com>
[ROCm/rocprofiler-systems commit: 1e13b590e7]
2309 خطوط
78 KiB
C++
2309 خطوط
78 KiB
C++
|
|
#include <climits>
|
|
#include <ctype.h>
|
|
#include <iostream>
|
|
#include <math.h>
|
|
#include <stdio.h>
|
|
#include <stdlib.h>
|
|
#include <string.h>
|
|
#include <sys/time.h>
|
|
#include <time.h>
|
|
#include <unistd.h>
|
|
#include <vector>
|
|
|
|
#include "lulesh.h"
|
|
|
|
#include "causal.hpp"
|
|
|
|
static Kokkos::View<Real_t*> buffer;
|
|
static size_t buffer_size;
|
|
static size_t buffer_offset;
|
|
static int do_atomic;
|
|
|
|
void
|
|
ResizeBuffer(const size_t size)
|
|
{
|
|
buffer_offset = 0;
|
|
if(size / sizeof(Real_t) + 1 > buffer_size)
|
|
{
|
|
buffer_size = size / sizeof(Real_t) + 1;
|
|
buffer = Kokkos::View<Real_t*>("Buffer", buffer_size);
|
|
}
|
|
}
|
|
|
|
template <class Type>
|
|
Type*
|
|
AllocateFromBuffer(const Index_t& count)
|
|
{
|
|
const Index_t offset = (count * sizeof(Type) + sizeof(Real_t) - 1) / sizeof(Real_t);
|
|
Real_t* ptr = buffer.data() + buffer_offset;
|
|
buffer_offset += ((offset + 511) / 512) * 512;
|
|
return static_cast<Type*>(ptr);
|
|
}
|
|
|
|
static inline void
|
|
TimeIncrement(Domain& domain)
|
|
{
|
|
Real_t targetdt = domain.stoptime() - domain.time();
|
|
|
|
if((domain.dtfixed() <= Real_t(0.0)) && (domain.cycle() != Int_t(0)))
|
|
{
|
|
// CAUSAL_BEGIN("TimeIncrement_Iteration")
|
|
Real_t ratio;
|
|
Real_t olddt = domain.deltatime();
|
|
|
|
Real_t gnewdt = Real_t(1.0e+20);
|
|
Real_t newdt;
|
|
if(domain.dtcourant() < gnewdt)
|
|
{
|
|
gnewdt = domain.dtcourant() / Real_t(2.0);
|
|
}
|
|
if(domain.dthydro() < gnewdt)
|
|
{
|
|
gnewdt = domain.dthydro() * Real_t(2.0) / Real_t(3.0);
|
|
}
|
|
|
|
#if USE_MPI
|
|
MPI_Allreduce(&gnewdt, &newdt, 1,
|
|
((sizeof(Real_t) == 4) ? MPI_FLOAT : MPI_DOUBLE), MPI_MIN,
|
|
MPI_COMM_WORLD);
|
|
#else
|
|
newdt = gnewdt;
|
|
#endif
|
|
|
|
ratio = newdt / olddt;
|
|
if(ratio >= Real_t(1.0))
|
|
{
|
|
if(ratio < domain.deltatimemultlb())
|
|
{
|
|
newdt = olddt;
|
|
}
|
|
else if(ratio > domain.deltatimemultub())
|
|
{
|
|
newdt = olddt * domain.deltatimemultub();
|
|
}
|
|
}
|
|
|
|
if(newdt > domain.dtmax())
|
|
{
|
|
newdt = domain.dtmax();
|
|
}
|
|
domain.deltatime() = newdt;
|
|
CAUSAL_PROGRESS_NAMED("TimeIncrement_Iteration");
|
|
// CAUSAL_END("TimeIncrement_Iteration")
|
|
}
|
|
|
|
if((targetdt > domain.deltatime()) &&
|
|
(targetdt < (Real_t(4.0) * domain.deltatime() / Real_t(3.0))))
|
|
{
|
|
targetdt = Real_t(2.0) * domain.deltatime() / Real_t(3.0);
|
|
}
|
|
|
|
if(targetdt < domain.deltatime())
|
|
{
|
|
domain.deltatime() = targetdt;
|
|
}
|
|
|
|
domain.time() += domain.deltatime();
|
|
|
|
++domain.cycle();
|
|
}
|
|
|
|
KOKKOS_INLINE_FUNCTION void
|
|
CollectDomainNodesToElemNodes(const Domain& domain, const Index_t* elemToNode,
|
|
Real_t elemX[8], Real_t elemY[8], Real_t elemZ[8])
|
|
{
|
|
Index_t nd0i = elemToNode[0];
|
|
Index_t nd1i = elemToNode[1];
|
|
Index_t nd2i = elemToNode[2];
|
|
Index_t nd3i = elemToNode[3];
|
|
Index_t nd4i = elemToNode[4];
|
|
Index_t nd5i = elemToNode[5];
|
|
Index_t nd6i = elemToNode[6];
|
|
Index_t nd7i = elemToNode[7];
|
|
|
|
elemX[0] = domain.c_x(nd0i);
|
|
elemX[1] = domain.c_x(nd1i);
|
|
elemX[2] = domain.c_x(nd2i);
|
|
elemX[3] = domain.c_x(nd3i);
|
|
elemX[4] = domain.c_x(nd4i);
|
|
elemX[5] = domain.c_x(nd5i);
|
|
elemX[6] = domain.c_x(nd6i);
|
|
elemX[7] = domain.c_x(nd7i);
|
|
|
|
elemY[0] = domain.c_y(nd0i);
|
|
elemY[1] = domain.c_y(nd1i);
|
|
elemY[2] = domain.c_y(nd2i);
|
|
elemY[3] = domain.c_y(nd3i);
|
|
elemY[4] = domain.c_y(nd4i);
|
|
elemY[5] = domain.c_y(nd5i);
|
|
elemY[6] = domain.c_y(nd6i);
|
|
elemY[7] = domain.c_y(nd7i);
|
|
|
|
elemZ[0] = domain.c_z(nd0i);
|
|
elemZ[1] = domain.c_z(nd1i);
|
|
elemZ[2] = domain.c_z(nd2i);
|
|
elemZ[3] = domain.c_z(nd3i);
|
|
elemZ[4] = domain.c_z(nd4i);
|
|
elemZ[5] = domain.c_z(nd5i);
|
|
elemZ[6] = domain.c_z(nd6i);
|
|
elemZ[7] = domain.c_z(nd7i);
|
|
}
|
|
|
|
static inline void
|
|
InitStressTermsForElems(Domain& domain, Real_t* sigxx, Real_t* sigyy, Real_t* sigzz,
|
|
Index_t numElem)
|
|
{
|
|
Kokkos::parallel_for(
|
|
"InitStressTermsForElems", numElem, KOKKOS_LAMBDA(const Index_t& i) {
|
|
sigxx[i] = sigyy[i] = sigzz[i] = -domain.p(i) - domain.q(i);
|
|
});
|
|
}
|
|
|
|
KOKKOS_INLINE_FUNCTION void
|
|
CalcElemShapeFunctionDerivatives(Real_t const x[], Real_t const y[], Real_t const z[],
|
|
Real_t b[][8], Real_t* const volume)
|
|
{
|
|
const Real_t x0 = x[0];
|
|
const Real_t x1 = x[1];
|
|
const Real_t x2 = x[2];
|
|
const Real_t x3 = x[3];
|
|
const Real_t x4 = x[4];
|
|
const Real_t x5 = x[5];
|
|
const Real_t x6 = x[6];
|
|
const Real_t x7 = x[7];
|
|
|
|
const Real_t y0 = y[0];
|
|
const Real_t y1 = y[1];
|
|
const Real_t y2 = y[2];
|
|
const Real_t y3 = y[3];
|
|
const Real_t y4 = y[4];
|
|
const Real_t y5 = y[5];
|
|
const Real_t y6 = y[6];
|
|
const Real_t y7 = y[7];
|
|
|
|
const Real_t z0 = z[0];
|
|
const Real_t z1 = z[1];
|
|
const Real_t z2 = z[2];
|
|
const Real_t z3 = z[3];
|
|
const Real_t z4 = z[4];
|
|
const Real_t z5 = z[5];
|
|
const Real_t z6 = z[6];
|
|
const Real_t z7 = z[7];
|
|
|
|
Real_t fjxxi, fjxet, fjxze;
|
|
Real_t fjyxi, fjyet, fjyze;
|
|
Real_t fjzxi, fjzet, fjzze;
|
|
Real_t cjxxi, cjxet, cjxze;
|
|
Real_t cjyxi, cjyet, cjyze;
|
|
Real_t cjzxi, cjzet, cjzze;
|
|
|
|
fjxxi = Real_t(.125) * ((x6 - x0) + (x5 - x3) - (x7 - x1) - (x4 - x2));
|
|
fjxet = Real_t(.125) * ((x6 - x0) - (x5 - x3) + (x7 - x1) - (x4 - x2));
|
|
fjxze = Real_t(.125) * ((x6 - x0) + (x5 - x3) + (x7 - x1) + (x4 - x2));
|
|
|
|
fjyxi = Real_t(.125) * ((y6 - y0) + (y5 - y3) - (y7 - y1) - (y4 - y2));
|
|
fjyet = Real_t(.125) * ((y6 - y0) - (y5 - y3) + (y7 - y1) - (y4 - y2));
|
|
fjyze = Real_t(.125) * ((y6 - y0) + (y5 - y3) + (y7 - y1) + (y4 - y2));
|
|
|
|
fjzxi = Real_t(.125) * ((z6 - z0) + (z5 - z3) - (z7 - z1) - (z4 - z2));
|
|
fjzet = Real_t(.125) * ((z6 - z0) - (z5 - z3) + (z7 - z1) - (z4 - z2));
|
|
fjzze = Real_t(.125) * ((z6 - z0) + (z5 - z3) + (z7 - z1) + (z4 - z2));
|
|
|
|
cjxxi = (fjyet * fjzze) - (fjzet * fjyze);
|
|
cjxet = -(fjyxi * fjzze) + (fjzxi * fjyze);
|
|
cjxze = (fjyxi * fjzet) - (fjzxi * fjyet);
|
|
|
|
cjyxi = -(fjxet * fjzze) + (fjzet * fjxze);
|
|
cjyet = (fjxxi * fjzze) - (fjzxi * fjxze);
|
|
cjyze = -(fjxxi * fjzet) + (fjzxi * fjxet);
|
|
|
|
cjzxi = (fjxet * fjyze) - (fjyet * fjxze);
|
|
cjzet = -(fjxxi * fjyze) + (fjyxi * fjxze);
|
|
cjzze = (fjxxi * fjyet) - (fjyxi * fjxet);
|
|
|
|
b[0][0] = -cjxxi - cjxet - cjxze;
|
|
b[0][1] = cjxxi - cjxet - cjxze;
|
|
b[0][2] = cjxxi + cjxet - cjxze;
|
|
b[0][3] = -cjxxi + cjxet - cjxze;
|
|
b[0][4] = -b[0][2];
|
|
b[0][5] = -b[0][3];
|
|
b[0][6] = -b[0][0];
|
|
b[0][7] = -b[0][1];
|
|
|
|
b[1][0] = -cjyxi - cjyet - cjyze;
|
|
b[1][1] = cjyxi - cjyet - cjyze;
|
|
b[1][2] = cjyxi + cjyet - cjyze;
|
|
b[1][3] = -cjyxi + cjyet - cjyze;
|
|
b[1][4] = -b[1][2];
|
|
b[1][5] = -b[1][3];
|
|
b[1][6] = -b[1][0];
|
|
b[1][7] = -b[1][1];
|
|
|
|
b[2][0] = -cjzxi - cjzet - cjzze;
|
|
b[2][1] = cjzxi - cjzet - cjzze;
|
|
b[2][2] = cjzxi + cjzet - cjzze;
|
|
b[2][3] = -cjzxi + cjzet - cjzze;
|
|
b[2][4] = -b[2][2];
|
|
b[2][5] = -b[2][3];
|
|
b[2][6] = -b[2][0];
|
|
b[2][7] = -b[2][1];
|
|
|
|
*volume = Real_t(8.) * (fjxet * cjxet + fjyet * cjyet + fjzet * cjzet);
|
|
}
|
|
|
|
KOKKOS_INLINE_FUNCTION void
|
|
SumElemFaceNormal(Real_t* normalX0, Real_t* normalY0, Real_t* normalZ0, Real_t* normalX1,
|
|
Real_t* normalY1, Real_t* normalZ1, Real_t* normalX2, Real_t* normalY2,
|
|
Real_t* normalZ2, Real_t* normalX3, Real_t* normalY3, Real_t* normalZ3,
|
|
const Real_t x0, const Real_t y0, const Real_t z0, const Real_t x1,
|
|
const Real_t y1, const Real_t z1, const Real_t x2, const Real_t y2,
|
|
const Real_t z2, const Real_t x3, const Real_t y3, const Real_t z3)
|
|
{
|
|
Real_t bisectX0 = Real_t(0.5) * (x3 + x2 - x1 - x0);
|
|
Real_t bisectY0 = Real_t(0.5) * (y3 + y2 - y1 - y0);
|
|
Real_t bisectZ0 = Real_t(0.5) * (z3 + z2 - z1 - z0);
|
|
Real_t bisectX1 = Real_t(0.5) * (x2 + x1 - x3 - x0);
|
|
Real_t bisectY1 = Real_t(0.5) * (y2 + y1 - y3 - y0);
|
|
Real_t bisectZ1 = Real_t(0.5) * (z2 + z1 - z3 - z0);
|
|
Real_t areaX = Real_t(0.25) * (bisectY0 * bisectZ1 - bisectZ0 * bisectY1);
|
|
Real_t areaY = Real_t(0.25) * (bisectZ0 * bisectX1 - bisectX0 * bisectZ1);
|
|
Real_t areaZ = Real_t(0.25) * (bisectX0 * bisectY1 - bisectY0 * bisectX1);
|
|
|
|
*normalX0 += areaX;
|
|
*normalX1 += areaX;
|
|
*normalX2 += areaX;
|
|
*normalX3 += areaX;
|
|
|
|
*normalY0 += areaY;
|
|
*normalY1 += areaY;
|
|
*normalY2 += areaY;
|
|
*normalY3 += areaY;
|
|
|
|
*normalZ0 += areaZ;
|
|
*normalZ1 += areaZ;
|
|
*normalZ2 += areaZ;
|
|
*normalZ3 += areaZ;
|
|
}
|
|
|
|
KOKKOS_INLINE_FUNCTION void
|
|
CalcElemNodeNormals(Real_t pfx[8], Real_t pfy[8], Real_t pfz[8], const Real_t x[8],
|
|
const Real_t y[8], const Real_t z[8])
|
|
{
|
|
for(Index_t i = 0; i < 8; ++i)
|
|
{
|
|
pfx[i] = Real_t(0.0);
|
|
pfy[i] = Real_t(0.0);
|
|
pfz[i] = Real_t(0.0);
|
|
}
|
|
|
|
SumElemFaceNormal(&pfx[0], &pfy[0], &pfz[0], &pfx[1], &pfy[1], &pfz[1], &pfx[2],
|
|
&pfy[2], &pfz[2], &pfx[3], &pfy[3], &pfz[3], x[0], y[0], z[0], x[1],
|
|
y[1], z[1], x[2], y[2], z[2], x[3], y[3], z[3]);
|
|
|
|
SumElemFaceNormal(&pfx[0], &pfy[0], &pfz[0], &pfx[4], &pfy[4], &pfz[4], &pfx[5],
|
|
&pfy[5], &pfz[5], &pfx[1], &pfy[1], &pfz[1], x[0], y[0], z[0], x[4],
|
|
y[4], z[4], x[5], y[5], z[5], x[1], y[1], z[1]);
|
|
|
|
SumElemFaceNormal(&pfx[1], &pfy[1], &pfz[1], &pfx[5], &pfy[5], &pfz[5], &pfx[6],
|
|
&pfy[6], &pfz[6], &pfx[2], &pfy[2], &pfz[2], x[1], y[1], z[1], x[5],
|
|
y[5], z[5], x[6], y[6], z[6], x[2], y[2], z[2]);
|
|
|
|
SumElemFaceNormal(&pfx[2], &pfy[2], &pfz[2], &pfx[6], &pfy[6], &pfz[6], &pfx[7],
|
|
&pfy[7], &pfz[7], &pfx[3], &pfy[3], &pfz[3], x[2], y[2], z[2], x[6],
|
|
y[6], z[6], x[7], y[7], z[7], x[3], y[3], z[3]);
|
|
|
|
SumElemFaceNormal(&pfx[3], &pfy[3], &pfz[3], &pfx[7], &pfy[7], &pfz[7], &pfx[4],
|
|
&pfy[4], &pfz[4], &pfx[0], &pfy[0], &pfz[0], x[3], y[3], z[3], x[7],
|
|
y[7], z[7], x[4], y[4], z[4], x[0], y[0], z[0]);
|
|
|
|
SumElemFaceNormal(&pfx[4], &pfy[4], &pfz[4], &pfx[7], &pfy[7], &pfz[7], &pfx[6],
|
|
&pfy[6], &pfz[6], &pfx[5], &pfy[5], &pfz[5], x[4], y[4], z[4], x[7],
|
|
y[7], z[7], x[6], y[6], z[6], x[5], y[5], z[5]);
|
|
}
|
|
|
|
KOKKOS_INLINE_FUNCTION void
|
|
SumElemStressesToNodeForces(const Real_t B[][8], const Real_t stress_xx,
|
|
const Real_t stress_yy, const Real_t stress_zz, Real_t fx[],
|
|
Real_t fy[], Real_t fz[])
|
|
{
|
|
for(Index_t i = 0; i < 8; i++)
|
|
{
|
|
fx[i] = -(stress_xx * B[0][i]);
|
|
fy[i] = -(stress_yy * B[1][i]);
|
|
fz[i] = -(stress_zz * B[2][i]);
|
|
}
|
|
}
|
|
|
|
static inline void
|
|
IntegrateStressForElems(Domain& domain, Real_t* sigxx, Real_t* sigyy, Real_t* sigzz,
|
|
Real_t* determ, Index_t numElem, Index_t numNode)
|
|
{
|
|
Index_t numElem8 = numElem * 8;
|
|
ResizeBuffer((numElem8 * sizeof(Real_t) + 4096) * 3);
|
|
Real_t* fx_elem = AllocateFromBuffer<Real_t>(numElem8);
|
|
Real_t* fy_elem = AllocateFromBuffer<Real_t>(numElem8);
|
|
Real_t* fz_elem = AllocateFromBuffer<Real_t>(numElem8);
|
|
|
|
Kokkos::parallel_for(
|
|
"IntegrateStressForElems A", numElem, KOKKOS_LAMBDA(const int k) {
|
|
const Index_t* const elemToNode = &domain.nodelist(k, 0);
|
|
Real_t B[3][8];
|
|
Real_t x_local[8];
|
|
Real_t y_local[8];
|
|
Real_t z_local[8];
|
|
|
|
CollectDomainNodesToElemNodes(domain, elemToNode, x_local, y_local, z_local);
|
|
|
|
CalcElemShapeFunctionDerivatives(x_local, y_local, z_local, B, &determ[k]);
|
|
|
|
CalcElemNodeNormals(B[0], B[1], B[2], x_local, y_local, z_local);
|
|
|
|
SumElemStressesToNodeForces(B, sigxx[k], sigyy[k], sigzz[k], &fx_elem[k * 8],
|
|
&fy_elem[k * 8], &fz_elem[k * 8]);
|
|
});
|
|
|
|
int team_size = 1;
|
|
if(Kokkos::DefaultExecutionSpace().concurrency() > 1024) team_size = 128;
|
|
|
|
Kokkos::parallel_for(
|
|
"IntegrateStressForElems B",
|
|
Kokkos::TeamPolicy<>((numNode + 127) / 128, team_size, 2),
|
|
KOKKOS_LAMBDA(const typename Kokkos::TeamPolicy<>::member_type& team) {
|
|
const Index_t gnode_begin = team.league_rank() * 128;
|
|
const Index_t gnode_end =
|
|
(gnode_begin + 128 < numNode) ? gnode_begin + 128 : numNode;
|
|
|
|
Kokkos::parallel_for(
|
|
Kokkos::TeamThreadRange(team, gnode_begin, gnode_end),
|
|
[&](const Index_t& gnode) {
|
|
Index_t count = domain.nodeElemCount(gnode);
|
|
Index_t* cornerList = domain.nodeElemCornerList(gnode);
|
|
reduce_double3 f_tmp;
|
|
Kokkos::parallel_reduce(
|
|
Kokkos::ThreadVectorRange(team, count),
|
|
[&](const Index_t& i,
|
|
reduce_double3& tmp) { // vectorized with ivdep
|
|
Index_t elem = cornerList[i];
|
|
tmp.x += fx_elem[elem];
|
|
tmp.y += fy_elem[elem];
|
|
tmp.z += fz_elem[elem];
|
|
},
|
|
f_tmp);
|
|
Kokkos::single(Kokkos::PerThread(team), [&]() {
|
|
domain.fx(gnode) += f_tmp.x;
|
|
domain.fy(gnode) += f_tmp.y;
|
|
domain.fz(gnode) += f_tmp.z;
|
|
});
|
|
});
|
|
});
|
|
}
|
|
|
|
KOKKOS_INLINE_FUNCTION void
|
|
VoluDer(const Real_t x0, const Real_t x1, const Real_t x2, const Real_t x3,
|
|
const Real_t x4, const Real_t x5, const Real_t y0, const Real_t y1,
|
|
const Real_t y2, const Real_t y3, const Real_t y4, const Real_t y5,
|
|
const Real_t z0, const Real_t z1, const Real_t z2, const Real_t z3,
|
|
const Real_t z4, const Real_t z5, Real_t& dvdx, Real_t& dvdy, Real_t& dvdz)
|
|
{
|
|
const Real_t twelfth = Real_t(1.0) / Real_t(12.0);
|
|
|
|
dvdx = (y1 + y2) * (z0 + z1) - (y0 + y1) * (z1 + z2) + (y0 + y4) * (z3 + z4) -
|
|
(y3 + y4) * (z0 + z4) - (y2 + y5) * (z3 + z5) + (y3 + y5) * (z2 + z5);
|
|
dvdy = -(x1 + x2) * (z0 + z1) + (x0 + x1) * (z1 + z2) - (x0 + x4) * (z3 + z4) +
|
|
(x3 + x4) * (z0 + z4) + (x2 + x5) * (z3 + z5) - (x3 + x5) * (z2 + z5);
|
|
|
|
dvdz = -(y1 + y2) * (x0 + x1) + (y0 + y1) * (x1 + x2) - (y0 + y4) * (x3 + x4) +
|
|
(y3 + y4) * (x0 + x4) + (y2 + y5) * (x3 + x5) - (y3 + y5) * (x2 + x5);
|
|
|
|
dvdx *= twelfth;
|
|
dvdy *= twelfth;
|
|
dvdz *= twelfth;
|
|
}
|
|
|
|
KOKKOS_INLINE_FUNCTION
|
|
void
|
|
CalcElemVolumeDerivative(
|
|
const Int_t& i,
|
|
const Kokkos::View<Real_t**, Kokkos::MemoryTraits<Kokkos::Unmanaged>>& dvdx,
|
|
const Kokkos::View<Real_t**, Kokkos::MemoryTraits<Kokkos::Unmanaged>>& dvdy,
|
|
const Kokkos::View<Real_t**, Kokkos::MemoryTraits<Kokkos::Unmanaged>>& dvdz,
|
|
const Real_t x[8], const Real_t y[8], const Real_t z[8])
|
|
{
|
|
#pragma nounroll
|
|
for(int j = 0; j < 4; j++)
|
|
{
|
|
VoluDer(x[(j + 1) % 4], x[(j + 2) % 4], x[(j + 3) % 4], x[(j + 0) % 4 + 4],
|
|
x[(j + 1) % 4 + 4], x[(j + 3) % 4 + 4], y[(j + 1) % 4], y[(j + 2) % 4],
|
|
y[(j + 3) % 4], y[(j + 0) % 4 + 4], y[(j + 1) % 4 + 4],
|
|
y[(j + 3) % 4 + 4], z[(j + 1) % 4], z[(j + 2) % 4], z[(j + 3) % 4],
|
|
z[(j + 0) % 4 + 4], z[(j + 1) % 4 + 4], z[(j + 3) % 4 + 4], dvdx(i, j),
|
|
dvdy(i, j), dvdz(i, j));
|
|
VoluDer(x[(j + 3) % 4 + 4], x[(j + 2) % 4 + 4], x[(j + 1) % 4 + 4],
|
|
x[(j + 0) % 4], x[(j + 3) % 4], x[(j + 1) % 4], y[(j + 3) % 4 + 4],
|
|
y[(j + 2) % 4 + 4], y[(j + 1) % 4 + 4], y[(j + 0) % 4], y[(j + 3) % 4],
|
|
y[(j + 1) % 4], z[(j + 3) % 4 + 4], z[(j + 2) % 4 + 4],
|
|
z[(j + 1) % 4 + 4], z[(j + 0) % 4], z[(j + 3) % 4], z[(j + 1) % 4],
|
|
dvdx(i, j + 4), dvdy(i, j + 4), dvdz(i, j + 4));
|
|
}
|
|
}
|
|
|
|
KOKKOS_INLINE_FUNCTION
|
|
void
|
|
CalcElemFBHourglassForce(const Real_t* xd, const Real_t hourgam[][8],
|
|
const Real_t& coefficient, Real_t* hgfx)
|
|
{
|
|
Real_t hxx[4];
|
|
for(Index_t i = 0; i < 4; i++)
|
|
{
|
|
hxx[i] = hourgam[i][0] * xd[0] + hourgam[i][1] * xd[1] + hourgam[i][2] * xd[2] +
|
|
hourgam[i][3] * xd[3] + hourgam[i][4] * xd[4] + hourgam[i][5] * xd[5] +
|
|
hourgam[i][6] * xd[6] + hourgam[i][7] * xd[7];
|
|
}
|
|
for(Index_t i = 0; i < 8; i++)
|
|
{
|
|
hgfx[i] = coefficient * (hourgam[0][i] * hxx[0] + hourgam[1][i] * hxx[1] +
|
|
hourgam[2][i] * hxx[2] + hourgam[3][i] * hxx[3]);
|
|
}
|
|
}
|
|
|
|
struct Gamma
|
|
{
|
|
Real_t gamma[4][8];
|
|
Gamma()
|
|
{
|
|
gamma[0][0] = Real_t(1.);
|
|
gamma[0][1] = Real_t(1.);
|
|
gamma[0][2] = Real_t(-1.);
|
|
gamma[0][3] = Real_t(-1.);
|
|
gamma[0][4] = Real_t(-1.);
|
|
gamma[0][5] = Real_t(-1.);
|
|
gamma[0][6] = Real_t(1.);
|
|
gamma[0][7] = Real_t(1.);
|
|
gamma[1][0] = Real_t(1.);
|
|
gamma[1][1] = Real_t(-1.);
|
|
gamma[1][2] = Real_t(-1.);
|
|
gamma[1][3] = Real_t(1.);
|
|
gamma[1][4] = Real_t(-1.);
|
|
gamma[1][5] = Real_t(1.);
|
|
gamma[1][6] = Real_t(1.);
|
|
gamma[1][7] = Real_t(-1.);
|
|
gamma[2][0] = Real_t(1.);
|
|
gamma[2][1] = Real_t(-1.);
|
|
gamma[2][2] = Real_t(1.);
|
|
gamma[2][3] = Real_t(-1.);
|
|
gamma[2][4] = Real_t(1.);
|
|
gamma[2][5] = Real_t(-1.);
|
|
gamma[2][6] = Real_t(1.);
|
|
gamma[2][7] = Real_t(-1.);
|
|
gamma[3][0] = Real_t(-1.);
|
|
gamma[3][1] = Real_t(1.);
|
|
gamma[3][2] = Real_t(-1.);
|
|
gamma[3][3] = Real_t(1.);
|
|
gamma[3][4] = Real_t(1.);
|
|
gamma[3][5] = Real_t(-1.);
|
|
gamma[3][6] = Real_t(1.);
|
|
gamma[3][7] = Real_t(-1.);
|
|
}
|
|
};
|
|
|
|
static inline void
|
|
CalcFBHourglassForceForElems(
|
|
Domain& domain, Real_t* determ,
|
|
const Kokkos::View<const Real_t**, Kokkos::MemoryTraits<Kokkos::Unmanaged>> x8n,
|
|
const Kokkos::View<const Real_t**, Kokkos::MemoryTraits<Kokkos::Unmanaged>> y8n,
|
|
const Kokkos::View<const Real_t**, Kokkos::MemoryTraits<Kokkos::Unmanaged>> z8n,
|
|
const Kokkos::View<const Real_t**, Kokkos::MemoryTraits<Kokkos::Unmanaged>> dvdx,
|
|
const Kokkos::View<const Real_t**, Kokkos::MemoryTraits<Kokkos::Unmanaged>> dvdy,
|
|
const Kokkos::View<const Real_t**, Kokkos::MemoryTraits<Kokkos::Unmanaged>> dvdz,
|
|
Real_t hourg, Index_t numElem, Index_t numNode)
|
|
{
|
|
Index_t numElem8 = numElem * 8;
|
|
|
|
Real_t* fx_elem;
|
|
Real_t* fy_elem;
|
|
Real_t* fz_elem;
|
|
|
|
if(do_atomic == 0)
|
|
{
|
|
fx_elem = AllocateFromBuffer<Real_t>(numElem8);
|
|
fy_elem = AllocateFromBuffer<Real_t>(numElem8);
|
|
fz_elem = AllocateFromBuffer<Real_t>(numElem8);
|
|
}
|
|
|
|
Gamma G;
|
|
|
|
Int_t do_atomic_dev = do_atomic;
|
|
|
|
Kokkos::parallel_for(
|
|
"CalcFBHourglassForceForElems A", numElem, KOKKOS_LAMBDA(const int& i2) {
|
|
Real_t *fx_local, *fy_local, *fz_local;
|
|
Real_t hgfx[8];
|
|
|
|
Real_t hourgam[4][8];
|
|
Real_t xd1[8];
|
|
|
|
const Index_t* elemToNode = &domain.nodelist(i2, 0);
|
|
Index_t i3 = 8 * i2;
|
|
Real_t volinv = Real_t(1.0) / determ[i2];
|
|
|
|
for(Index_t i1 = 0; i1 < 4; ++i1)
|
|
{
|
|
Real_t hourmodx = 0.0;
|
|
for(int j = 0; j < 8; j++)
|
|
hourmodx += x8n(i2, j) * G.gamma[i1][j];
|
|
|
|
Real_t hourmody = 0.0;
|
|
for(int j = 0; j < 8; j++)
|
|
hourmody += y8n(i2, j) * G.gamma[i1][j];
|
|
|
|
Real_t hourmodz = 0.0;
|
|
for(int j = 0; j < 8; j++)
|
|
hourmodz += z8n(i2, j) * G.gamma[i1][j];
|
|
|
|
#pragma ivdep
|
|
for(int j = 0; j < 8; j++)
|
|
hourgam[i1][j] = G.gamma[i1][j] - volinv * (dvdx(i2, j) * hourmodx +
|
|
dvdy(i2, j) * hourmody +
|
|
dvdz(i2, j) * hourmodz);
|
|
}
|
|
|
|
const Real_t ss1 = domain.ss(i2);
|
|
const Real_t mass1 = domain.elemMass(i2);
|
|
const Real_t volume13 = CBRT(determ[i2]);
|
|
|
|
const Index_t n0si2 = elemToNode[0];
|
|
const Index_t n1si2 = elemToNode[1];
|
|
const Index_t n2si2 = elemToNode[2];
|
|
const Index_t n3si2 = elemToNode[3];
|
|
const Index_t n4si2 = elemToNode[4];
|
|
const Index_t n5si2 = elemToNode[5];
|
|
const Index_t n6si2 = elemToNode[6];
|
|
const Index_t n7si2 = elemToNode[7];
|
|
|
|
const Real_t coefficient = -hourg * Real_t(0.01) * ss1 * mass1 / volume13;
|
|
|
|
xd1[0] = domain.xd(n0si2);
|
|
xd1[1] = domain.xd(n1si2);
|
|
xd1[2] = domain.xd(n2si2);
|
|
xd1[3] = domain.xd(n3si2);
|
|
xd1[4] = domain.xd(n4si2);
|
|
xd1[5] = domain.xd(n5si2);
|
|
xd1[6] = domain.xd(n6si2);
|
|
xd1[7] = domain.xd(n7si2);
|
|
|
|
CalcElemFBHourglassForce(xd1, hourgam, coefficient, hgfx);
|
|
|
|
if(!do_atomic_dev)
|
|
{
|
|
fx_local = &fx_elem[i3];
|
|
fx_local[0] = hgfx[0];
|
|
fx_local[1] = hgfx[1];
|
|
fx_local[2] = hgfx[2];
|
|
fx_local[3] = hgfx[3];
|
|
fx_local[4] = hgfx[4];
|
|
fx_local[5] = hgfx[5];
|
|
fx_local[6] = hgfx[6];
|
|
fx_local[7] = hgfx[7];
|
|
}
|
|
else
|
|
{
|
|
Kokkos::atomic_add(&domain.fx(n0si2), hgfx[0]);
|
|
Kokkos::atomic_add(&domain.fx(n1si2), hgfx[1]);
|
|
Kokkos::atomic_add(&domain.fx(n2si2), hgfx[2]);
|
|
Kokkos::atomic_add(&domain.fx(n3si2), hgfx[3]);
|
|
Kokkos::atomic_add(&domain.fx(n4si2), hgfx[4]);
|
|
Kokkos::atomic_add(&domain.fx(n5si2), hgfx[5]);
|
|
Kokkos::atomic_add(&domain.fx(n6si2), hgfx[6]);
|
|
Kokkos::atomic_add(&domain.fx(n7si2), hgfx[7]);
|
|
}
|
|
|
|
xd1[0] = domain.yd(n0si2);
|
|
xd1[1] = domain.yd(n1si2);
|
|
xd1[2] = domain.yd(n2si2);
|
|
xd1[3] = domain.yd(n3si2);
|
|
xd1[4] = domain.yd(n4si2);
|
|
xd1[5] = domain.yd(n5si2);
|
|
xd1[6] = domain.yd(n6si2);
|
|
xd1[7] = domain.yd(n7si2);
|
|
|
|
CalcElemFBHourglassForce(xd1, hourgam, coefficient, hgfx);
|
|
|
|
if(!do_atomic_dev)
|
|
{
|
|
fy_local = &fy_elem[i3];
|
|
fy_local[0] = hgfx[0];
|
|
fy_local[1] = hgfx[1];
|
|
fy_local[2] = hgfx[2];
|
|
fy_local[3] = hgfx[3];
|
|
fy_local[4] = hgfx[4];
|
|
fy_local[5] = hgfx[5];
|
|
fy_local[6] = hgfx[6];
|
|
fy_local[7] = hgfx[7];
|
|
}
|
|
else
|
|
{
|
|
Kokkos::atomic_add(&domain.fy(n0si2), hgfx[0]);
|
|
Kokkos::atomic_add(&domain.fy(n1si2), hgfx[1]);
|
|
Kokkos::atomic_add(&domain.fy(n2si2), hgfx[2]);
|
|
Kokkos::atomic_add(&domain.fy(n3si2), hgfx[3]);
|
|
Kokkos::atomic_add(&domain.fy(n4si2), hgfx[4]);
|
|
Kokkos::atomic_add(&domain.fy(n5si2), hgfx[5]);
|
|
Kokkos::atomic_add(&domain.fy(n6si2), hgfx[6]);
|
|
Kokkos::atomic_add(&domain.fy(n7si2), hgfx[7]);
|
|
}
|
|
|
|
xd1[0] = domain.zd(n0si2);
|
|
xd1[1] = domain.zd(n1si2);
|
|
xd1[2] = domain.zd(n2si2);
|
|
xd1[3] = domain.zd(n3si2);
|
|
xd1[4] = domain.zd(n4si2);
|
|
xd1[5] = domain.zd(n5si2);
|
|
xd1[6] = domain.zd(n6si2);
|
|
xd1[7] = domain.zd(n7si2);
|
|
|
|
CalcElemFBHourglassForce(xd1, hourgam, coefficient, hgfx);
|
|
|
|
if(!do_atomic_dev)
|
|
{
|
|
fz_local = &fz_elem[i3];
|
|
fz_local[0] = hgfx[0];
|
|
fz_local[1] = hgfx[1];
|
|
fz_local[2] = hgfx[2];
|
|
fz_local[3] = hgfx[3];
|
|
fz_local[4] = hgfx[4];
|
|
fz_local[5] = hgfx[5];
|
|
fz_local[6] = hgfx[6];
|
|
fz_local[7] = hgfx[7];
|
|
}
|
|
else
|
|
{
|
|
Kokkos::atomic_add(&domain.fz(n0si2), hgfx[0]);
|
|
Kokkos::atomic_add(&domain.fz(n1si2), hgfx[1]);
|
|
Kokkos::atomic_add(&domain.fz(n2si2), hgfx[2]);
|
|
Kokkos::atomic_add(&domain.fz(n3si2), hgfx[3]);
|
|
Kokkos::atomic_add(&domain.fz(n4si2), hgfx[4]);
|
|
Kokkos::atomic_add(&domain.fz(n5si2), hgfx[5]);
|
|
Kokkos::atomic_add(&domain.fz(n6si2), hgfx[6]);
|
|
Kokkos::atomic_add(&domain.fz(n7si2), hgfx[7]);
|
|
}
|
|
});
|
|
|
|
if(!do_atomic)
|
|
{
|
|
int team_size = 1;
|
|
if(Kokkos::DefaultExecutionSpace().concurrency() > 1024) team_size = 128;
|
|
|
|
Kokkos::parallel_for(
|
|
"CalcFBHourglassForceForElems B",
|
|
Kokkos::TeamPolicy<>((numNode + 127) / 128, team_size, 2),
|
|
KOKKOS_LAMBDA(const typename Kokkos::TeamPolicy<>::member_type& team) {
|
|
const Index_t gnode_begin = team.league_rank() * 128;
|
|
const Index_t gnode_end =
|
|
(gnode_begin + 128 < numNode) ? gnode_begin + 128 : numNode;
|
|
|
|
Kokkos::parallel_for(
|
|
Kokkos::TeamThreadRange(team, gnode_begin, gnode_end),
|
|
[&](const Index_t& gnode) {
|
|
Index_t count = domain.nodeElemCount(gnode);
|
|
Index_t* cornerList = domain.nodeElemCornerList(gnode);
|
|
reduce_double3 f_tmp;
|
|
Kokkos::parallel_reduce(
|
|
Kokkos::ThreadVectorRange(team, count),
|
|
[&](const Index_t& i,
|
|
reduce_double3& tmp) { // vectorized with ivdep
|
|
Index_t elem = cornerList[i];
|
|
tmp.x += fx_elem[elem];
|
|
tmp.y += fy_elem[elem];
|
|
tmp.z += fz_elem[elem];
|
|
},
|
|
f_tmp);
|
|
Kokkos::single(Kokkos::PerThread(team), [&]() {
|
|
domain.fx(gnode) += f_tmp.x;
|
|
domain.fy(gnode) += f_tmp.y;
|
|
domain.fz(gnode) += f_tmp.z;
|
|
});
|
|
});
|
|
});
|
|
}
|
|
}
|
|
|
|
static inline void
|
|
CalcHourglassControlForElems(Domain& domain, Real_t determ[], Real_t hgcoef)
|
|
{
|
|
Index_t numElem = domain.numElem();
|
|
Index_t numElem8 = numElem * 8;
|
|
ResizeBuffer((numElem8 * sizeof(Real_t) + 4096) * (do_atomic ? 6 : 9));
|
|
|
|
Real_t* dvdx = AllocateFromBuffer<Real_t>(numElem8);
|
|
Real_t* dvdy = AllocateFromBuffer<Real_t>(numElem8);
|
|
Real_t* dvdz = AllocateFromBuffer<Real_t>(numElem8);
|
|
Real_t* x8n = AllocateFromBuffer<Real_t>(numElem8);
|
|
Real_t* y8n = AllocateFromBuffer<Real_t>(numElem8);
|
|
Real_t* z8n = AllocateFromBuffer<Real_t>(numElem8);
|
|
Kokkos::View<Real_t**, Kokkos::MemoryTraits<Kokkos::Unmanaged>> v_x8n(x8n, numElem,
|
|
8);
|
|
Kokkos::View<Real_t**, Kokkos::MemoryTraits<Kokkos::Unmanaged>> v_y8n(y8n, numElem,
|
|
8);
|
|
Kokkos::View<Real_t**, Kokkos::MemoryTraits<Kokkos::Unmanaged>> v_z8n(z8n, numElem,
|
|
8);
|
|
Kokkos::View<Real_t**, Kokkos::MemoryTraits<Kokkos::Unmanaged>> v_dvdx(dvdx, numElem,
|
|
8);
|
|
Kokkos::View<Real_t**, Kokkos::MemoryTraits<Kokkos::Unmanaged>> v_dvdy(dvdy, numElem,
|
|
8);
|
|
Kokkos::View<Real_t**, Kokkos::MemoryTraits<Kokkos::Unmanaged>> v_dvdz(dvdz, numElem,
|
|
8);
|
|
|
|
int error = 0;
|
|
Kokkos::parallel_reduce(
|
|
"CalcHourglassControlForElems", numElem,
|
|
KOKKOS_LAMBDA(const int i, int& err) {
|
|
Real_t x1[8], y1[8], z1[8];
|
|
|
|
Index_t* elemToNode = &domain.nodelist(i, 0);
|
|
CollectDomainNodesToElemNodes(domain, elemToNode, x1, y1, z1);
|
|
|
|
CalcElemVolumeDerivative(i, v_dvdx, v_dvdy, v_dvdz, x1, y1, z1);
|
|
|
|
for(Index_t ii = 0; ii < 8; ++ii)
|
|
{
|
|
v_x8n(i, ii) = x1[ii];
|
|
v_y8n(i, ii) = y1[ii];
|
|
v_z8n(i, ii) = z1[ii];
|
|
}
|
|
|
|
determ[i] = domain.volo(i) * domain.v(i);
|
|
|
|
if(domain.v(i) <= Real_t(0.0))
|
|
{
|
|
err++;
|
|
}
|
|
},
|
|
error);
|
|
|
|
if(error)
|
|
#if USE_MPI
|
|
MPI_Abort(MPI_COMM_WORLD, VolumeError);
|
|
#else
|
|
exit(VolumeError);
|
|
#endif
|
|
|
|
if(hgcoef > Real_t(0.))
|
|
{
|
|
CalcFBHourglassForceForElems(domain, determ, v_x8n, v_y8n, v_z8n, v_dvdx, v_dvdy,
|
|
v_dvdz, hgcoef, numElem, domain.numNode());
|
|
}
|
|
|
|
return;
|
|
}
|
|
|
|
static inline void
|
|
CalcVolumeForceForElems(Domain& domain)
|
|
{
|
|
Index_t numElem = domain.numElem();
|
|
if(numElem != 0)
|
|
{
|
|
Real_t hgcoef = domain.hgcoef();
|
|
Kokkos::View<Real_t*> sigxx("sigxx", numElem);
|
|
Kokkos::View<Real_t*> sigyy("sigyy", numElem);
|
|
Kokkos::View<Real_t*> sigzz("sigzz", numElem);
|
|
Kokkos::View<Real_t*> determ("determ", numElem);
|
|
|
|
InitStressTermsForElems(domain, sigxx.data(), sigyy.data(), sigzz.data(),
|
|
numElem);
|
|
|
|
IntegrateStressForElems(domain, sigxx.data(), sigyy.data(), sigzz.data(),
|
|
determ.data(), numElem, domain.numNode());
|
|
|
|
// check for negative element volume
|
|
int error = 0;
|
|
Kokkos::parallel_reduce(
|
|
"CalcVolumeForceForElems", numElem,
|
|
KOKKOS_LAMBDA(const int k, int& err) {
|
|
if(determ[k] <= Real_t(0.0))
|
|
{
|
|
err++;
|
|
}
|
|
},
|
|
error);
|
|
|
|
if(error)
|
|
#if USE_MPI
|
|
MPI_Abort(MPI_COMM_WORLD, VolumeError);
|
|
#else
|
|
exit(VolumeError);
|
|
#endif
|
|
|
|
CalcHourglassControlForElems(domain, determ.data(), hgcoef);
|
|
}
|
|
}
|
|
|
|
static inline void
|
|
CalcForceForNodes(Domain& domain)
|
|
{
|
|
Index_t numNode = domain.numNode();
|
|
|
|
#if USE_MPI
|
|
CommRecv(domain, MSG_COMM_SBN, 3, domain.sizeX() + 1, domain.sizeY() + 1,
|
|
domain.sizeZ() + 1, true, false);
|
|
#endif
|
|
|
|
Kokkos::parallel_for(
|
|
"CalcForceForNodes", numNode, KOKKOS_LAMBDA(const int i) {
|
|
domain.fx(i) = Real_t(0.0);
|
|
domain.fy(i) = Real_t(0.0);
|
|
domain.fz(i) = Real_t(0.0);
|
|
});
|
|
|
|
CalcVolumeForceForElems(domain);
|
|
|
|
#if USE_MPI
|
|
Domain_member fieldData[3];
|
|
fieldData[0] = &Domain::fx;
|
|
fieldData[1] = &Domain::fy;
|
|
fieldData[2] = &Domain::fz;
|
|
|
|
CommSend(domain, MSG_COMM_SBN, 3, fieldData, domain.sizeX() + 1, domain.sizeY() + 1,
|
|
domain.sizeZ() + 1, true, false);
|
|
CommSBN(domain, 3, fieldData);
|
|
#endif
|
|
}
|
|
|
|
static inline void
|
|
CalcAccelerationForNodes(Domain& domain, Index_t numNode)
|
|
{
|
|
Kokkos::parallel_for(
|
|
"CalcAccelerationForNodes", numNode, KOKKOS_LAMBDA(const int i) {
|
|
domain.xdd(i) = domain.fx(i) / domain.nodalMass(i);
|
|
domain.ydd(i) = domain.fy(i) / domain.nodalMass(i);
|
|
domain.zdd(i) = domain.fz(i) / domain.nodalMass(i);
|
|
});
|
|
}
|
|
|
|
static inline void
|
|
ApplyAccelerationBoundaryConditionsForNodes(Domain& domain)
|
|
{
|
|
Index_t size = domain.sizeX();
|
|
Index_t numNodeBC = (size + 1) * (size + 1);
|
|
|
|
if(!domain.symmXempty() != 0)
|
|
{
|
|
Kokkos::parallel_for(
|
|
"ApplyAccelerationBoundaryConditionsForNodes A", numNodeBC,
|
|
KOKKOS_LAMBDA(const int i) { domain.xdd(domain.symmX(i)) = Real_t(0.0); });
|
|
}
|
|
|
|
if(!domain.symmYempty() != 0)
|
|
{
|
|
Kokkos::parallel_for(
|
|
"ApplyAccelerationBoundaryConditionsForNodes B", numNodeBC,
|
|
KOKKOS_LAMBDA(const int i) { domain.ydd(domain.symmY(i)) = Real_t(0.0); });
|
|
}
|
|
|
|
if(!domain.symmZempty() != 0)
|
|
{
|
|
Kokkos::parallel_for(
|
|
"ApplyAccelerationBoundaryConditionsForNodes C", numNodeBC,
|
|
KOKKOS_LAMBDA(const int i) { domain.zdd(domain.symmZ(i)) = Real_t(0.0); });
|
|
}
|
|
}
|
|
|
|
static inline void
|
|
CalcVelocityForNodes(Domain& domain, const Real_t dt, const Real_t u_cut, Index_t numNode)
|
|
{
|
|
Kokkos::parallel_for(
|
|
"CalcVelocityForNodes", numNode, KOKKOS_LAMBDA(const int i) {
|
|
Real_t xdtmp, ydtmp, zdtmp;
|
|
|
|
xdtmp = domain.xd(i) + domain.xdd(i) * dt;
|
|
if(FABS(xdtmp) < u_cut) xdtmp = Real_t(0.0);
|
|
domain.xd(i) = xdtmp;
|
|
|
|
ydtmp = domain.yd(i) + domain.ydd(i) * dt;
|
|
if(FABS(ydtmp) < u_cut) ydtmp = Real_t(0.0);
|
|
domain.yd(i) = ydtmp;
|
|
|
|
zdtmp = domain.zd(i) + domain.zdd(i) * dt;
|
|
if(FABS(zdtmp) < u_cut) zdtmp = Real_t(0.0);
|
|
domain.zd(i) = zdtmp;
|
|
});
|
|
}
|
|
|
|
static inline void
|
|
CalcPositionForNodes(Domain& domain, const Real_t dt, Index_t numNode)
|
|
{
|
|
Kokkos::parallel_for(
|
|
"CalcPositionForNodes", numNode, KOKKOS_LAMBDA(const int i) {
|
|
domain.x(i) += domain.xd(i) * dt;
|
|
domain.y(i) += domain.yd(i) * dt;
|
|
domain.z(i) += domain.zd(i) * dt;
|
|
});
|
|
}
|
|
|
|
static inline void
|
|
LagrangeNodal(Domain& domain)
|
|
{
|
|
#ifdef SEDOV_SYNC_POS_VEL_EARLY
|
|
Domain_member fieldData[6];
|
|
#endif
|
|
|
|
const Real_t delt = domain.deltatime();
|
|
Real_t u_cut = domain.u_cut();
|
|
|
|
CalcForceForNodes(domain);
|
|
|
|
#if USE_MPI
|
|
# ifdef SEDOV_SYNC_POS_VEL_EARLY
|
|
CommRecv(domain, MSG_SYNC_POS_VEL, 6, domain.sizeX() + 1, domain.sizeY() + 1,
|
|
domain.sizeZ() + 1, false, false);
|
|
# endif
|
|
#endif
|
|
|
|
CalcAccelerationForNodes(domain, domain.numNode());
|
|
|
|
ApplyAccelerationBoundaryConditionsForNodes(domain);
|
|
|
|
CalcVelocityForNodes(domain, delt, u_cut, domain.numNode());
|
|
|
|
CalcPositionForNodes(domain, delt, domain.numNode());
|
|
#if USE_MPI
|
|
# ifdef SEDOV_SYNC_POS_VEL_EARLY
|
|
fieldData[0] = &Domain::x;
|
|
fieldData[1] = &Domain::y;
|
|
fieldData[2] = &Domain::z;
|
|
fieldData[3] = &Domain::xd;
|
|
fieldData[4] = &Domain::yd;
|
|
fieldData[5] = &Domain::zd;
|
|
|
|
CommSend(domain, MSG_SYNC_POS_VEL, 6, fieldData, domain.sizeX() + 1,
|
|
domain.sizeY() + 1, domain.sizeZ() + 1, false, false);
|
|
CommSyncPosVel(domain);
|
|
# endif
|
|
#endif
|
|
|
|
return;
|
|
}
|
|
|
|
KOKKOS_INLINE_FUNCTION Real_t
|
|
CalcElemVolume(const Real_t x0, const Real_t x1, const Real_t x2, const Real_t x3,
|
|
const Real_t x4, const Real_t x5, const Real_t x6, const Real_t x7,
|
|
const Real_t y0, const Real_t y1, const Real_t y2, const Real_t y3,
|
|
const Real_t y4, const Real_t y5, const Real_t y6, const Real_t y7,
|
|
const Real_t z0, const Real_t z1, const Real_t z2, const Real_t z3,
|
|
const Real_t z4, const Real_t z5, const Real_t z6, const Real_t z7)
|
|
{
|
|
Real_t twelveth = Real_t(1.0) / Real_t(12.0);
|
|
|
|
Real_t dx61 = x6 - x1;
|
|
Real_t dy61 = y6 - y1;
|
|
Real_t dz61 = z6 - z1;
|
|
|
|
Real_t dx70 = x7 - x0;
|
|
Real_t dy70 = y7 - y0;
|
|
Real_t dz70 = z7 - z0;
|
|
|
|
Real_t dx63 = x6 - x3;
|
|
Real_t dy63 = y6 - y3;
|
|
Real_t dz63 = z6 - z3;
|
|
|
|
Real_t dx20 = x2 - x0;
|
|
Real_t dy20 = y2 - y0;
|
|
Real_t dz20 = z2 - z0;
|
|
|
|
Real_t dx50 = x5 - x0;
|
|
Real_t dy50 = y5 - y0;
|
|
Real_t dz50 = z5 - z0;
|
|
|
|
Real_t dx64 = x6 - x4;
|
|
Real_t dy64 = y6 - y4;
|
|
Real_t dz64 = z6 - z4;
|
|
|
|
Real_t dx31 = x3 - x1;
|
|
Real_t dy31 = y3 - y1;
|
|
Real_t dz31 = z3 - z1;
|
|
|
|
Real_t dx72 = x7 - x2;
|
|
Real_t dy72 = y7 - y2;
|
|
Real_t dz72 = z7 - z2;
|
|
|
|
Real_t dx43 = x4 - x3;
|
|
Real_t dy43 = y4 - y3;
|
|
Real_t dz43 = z4 - z3;
|
|
|
|
Real_t dx57 = x5 - x7;
|
|
Real_t dy57 = y5 - y7;
|
|
Real_t dz57 = z5 - z7;
|
|
|
|
Real_t dx14 = x1 - x4;
|
|
Real_t dy14 = y1 - y4;
|
|
Real_t dz14 = z1 - z4;
|
|
|
|
Real_t dx25 = x2 - x5;
|
|
Real_t dy25 = y2 - y5;
|
|
Real_t dz25 = z2 - z5;
|
|
|
|
#define TRIPLE_PRODUCT(x1, y1, z1, x2, y2, z2, x3, y3, z3) \
|
|
((x1) * ((y2) * (z3) - (z2) * (y3)) + (x2) * ((z1) * (y3) - (y1) * (z3)) + \
|
|
(x3) * ((y1) * (z2) - (z1) * (y2)))
|
|
|
|
Real_t volume = TRIPLE_PRODUCT(dx31 + dx72, dx63, dx20, dy31 + dy72, dy63, dy20,
|
|
dz31 + dz72, dz63, dz20) +
|
|
TRIPLE_PRODUCT(dx43 + dx57, dx64, dx70, dy43 + dy57, dy64, dy70,
|
|
dz43 + dz57, dz64, dz70) +
|
|
TRIPLE_PRODUCT(dx14 + dx25, dx61, dx50, dy14 + dy25, dy61, dy50,
|
|
dz14 + dz25, dz61, dz50);
|
|
|
|
#undef TRIPLE_PRODUCT
|
|
|
|
volume *= twelveth;
|
|
|
|
return volume;
|
|
}
|
|
|
|
KOKKOS_INLINE_FUNCTION
|
|
Real_t
|
|
CalcElemVolume(const Real_t x[8], const Real_t y[8], const Real_t z[8])
|
|
{
|
|
return CalcElemVolume(x[0], x[1], x[2], x[3], x[4], x[5], x[6], x[7], y[0], y[1],
|
|
y[2], y[3], y[4], y[5], y[6], y[7], z[0], z[1], z[2], z[3],
|
|
z[4], z[5], z[6], z[7]);
|
|
}
|
|
|
|
KOKKOS_INLINE_FUNCTION
|
|
Real_t
|
|
AreaFace(const Real_t x0, const Real_t x1, const Real_t x2, const Real_t x3,
|
|
const Real_t y0, const Real_t y1, const Real_t y2, const Real_t y3,
|
|
const Real_t z0, const Real_t z1, const Real_t z2, const Real_t z3)
|
|
{
|
|
Real_t fx = (x2 - x0) - (x3 - x1);
|
|
Real_t fy = (y2 - y0) - (y3 - y1);
|
|
Real_t fz = (z2 - z0) - (z3 - z1);
|
|
Real_t gx = (x2 - x0) + (x3 - x1);
|
|
Real_t gy = (y2 - y0) + (y3 - y1);
|
|
Real_t gz = (z2 - z0) + (z3 - z1);
|
|
Real_t area = (fx * fx + fy * fy + fz * fz) * (gx * gx + gy * gy + gz * gz) -
|
|
(fx * gx + fy * gy + fz * gz) * (fx * gx + fy * gy + fz * gz);
|
|
return area;
|
|
}
|
|
|
|
KOKKOS_INLINE_FUNCTION Real_t
|
|
CalcElemCharacteristicLength(const Real_t x[8], const Real_t y[8], const Real_t z[8],
|
|
const Real_t volume)
|
|
{
|
|
Real_t a, charLength = Real_t(0.0);
|
|
|
|
a = AreaFace(x[0], x[1], x[2], x[3], y[0], y[1], y[2], y[3], z[0], z[1], z[2], z[3]);
|
|
charLength = MAX(a, charLength);
|
|
|
|
a = AreaFace(x[4], x[5], x[6], x[7], y[4], y[5], y[6], y[7], z[4], z[5], z[6], z[7]);
|
|
charLength = MAX(a, charLength);
|
|
|
|
a = AreaFace(x[0], x[1], x[5], x[4], y[0], y[1], y[5], y[4], z[0], z[1], z[5], z[4]);
|
|
charLength = MAX(a, charLength);
|
|
|
|
a = AreaFace(x[1], x[2], x[6], x[5], y[1], y[2], y[6], y[5], z[1], z[2], z[6], z[5]);
|
|
charLength = MAX(a, charLength);
|
|
|
|
a = AreaFace(x[2], x[3], x[7], x[6], y[2], y[3], y[7], y[6], z[2], z[3], z[7], z[6]);
|
|
charLength = MAX(a, charLength);
|
|
|
|
a = AreaFace(x[3], x[0], x[4], x[7], y[3], y[0], y[4], y[7], z[3], z[0], z[4], z[7]);
|
|
charLength = MAX(a, charLength);
|
|
|
|
charLength = Real_t(4.0) * volume / SQRT(charLength);
|
|
|
|
return charLength;
|
|
}
|
|
|
|
KOKKOS_INLINE_FUNCTION void
|
|
CalcElemVelocityGradient(const Real_t* const xvel, const Real_t* const yvel,
|
|
const Real_t* const zvel, const Real_t b[][8], const Real_t detJ,
|
|
Real_t* const d)
|
|
{
|
|
const Real_t inv_detJ = Real_t(1.0) / detJ;
|
|
Real_t dyddx, dxddy, dzddx, dxddz, dzddy, dyddz;
|
|
const Real_t* const pfx = b[0];
|
|
const Real_t* const pfy = b[1];
|
|
const Real_t* const pfz = b[2];
|
|
|
|
d[0] = inv_detJ * (pfx[0] * (xvel[0] - xvel[6]) + pfx[1] * (xvel[1] - xvel[7]) +
|
|
pfx[2] * (xvel[2] - xvel[4]) + pfx[3] * (xvel[3] - xvel[5]));
|
|
|
|
d[1] = inv_detJ * (pfy[0] * (yvel[0] - yvel[6]) + pfy[1] * (yvel[1] - yvel[7]) +
|
|
pfy[2] * (yvel[2] - yvel[4]) + pfy[3] * (yvel[3] - yvel[5]));
|
|
|
|
d[2] = inv_detJ * (pfz[0] * (zvel[0] - zvel[6]) + pfz[1] * (zvel[1] - zvel[7]) +
|
|
pfz[2] * (zvel[2] - zvel[4]) + pfz[3] * (zvel[3] - zvel[5]));
|
|
|
|
dyddx = inv_detJ * (pfx[0] * (yvel[0] - yvel[6]) + pfx[1] * (yvel[1] - yvel[7]) +
|
|
pfx[2] * (yvel[2] - yvel[4]) + pfx[3] * (yvel[3] - yvel[5]));
|
|
|
|
dxddy = inv_detJ * (pfy[0] * (xvel[0] - xvel[6]) + pfy[1] * (xvel[1] - xvel[7]) +
|
|
pfy[2] * (xvel[2] - xvel[4]) + pfy[3] * (xvel[3] - xvel[5]));
|
|
|
|
dzddx = inv_detJ * (pfx[0] * (zvel[0] - zvel[6]) + pfx[1] * (zvel[1] - zvel[7]) +
|
|
pfx[2] * (zvel[2] - zvel[4]) + pfx[3] * (zvel[3] - zvel[5]));
|
|
|
|
dxddz = inv_detJ * (pfz[0] * (xvel[0] - xvel[6]) + pfz[1] * (xvel[1] - xvel[7]) +
|
|
pfz[2] * (xvel[2] - xvel[4]) + pfz[3] * (xvel[3] - xvel[5]));
|
|
|
|
dzddy = inv_detJ * (pfy[0] * (zvel[0] - zvel[6]) + pfy[1] * (zvel[1] - zvel[7]) +
|
|
pfy[2] * (zvel[2] - zvel[4]) + pfy[3] * (zvel[3] - zvel[5]));
|
|
|
|
dyddz = inv_detJ * (pfz[0] * (yvel[0] - yvel[6]) + pfz[1] * (yvel[1] - yvel[7]) +
|
|
pfz[2] * (yvel[2] - yvel[4]) + pfz[3] * (yvel[3] - yvel[5]));
|
|
d[5] = Real_t(.5) * (dxddy + dyddx);
|
|
d[4] = Real_t(.5) * (dxddz + dzddx);
|
|
d[3] = Real_t(.5) * (dzddy + dyddz);
|
|
}
|
|
|
|
void
|
|
CalcKinematicsForElems(Domain& domain, Real_t deltaTime, Index_t numElem)
|
|
{
|
|
Kokkos::parallel_for(
|
|
"CalcKinematicsForElems", numElem, KOKKOS_LAMBDA(const int k) {
|
|
Real_t B[3][8];
|
|
Real_t D[6];
|
|
Real_t x_local[8];
|
|
Real_t y_local[8];
|
|
Real_t z_local[8];
|
|
Real_t xd_local[8];
|
|
Real_t yd_local[8];
|
|
Real_t zd_local[8];
|
|
Real_t detJ = Real_t(0.0);
|
|
|
|
Real_t volume;
|
|
Real_t relativeVolume;
|
|
const Index_t* const elemToNode = &domain.nodelist(k, 0);
|
|
|
|
CollectDomainNodesToElemNodes(domain, elemToNode, x_local, y_local, z_local);
|
|
|
|
volume = CalcElemVolume(x_local, y_local, z_local);
|
|
relativeVolume = volume / domain.volo(k);
|
|
domain.vnew(k) = relativeVolume;
|
|
domain.delv(k) = relativeVolume - domain.v(k);
|
|
|
|
domain.arealg(k) =
|
|
CalcElemCharacteristicLength(x_local, y_local, z_local, volume);
|
|
|
|
for(Index_t lnode = 0; lnode < 8; ++lnode)
|
|
{
|
|
Index_t gnode = elemToNode[lnode];
|
|
xd_local[lnode] = domain.c_xd(gnode);
|
|
yd_local[lnode] = domain.c_yd(gnode);
|
|
zd_local[lnode] = domain.c_zd(gnode);
|
|
}
|
|
|
|
Real_t dt2 = Real_t(0.5) * deltaTime;
|
|
for(Index_t j = 0; j < 8; ++j)
|
|
{
|
|
x_local[j] -= dt2 * xd_local[j];
|
|
y_local[j] -= dt2 * yd_local[j];
|
|
z_local[j] -= dt2 * zd_local[j];
|
|
}
|
|
|
|
CalcElemShapeFunctionDerivatives(x_local, y_local, z_local, B, &detJ);
|
|
|
|
CalcElemVelocityGradient(xd_local, yd_local, zd_local, B, detJ, D);
|
|
|
|
domain.dxx(k) = D[0];
|
|
domain.dyy(k) = D[1];
|
|
domain.dzz(k) = D[2];
|
|
});
|
|
}
|
|
|
|
static inline void
|
|
CalcLagrangeElements(Domain& domain)
|
|
{
|
|
Index_t numElem = domain.numElem();
|
|
if(numElem > 0)
|
|
{
|
|
const Real_t deltatime = domain.deltatime();
|
|
|
|
domain.AllocateStrains(numElem);
|
|
|
|
CalcKinematicsForElems(domain, deltatime, numElem);
|
|
|
|
int error = 0;
|
|
Kokkos::parallel_reduce(
|
|
"CalcLagrangeElements", numElem,
|
|
KOKKOS_LAMBDA(const int k, int& err) {
|
|
Real_t vdov = domain.dxx(k) + domain.dyy(k) + domain.dzz(k);
|
|
Real_t vdovthird = vdov / Real_t(3.0);
|
|
|
|
domain.vdov(k) = vdov;
|
|
domain.dxx(k) -= vdovthird;
|
|
domain.dyy(k) -= vdovthird;
|
|
domain.dzz(k) -= vdovthird;
|
|
|
|
if(domain.vnew(k) <= Real_t(0.0))
|
|
{
|
|
err++;
|
|
}
|
|
},
|
|
error);
|
|
|
|
if(error)
|
|
#if USE_MPI
|
|
MPI_Abort(MPI_COMM_WORLD, VolumeError);
|
|
#else
|
|
exit(VolumeError);
|
|
#endif
|
|
|
|
domain.DeallocateStrains();
|
|
}
|
|
}
|
|
|
|
static inline void
|
|
CalcMonotonicQGradientsForElems(Domain& domain)
|
|
{
|
|
Index_t numElem = domain.numElem();
|
|
|
|
Kokkos::parallel_for(
|
|
"CalcMonotonicQGradientsForElems", numElem, KOKKOS_LAMBDA(const int i) {
|
|
const Real_t ptiny = Real_t(1.e-36);
|
|
Real_t ax, ay, az;
|
|
Real_t dxv, dyv, dzv;
|
|
|
|
const Index_t* elemToNode = &domain.nodelist(i, 0);
|
|
Index_t n0 = elemToNode[0];
|
|
Index_t n1 = elemToNode[1];
|
|
Index_t n2 = elemToNode[2];
|
|
Index_t n3 = elemToNode[3];
|
|
Index_t n4 = elemToNode[4];
|
|
Index_t n5 = elemToNode[5];
|
|
Index_t n6 = elemToNode[6];
|
|
Index_t n7 = elemToNode[7];
|
|
|
|
Real_t x0 = domain.x(n0);
|
|
Real_t x1 = domain.x(n1);
|
|
Real_t x2 = domain.x(n2);
|
|
Real_t x3 = domain.x(n3);
|
|
Real_t x4 = domain.x(n4);
|
|
Real_t x5 = domain.x(n5);
|
|
Real_t x6 = domain.x(n6);
|
|
Real_t x7 = domain.x(n7);
|
|
|
|
Real_t y0 = domain.y(n0);
|
|
Real_t y1 = domain.y(n1);
|
|
Real_t y2 = domain.y(n2);
|
|
Real_t y3 = domain.y(n3);
|
|
Real_t y4 = domain.y(n4);
|
|
Real_t y5 = domain.y(n5);
|
|
Real_t y6 = domain.y(n6);
|
|
Real_t y7 = domain.y(n7);
|
|
|
|
Real_t z0 = domain.z(n0);
|
|
Real_t z1 = domain.z(n1);
|
|
Real_t z2 = domain.z(n2);
|
|
Real_t z3 = domain.z(n3);
|
|
Real_t z4 = domain.z(n4);
|
|
Real_t z5 = domain.z(n5);
|
|
Real_t z6 = domain.z(n6);
|
|
Real_t z7 = domain.z(n7);
|
|
|
|
Real_t xv0 = domain.xd(n0);
|
|
Real_t xv1 = domain.xd(n1);
|
|
Real_t xv2 = domain.xd(n2);
|
|
Real_t xv3 = domain.xd(n3);
|
|
Real_t xv4 = domain.xd(n4);
|
|
Real_t xv5 = domain.xd(n5);
|
|
Real_t xv6 = domain.xd(n6);
|
|
Real_t xv7 = domain.xd(n7);
|
|
|
|
Real_t yv0 = domain.yd(n0);
|
|
Real_t yv1 = domain.yd(n1);
|
|
Real_t yv2 = domain.yd(n2);
|
|
Real_t yv3 = domain.yd(n3);
|
|
Real_t yv4 = domain.yd(n4);
|
|
Real_t yv5 = domain.yd(n5);
|
|
Real_t yv6 = domain.yd(n6);
|
|
Real_t yv7 = domain.yd(n7);
|
|
|
|
Real_t zv0 = domain.zd(n0);
|
|
Real_t zv1 = domain.zd(n1);
|
|
Real_t zv2 = domain.zd(n2);
|
|
Real_t zv3 = domain.zd(n3);
|
|
Real_t zv4 = domain.zd(n4);
|
|
Real_t zv5 = domain.zd(n5);
|
|
Real_t zv6 = domain.zd(n6);
|
|
Real_t zv7 = domain.zd(n7);
|
|
|
|
Real_t vol = domain.volo(i) * domain.vnew(i);
|
|
Real_t norm = Real_t(1.0) / (vol + ptiny);
|
|
|
|
Real_t dxj = Real_t(-0.25) * ((x0 + x1 + x5 + x4) - (x3 + x2 + x6 + x7));
|
|
Real_t dyj = Real_t(-0.25) * ((y0 + y1 + y5 + y4) - (y3 + y2 + y6 + y7));
|
|
Real_t dzj = Real_t(-0.25) * ((z0 + z1 + z5 + z4) - (z3 + z2 + z6 + z7));
|
|
|
|
Real_t dxi = Real_t(0.25) * ((x1 + x2 + x6 + x5) - (x0 + x3 + x7 + x4));
|
|
Real_t dyi = Real_t(0.25) * ((y1 + y2 + y6 + y5) - (y0 + y3 + y7 + y4));
|
|
Real_t dzi = Real_t(0.25) * ((z1 + z2 + z6 + z5) - (z0 + z3 + z7 + z4));
|
|
|
|
Real_t dxk = Real_t(0.25) * ((x4 + x5 + x6 + x7) - (x0 + x1 + x2 + x3));
|
|
Real_t dyk = Real_t(0.25) * ((y4 + y5 + y6 + y7) - (y0 + y1 + y2 + y3));
|
|
Real_t dzk = Real_t(0.25) * ((z4 + z5 + z6 + z7) - (z0 + z1 + z2 + z3));
|
|
|
|
ax = dyi * dzj - dzi * dyj;
|
|
ay = dzi * dxj - dxi * dzj;
|
|
az = dxi * dyj - dyi * dxj;
|
|
|
|
domain.delx_zeta(i) = vol / SQRT(ax * ax + ay * ay + az * az + ptiny);
|
|
|
|
ax *= norm;
|
|
ay *= norm;
|
|
az *= norm;
|
|
|
|
dxv = Real_t(0.25) * ((xv4 + xv5 + xv6 + xv7) - (xv0 + xv1 + xv2 + xv3));
|
|
dyv = Real_t(0.25) * ((yv4 + yv5 + yv6 + yv7) - (yv0 + yv1 + yv2 + yv3));
|
|
dzv = Real_t(0.25) * ((zv4 + zv5 + zv6 + zv7) - (zv0 + zv1 + zv2 + zv3));
|
|
|
|
domain.delv_zeta(i) = ax * dxv + ay * dyv + az * dzv;
|
|
|
|
ax = dyj * dzk - dzj * dyk;
|
|
ay = dzj * dxk - dxj * dzk;
|
|
az = dxj * dyk - dyj * dxk;
|
|
|
|
domain.delx_xi(i) = vol / SQRT(ax * ax + ay * ay + az * az + ptiny);
|
|
|
|
ax *= norm;
|
|
ay *= norm;
|
|
az *= norm;
|
|
|
|
dxv = Real_t(0.25) * ((xv1 + xv2 + xv6 + xv5) - (xv0 + xv3 + xv7 + xv4));
|
|
dyv = Real_t(0.25) * ((yv1 + yv2 + yv6 + yv5) - (yv0 + yv3 + yv7 + yv4));
|
|
dzv = Real_t(0.25) * ((zv1 + zv2 + zv6 + zv5) - (zv0 + zv3 + zv7 + zv4));
|
|
|
|
domain.delv_xi(i) = ax * dxv + ay * dyv + az * dzv;
|
|
|
|
ax = dyk * dzi - dzk * dyi;
|
|
ay = dzk * dxi - dxk * dzi;
|
|
az = dxk * dyi - dyk * dxi;
|
|
|
|
domain.delx_eta(i) = vol / SQRT(ax * ax + ay * ay + az * az + ptiny);
|
|
|
|
ax *= norm;
|
|
ay *= norm;
|
|
az *= norm;
|
|
|
|
dxv = Real_t(-0.25) * ((xv0 + xv1 + xv5 + xv4) - (xv3 + xv2 + xv6 + xv7));
|
|
dyv = Real_t(-0.25) * ((yv0 + yv1 + yv5 + yv4) - (yv3 + yv2 + yv6 + yv7));
|
|
dzv = Real_t(-0.25) * ((zv0 + zv1 + zv5 + zv4) - (zv3 + zv2 + zv6 + zv7));
|
|
|
|
domain.delv_eta(i) = ax * dxv + ay * dyv + az * dzv;
|
|
});
|
|
}
|
|
|
|
static inline void
|
|
CalcMonotonicQRegionForElems(Domain& domain, Int_t r, Real_t ptiny)
|
|
{
|
|
Real_t monoq_limiter_mult = domain.monoq_limiter_mult();
|
|
Real_t monoq_max_slope = domain.monoq_max_slope();
|
|
Real_t qlc_monoq = domain.qlc_monoq();
|
|
Real_t qqc_monoq = domain.qqc_monoq();
|
|
|
|
Kokkos::parallel_for(
|
|
"CalcMonotonicQRegionForElems", domain.regElemSize(r),
|
|
KOKKOS_LAMBDA(const int i) {
|
|
Index_t ielem = domain.regElemlist(r, i);
|
|
Real_t qlin, qquad;
|
|
Real_t phixi, phieta, phizeta;
|
|
Int_t bcMask = domain.elemBC(ielem);
|
|
Real_t delvm = 0.0, delvp = 0.0;
|
|
|
|
Real_t norm = Real_t(1.) / (domain.delv_xi(ielem) + ptiny);
|
|
|
|
switch(bcMask & XI_M)
|
|
{
|
|
case XI_M_COMM:
|
|
case 0: delvm = domain.delv_xi(domain.lxim(ielem)); break;
|
|
case XI_M_SYMM: delvm = domain.delv_xi(ielem); break;
|
|
case XI_M_FREE: delvm = Real_t(0.0); break;
|
|
default:
|
|
printf("Error in switch at %s line %d\n", __FILE__, __LINE__);
|
|
delvm = 0;
|
|
break;
|
|
}
|
|
switch(bcMask & XI_P)
|
|
{
|
|
case XI_P_COMM:
|
|
case 0: delvp = domain.delv_xi(domain.lxip(ielem)); break;
|
|
case XI_P_SYMM: delvp = domain.delv_xi(ielem); break;
|
|
case XI_P_FREE: delvp = Real_t(0.0); break;
|
|
default:
|
|
printf("Error in switch at %s line %d\n", __FILE__, __LINE__);
|
|
delvp = 0;
|
|
break;
|
|
}
|
|
|
|
delvm = delvm * norm;
|
|
delvp = delvp * norm;
|
|
|
|
phixi = Real_t(.5) * (delvm + delvp);
|
|
|
|
delvm *= monoq_limiter_mult;
|
|
delvp *= monoq_limiter_mult;
|
|
|
|
if(delvm < phixi) phixi = delvm;
|
|
if(delvp < phixi) phixi = delvp;
|
|
if(phixi < Real_t(0.)) phixi = Real_t(0.);
|
|
if(phixi > monoq_max_slope) phixi = monoq_max_slope;
|
|
|
|
norm = Real_t(1.) / (domain.delv_eta(ielem) + ptiny);
|
|
|
|
switch(bcMask & ETA_M)
|
|
{
|
|
case ETA_M_COMM:
|
|
case 0: delvm = domain.delv_eta(domain.letam(ielem)); break;
|
|
case ETA_M_SYMM: delvm = domain.delv_eta(ielem); break;
|
|
case ETA_M_FREE: delvm = Real_t(0.0); break;
|
|
default:
|
|
printf("Error in switch at %s line %d\n", __FILE__, __LINE__);
|
|
delvm = 0;
|
|
break;
|
|
}
|
|
switch(bcMask & ETA_P)
|
|
{
|
|
case ETA_P_COMM:
|
|
case 0: delvp = domain.delv_eta(domain.letap(ielem)); break;
|
|
case ETA_P_SYMM: delvp = domain.delv_eta(ielem); break;
|
|
case ETA_P_FREE: delvp = Real_t(0.0); break;
|
|
default:
|
|
printf("Error in switch at %s line %d\n", __FILE__, __LINE__);
|
|
delvp = 0;
|
|
break;
|
|
}
|
|
|
|
delvm = delvm * norm;
|
|
delvp = delvp * norm;
|
|
|
|
phieta = Real_t(.5) * (delvm + delvp);
|
|
|
|
delvm *= monoq_limiter_mult;
|
|
delvp *= monoq_limiter_mult;
|
|
|
|
if(delvm < phieta) phieta = delvm;
|
|
if(delvp < phieta) phieta = delvp;
|
|
if(phieta < Real_t(0.)) phieta = Real_t(0.);
|
|
if(phieta > monoq_max_slope) phieta = monoq_max_slope;
|
|
|
|
norm = Real_t(1.) / (domain.delv_zeta(ielem) + ptiny);
|
|
|
|
switch(bcMask & ZETA_M)
|
|
{
|
|
case ZETA_M_COMM:
|
|
case 0: delvm = domain.delv_zeta(domain.lzetam(ielem)); break;
|
|
case ZETA_M_SYMM: delvm = domain.delv_zeta(ielem); break;
|
|
case ZETA_M_FREE: delvm = Real_t(0.0); break;
|
|
default:
|
|
printf("Error in switch at %s line %d\n", __FILE__, __LINE__);
|
|
delvm = 0;
|
|
break;
|
|
}
|
|
switch(bcMask & ZETA_P)
|
|
{
|
|
case ZETA_P_COMM:
|
|
case 0: delvp = domain.delv_zeta(domain.lzetap(ielem)); break;
|
|
case ZETA_P_SYMM: delvp = domain.delv_zeta(ielem); break;
|
|
case ZETA_P_FREE: delvp = Real_t(0.0); break;
|
|
default:
|
|
printf("Error in switch at %s line %d\n", __FILE__, __LINE__);
|
|
delvp = 0;
|
|
break;
|
|
}
|
|
|
|
delvm = delvm * norm;
|
|
delvp = delvp * norm;
|
|
|
|
phizeta = Real_t(.5) * (delvm + delvp);
|
|
|
|
delvm *= monoq_limiter_mult;
|
|
delvp *= monoq_limiter_mult;
|
|
|
|
if(delvm < phizeta) phizeta = delvm;
|
|
if(delvp < phizeta) phizeta = delvp;
|
|
if(phizeta < Real_t(0.)) phizeta = Real_t(0.);
|
|
if(phizeta > monoq_max_slope) phizeta = monoq_max_slope;
|
|
|
|
if(domain.vdov(ielem) > Real_t(0.))
|
|
{
|
|
qlin = Real_t(0.);
|
|
qquad = Real_t(0.);
|
|
}
|
|
else
|
|
{
|
|
Real_t delvxxi = domain.delv_xi(ielem) * domain.delx_xi(ielem);
|
|
Real_t delvxeta = domain.delv_eta(ielem) * domain.delx_eta(ielem);
|
|
Real_t delvxzeta = domain.delv_zeta(ielem) * domain.delx_zeta(ielem);
|
|
|
|
if(delvxxi > Real_t(0.)) delvxxi = Real_t(0.);
|
|
if(delvxeta > Real_t(0.)) delvxeta = Real_t(0.);
|
|
if(delvxzeta > Real_t(0.)) delvxzeta = Real_t(0.);
|
|
|
|
Real_t rho =
|
|
domain.elemMass(ielem) / (domain.volo(ielem) * domain.vnew(ielem));
|
|
|
|
qlin =
|
|
-qlc_monoq * rho *
|
|
(delvxxi * (Real_t(1.) - phixi) + delvxeta * (Real_t(1.) - phieta) +
|
|
delvxzeta * (Real_t(1.) - phizeta));
|
|
|
|
qquad = qqc_monoq * rho *
|
|
(delvxxi * delvxxi * (Real_t(1.) - phixi * phixi) +
|
|
delvxeta * delvxeta * (Real_t(1.) - phieta * phieta) +
|
|
delvxzeta * delvxzeta * (Real_t(1.) - phizeta * phizeta));
|
|
}
|
|
|
|
domain.qq(ielem) = qquad;
|
|
domain.ql(ielem) = qlin;
|
|
});
|
|
}
|
|
|
|
static inline void
|
|
CalcMonotonicQForElems(Domain& domain)
|
|
{
|
|
const Real_t ptiny = Real_t(1.e-36);
|
|
|
|
for(Index_t r = 0; r < domain.numReg(); ++r)
|
|
{
|
|
if(domain.regElemSize(r) > 0)
|
|
{
|
|
CalcMonotonicQRegionForElems(domain, r, ptiny);
|
|
}
|
|
}
|
|
}
|
|
|
|
static inline void
|
|
CalcQForElems(Domain& domain)
|
|
{
|
|
Index_t numElem = domain.numElem();
|
|
|
|
if(numElem != 0)
|
|
{
|
|
Int_t allElem = numElem + /* local elem */
|
|
2 * domain.sizeX() * domain.sizeY() + /* plane ghosts */
|
|
2 * domain.sizeX() * domain.sizeZ() + /* row ghosts */
|
|
2 * domain.sizeY() * domain.sizeZ(); /* col ghosts */
|
|
|
|
domain.AllocateGradients(numElem, allElem);
|
|
|
|
#if USE_MPI
|
|
CommRecv(domain, MSG_MONOQ, 3, domain.sizeX(), domain.sizeY(), domain.sizeZ(),
|
|
true, true);
|
|
#endif
|
|
CalcMonotonicQGradientsForElems(domain);
|
|
|
|
#if USE_MPI
|
|
Domain_member fieldData[3];
|
|
|
|
fieldData[0] = &Domain::delv_xi;
|
|
fieldData[1] = &Domain::delv_eta;
|
|
fieldData[2] = &Domain::delv_zeta;
|
|
|
|
CommSend(domain, MSG_MONOQ, 3, fieldData, domain.sizeX(), domain.sizeY(),
|
|
domain.sizeZ(), true, true);
|
|
|
|
CommMonoQ(domain);
|
|
#endif
|
|
|
|
CalcMonotonicQForElems(domain);
|
|
|
|
domain.DeallocateGradients();
|
|
|
|
Index_t idx = 0;
|
|
Kokkos::parallel_reduce(
|
|
"CalcQForElems", numElem,
|
|
KOKKOS_LAMBDA(const Index_t& i, Index_t& count) {
|
|
if(domain.q(i) > domain.qstop())
|
|
{
|
|
count++;
|
|
}
|
|
},
|
|
idx);
|
|
|
|
if(idx > 0)
|
|
{
|
|
#if USE_MPI
|
|
MPI_Abort(MPI_COMM_WORLD, QStopError);
|
|
#else
|
|
exit(QStopError);
|
|
#endif
|
|
}
|
|
}
|
|
}
|
|
|
|
KOKKOS_INLINE_FUNCTION
|
|
void
|
|
CalcPressureForElem(Real_t& p_new_i, Real_t& bvc_i, Real_t& pbvc_i, const Real_t& e_old_i,
|
|
const Real_t& compression_i, const Real_t& vnewc_e,
|
|
const Real_t& pmin, const Real_t& p_cut, const Real_t& eosvmax)
|
|
{
|
|
const Real_t c1s = Real_t(2.0) / Real_t(3.0);
|
|
bvc_i = c1s * (compression_i + Real_t(1.));
|
|
|
|
pbvc_i = c1s;
|
|
|
|
p_new_i = bvc_i * e_old_i;
|
|
|
|
if(FABS(p_new_i) < p_cut) p_new_i = Real_t(0.0);
|
|
|
|
if(vnewc_e >= eosvmax) /* impossible condition here? */
|
|
p_new_i = Real_t(0.0);
|
|
|
|
if(p_new_i < pmin) p_new_i = pmin;
|
|
}
|
|
|
|
static inline void
|
|
CalcEnergyForElems(Real_t* p_new, Real_t* e_new, Real_t* q_new, Real_t* bvc, Real_t* pbvc,
|
|
Real_t* p_old, Real_t* e_old, Real_t* q_old, Real_t* compression,
|
|
Real_t* compHalfStep, Real_t* vnewc, Real_t* work, Real_t* delvc,
|
|
Real_t pmin, Real_t p_cut, Real_t e_cut, Real_t q_cut, Real_t emin,
|
|
Real_t* qq_old, Real_t* ql_old, Real_t rho0, Real_t eosvmax,
|
|
Index_t length, Domain& domain, Index_t r)
|
|
{
|
|
Kokkos::parallel_for(
|
|
"CalcEnergyForElems", length, KOKKOS_LAMBDA(const int i) {
|
|
const Real_t delvc_i = delvc[i];
|
|
const Real_t p_old_i = p_old[i];
|
|
const Real_t q_old_i = q_old[i];
|
|
Real_t e_new_i = e_old[i] - Real_t(0.5) * delvc_i * (p_old_i + q_old_i) +
|
|
Real_t(0.5) * work[i];
|
|
|
|
if(e_new_i < emin)
|
|
{
|
|
e_new_i = emin;
|
|
}
|
|
|
|
Real_t bvc_i, pbvc_i;
|
|
Real_t pHalfStep_i;
|
|
const Real_t vnewc_e = vnewc[domain.regElemlist(r, i)];
|
|
const Real_t compHalfStep_i = compHalfStep[i];
|
|
CalcPressureForElem(pHalfStep_i, bvc_i, pbvc_i, e_new_i, compHalfStep_i,
|
|
vnewc_e, pmin, p_cut, eosvmax);
|
|
|
|
Real_t vhalf = Real_t(1.) / (Real_t(1.) + compHalfStep_i);
|
|
|
|
Real_t q_new_i;
|
|
const Real_t ql_old_i = ql_old[i];
|
|
const Real_t qq_old_i = qq_old[i];
|
|
if(delvc_i > Real_t(0.))
|
|
{
|
|
q_new_i /* = qq_old[i] = ql_old[i] */ = Real_t(0.);
|
|
}
|
|
else
|
|
{
|
|
Real_t ssc =
|
|
(pbvc_i * e_new_i + vhalf * vhalf * bvc_i * pHalfStep_i) / rho0;
|
|
|
|
if(ssc <= Real_t(.1111111e-36))
|
|
{
|
|
ssc = Real_t(.3333333e-18);
|
|
}
|
|
else
|
|
{
|
|
ssc = SQRT(ssc);
|
|
}
|
|
|
|
q_new_i = (ssc * ql_old_i + qq_old_i);
|
|
}
|
|
|
|
e_new_i = e_new_i + Real_t(0.5) * delvc_i *
|
|
(Real_t(3.0) * (p_old_i + q_old_i) -
|
|
Real_t(4.0) * (pHalfStep_i + q_new_i));
|
|
|
|
e_new_i += Real_t(0.5) * work[i];
|
|
|
|
if(FABS(e_new_i) < e_cut)
|
|
{
|
|
e_new_i = Real_t(0.);
|
|
}
|
|
if(e_new_i < emin)
|
|
{
|
|
e_new_i = emin;
|
|
}
|
|
Real_t p_new_i;
|
|
const Real_t compression_i = compression[i];
|
|
CalcPressureForElem(p_new_i, bvc_i, pbvc_i, e_new_i, compression_i, vnewc_e,
|
|
pmin, p_cut, eosvmax);
|
|
|
|
const Real_t sixth = Real_t(1.0) / Real_t(6.0);
|
|
Real_t q_tilde;
|
|
|
|
if(delvc_i > Real_t(0.))
|
|
{
|
|
q_tilde = Real_t(0.);
|
|
}
|
|
else
|
|
{
|
|
Real_t ssc =
|
|
(pbvc_i * e_new_i + vnewc_e * vnewc_e * bvc_i * p_new_i) / rho0;
|
|
|
|
if(ssc <= Real_t(.1111111e-36))
|
|
{
|
|
ssc = Real_t(.3333333e-18);
|
|
}
|
|
else
|
|
{
|
|
ssc = SQRT(ssc);
|
|
}
|
|
|
|
q_tilde = (ssc * ql_old_i + qq_old_i);
|
|
}
|
|
|
|
e_new_i =
|
|
e_new_i - (Real_t(7.0) * (p_old_i + q_old_i) -
|
|
Real_t(8.0) * (pHalfStep_i + q_new_i) + (p_new_i + q_tilde)) *
|
|
delvc_i * sixth;
|
|
|
|
if(FABS(e_new_i) < e_cut)
|
|
{
|
|
e_new_i = Real_t(0.);
|
|
}
|
|
if(e_new_i < emin)
|
|
{
|
|
e_new_i = emin;
|
|
}
|
|
|
|
CalcPressureForElem(p_new_i, bvc_i, pbvc_i, e_new_i, compression_i, vnewc_e,
|
|
pmin, p_cut, eosvmax);
|
|
bvc[i] = bvc_i;
|
|
pbvc[i] = pbvc_i;
|
|
p_new[i] = p_new_i;
|
|
|
|
if(delvc_i <= Real_t(0.))
|
|
{
|
|
Real_t ssc =
|
|
(pbvc_i * e_new_i + vnewc_e * vnewc_e * bvc_i * p_new_i) / rho0;
|
|
|
|
if(ssc <= Real_t(.1111111e-36))
|
|
{
|
|
ssc = Real_t(.3333333e-18);
|
|
}
|
|
else
|
|
{
|
|
ssc = SQRT(ssc);
|
|
}
|
|
|
|
q_new_i = (ssc * ql_old_i + qq_old_i);
|
|
|
|
if(FABS(q_new_i) < q_cut) q_new_i = Real_t(0.);
|
|
}
|
|
q_new[i] = q_new_i;
|
|
e_new[i] = e_new_i;
|
|
});
|
|
|
|
return;
|
|
}
|
|
|
|
static inline void
|
|
CalcSoundSpeedForElems(Domain& domain, Real_t* vnewc, Real_t rho0, Real_t* enewc,
|
|
Real_t* pnewc, Real_t* pbvc, Real_t* bvc, Real_t ss4o3,
|
|
Index_t len, Index_t r)
|
|
{
|
|
Kokkos::parallel_for(
|
|
"CalcSoundSpeedForElems", len, KOKKOS_LAMBDA(const int i) {
|
|
Index_t ielem = domain.regElemlist(r, i);
|
|
Real_t ssTmp =
|
|
(pbvc[i] * enewc[i] + vnewc[ielem] * vnewc[ielem] * bvc[i] * pnewc[i]) /
|
|
rho0;
|
|
if(ssTmp <= Real_t(.1111111e-36))
|
|
{
|
|
ssTmp = Real_t(.3333333e-18);
|
|
}
|
|
else
|
|
{
|
|
ssTmp = SQRT(ssTmp);
|
|
}
|
|
domain.ss(ielem) = ssTmp;
|
|
});
|
|
}
|
|
|
|
static inline void
|
|
EvalEOSForElems(Domain& domain, Real_t* vnewc, Int_t numElemReg, Index_t r, Int_t rep)
|
|
{
|
|
Real_t e_cut = domain.e_cut();
|
|
Real_t p_cut = domain.p_cut();
|
|
Real_t ss4o3 = domain.ss4o3();
|
|
Real_t q_cut = domain.q_cut();
|
|
|
|
Real_t eosvmax = domain.eosvmax();
|
|
Real_t eosvmin = domain.eosvmin();
|
|
Real_t pmin = domain.pmin();
|
|
Real_t emin = domain.emin();
|
|
Real_t rho0 = domain.refdens();
|
|
|
|
ResizeBuffer((numElemReg * sizeof(Real_t) + 4096) * 16);
|
|
|
|
Real_t* e_old = AllocateFromBuffer<Real_t>(numElemReg);
|
|
Real_t* delvc = AllocateFromBuffer<Real_t>(numElemReg);
|
|
Real_t* p_old = AllocateFromBuffer<Real_t>(numElemReg);
|
|
Real_t* q_old = AllocateFromBuffer<Real_t>(numElemReg);
|
|
Real_t* compression = AllocateFromBuffer<Real_t>(numElemReg);
|
|
Real_t* compHalfStep = AllocateFromBuffer<Real_t>(numElemReg);
|
|
Real_t* qq_old = AllocateFromBuffer<Real_t>(numElemReg);
|
|
Real_t* ql_old = AllocateFromBuffer<Real_t>(numElemReg);
|
|
Real_t* work = AllocateFromBuffer<Real_t>(numElemReg);
|
|
Real_t* p_new = AllocateFromBuffer<Real_t>(numElemReg);
|
|
Real_t* e_new = AllocateFromBuffer<Real_t>(numElemReg);
|
|
Real_t* q_new = AllocateFromBuffer<Real_t>(numElemReg);
|
|
Real_t* bvc = AllocateFromBuffer<Real_t>(numElemReg);
|
|
Real_t* pbvc = AllocateFromBuffer<Real_t>(numElemReg);
|
|
|
|
for(Int_t j = 0; j < rep; j++)
|
|
{
|
|
Kokkos::parallel_for(
|
|
"EvalEOSForElems A", numElemReg, KOKKOS_LAMBDA(const int i) {
|
|
Index_t ielem = domain.regElemlist(r, i);
|
|
e_old[i] = domain.c_e(ielem);
|
|
delvc[i] = domain.c_delv(ielem);
|
|
p_old[i] = domain.c_p(ielem);
|
|
q_old[i] = domain.c_q(ielem);
|
|
qq_old[i] = domain.c_qq(ielem);
|
|
ql_old[i] = domain.c_ql(ielem);
|
|
const Real_t vnewc_ielem = vnewc[ielem];
|
|
Real_t vchalf;
|
|
compression[i] = Real_t(1.) / vnewc_ielem - Real_t(1.);
|
|
vchalf = vnewc_ielem - delvc[i] * Real_t(.5);
|
|
compHalfStep[i] = Real_t(1.) / vchalf - Real_t(1.);
|
|
|
|
if(eosvmin != Real_t(0.))
|
|
{
|
|
if(vnewc_ielem <= eosvmin)
|
|
{ /* impossible due to calling func? */
|
|
compHalfStep[i] = compression[i];
|
|
}
|
|
}
|
|
if(eosvmax != Real_t(0.))
|
|
{
|
|
if(vnewc_ielem >= eosvmax)
|
|
{ /* impossible due to calling func? */
|
|
p_old[i] = Real_t(0.);
|
|
compression[i] = Real_t(0.);
|
|
compHalfStep[i] = Real_t(0.);
|
|
}
|
|
}
|
|
work[i] = Real_t(0.);
|
|
});
|
|
|
|
CalcEnergyForElems(p_new, e_new, q_new, bvc, pbvc, p_old, e_old, q_old,
|
|
compression, compHalfStep, vnewc, work, delvc, pmin, p_cut,
|
|
e_cut, q_cut, emin, qq_old, ql_old, rho0, eosvmax, numElemReg,
|
|
domain, r);
|
|
}
|
|
|
|
Kokkos::parallel_for(
|
|
"EvalEOSForElems F", numElemReg, KOKKOS_LAMBDA(const int i) {
|
|
Index_t ielem = domain.regElemlist(r, i);
|
|
domain.p(ielem) = p_new[i];
|
|
domain.e(ielem) = e_new[i];
|
|
domain.q(ielem) = q_new[i];
|
|
});
|
|
|
|
CalcSoundSpeedForElems(domain, vnewc, rho0, e_new, p_new, pbvc, bvc, ss4o3,
|
|
numElemReg, r);
|
|
}
|
|
|
|
static inline void
|
|
ApplyMaterialPropertiesForElems(Domain& domain)
|
|
{
|
|
Index_t numElem = domain.numElem();
|
|
|
|
if(numElem != 0)
|
|
{
|
|
Real_t eosvmin = domain.eosvmin();
|
|
Real_t eosvmax = domain.eosvmax();
|
|
Kokkos::View<Real_t*> vnewc("vnewc", numElem);
|
|
|
|
Kokkos::parallel_for(
|
|
"ApplyMaterialPropertiesForElems A", numElem,
|
|
KOKKOS_LAMBDA(const int i) { vnewc[i] = domain.vnew(i); });
|
|
|
|
if(eosvmin != Real_t(0.))
|
|
{
|
|
Kokkos::parallel_for(
|
|
"ApplyMaterialPropertiesForElems B", numElem, KOKKOS_LAMBDA(const int i) {
|
|
if(vnewc[i] < eosvmin) vnewc[i] = eosvmin;
|
|
});
|
|
}
|
|
|
|
if(eosvmax != Real_t(0.))
|
|
{
|
|
Kokkos::parallel_for(
|
|
"ApplyMaterialPropertiesForElems C", numElem, KOKKOS_LAMBDA(const int i) {
|
|
if(vnewc[i] > eosvmax) vnewc[i] = eosvmax;
|
|
});
|
|
}
|
|
|
|
int error = 0;
|
|
Kokkos::parallel_reduce(
|
|
"ApplyMaterialPropertiesForElems", numElem,
|
|
KOKKOS_LAMBDA(const int i, int& err) {
|
|
Real_t vc = domain.v(i);
|
|
if(eosvmin != Real_t(0.))
|
|
{
|
|
if(vc < eosvmin) vc = eosvmin;
|
|
}
|
|
if(eosvmax != Real_t(0.))
|
|
{
|
|
if(vc > eosvmax) vc = eosvmax;
|
|
}
|
|
if(vc <= 0.)
|
|
{
|
|
err++;
|
|
}
|
|
},
|
|
error);
|
|
|
|
if(error != 0)
|
|
#if USE_MPI
|
|
MPI_Abort(MPI_COMM_WORLD, VolumeError);
|
|
#else
|
|
exit(VolumeError);
|
|
#endif
|
|
|
|
for(Int_t r = 0; r < domain.numReg(); r++)
|
|
{
|
|
Index_t numElemReg = domain.regElemSize(r);
|
|
// Index_t *regElemList = domain.regElemlist(r);
|
|
Int_t rep;
|
|
if(r < domain.numReg() / 2)
|
|
rep = 1;
|
|
else if(r < (domain.numReg() - (domain.numReg() + 15) / 20))
|
|
rep = 1 + domain.cost();
|
|
else
|
|
rep = 10 * (1 + domain.cost());
|
|
EvalEOSForElems(domain, vnewc.data(), numElemReg, r, rep);
|
|
}
|
|
}
|
|
}
|
|
|
|
static inline void
|
|
UpdateVolumesForElems(Domain& domain, Real_t v_cut, Index_t length)
|
|
{
|
|
if(length != 0)
|
|
{
|
|
Kokkos::parallel_for(
|
|
"UpdateVolumesForElems", length, KOKKOS_LAMBDA(const int i) {
|
|
Real_t tmpV = domain.vnew(i);
|
|
|
|
if(FABS(tmpV - Real_t(1.0)) < v_cut) tmpV = Real_t(1.0);
|
|
|
|
domain.v(i) = tmpV;
|
|
});
|
|
}
|
|
|
|
return;
|
|
}
|
|
|
|
static inline void
|
|
LagrangeElements(Domain& domain, Index_t numElem)
|
|
{
|
|
CalcLagrangeElements(domain);
|
|
|
|
CalcQForElems(domain);
|
|
|
|
ApplyMaterialPropertiesForElems(domain);
|
|
|
|
UpdateVolumesForElems(domain, domain.v_cut(), numElem);
|
|
}
|
|
|
|
static inline void
|
|
CalcCourantConstraintForElems(Domain& domain, Index_t length, Index_t r, Real_t qqc,
|
|
Real_t& dtcourant)
|
|
{
|
|
typedef Kokkos::View<Real_t*> view_real_t;
|
|
|
|
Real_t qqc2 = Real_t(64.0) * qqc * qqc;
|
|
Real_t dtcourant_tmp = dtcourant;
|
|
Index_t courant_elem = -1;
|
|
|
|
MinFinder result;
|
|
|
|
Kokkos::parallel_reduce(
|
|
"CalcCourantConstraintForElems", length,
|
|
KOKKOS_LAMBDA(const int i, MinFinder& minf) {
|
|
Index_t indx = domain.regElemlist(r, i);
|
|
Real_t dtf = domain.ss(indx) * domain.ss(indx);
|
|
|
|
if(domain.vdov(indx) < Real_t(0.))
|
|
{
|
|
dtf = dtf + qqc2 * domain.arealg(indx) * domain.arealg(indx) *
|
|
domain.vdov(indx) * domain.vdov(indx);
|
|
}
|
|
|
|
dtf = SQRT(dtf);
|
|
dtf = domain.arealg(indx) / dtf;
|
|
|
|
MinFinder tmp(dtf, i);
|
|
if(domain.vdov(indx) != Real_t(0.))
|
|
{
|
|
minf += tmp;
|
|
}
|
|
},
|
|
result);
|
|
|
|
dtcourant_tmp = result.val;
|
|
|
|
if(dtcourant_tmp > dtcourant)
|
|
{
|
|
dtcourant_tmp = dtcourant;
|
|
}
|
|
|
|
courant_elem = result.i;
|
|
|
|
if(courant_elem != -1)
|
|
{
|
|
dtcourant = dtcourant_tmp;
|
|
}
|
|
|
|
return;
|
|
}
|
|
|
|
static inline void
|
|
CalcHydroConstraintForElems(Domain& domain, Index_t length, Index_t r, Real_t dvovmax,
|
|
Real_t& dthydro)
|
|
{
|
|
typedef Kokkos::View<Real_t*> view_real_t;
|
|
|
|
Real_t dthydro_tmp = dthydro;
|
|
Index_t hydro_elem = -1;
|
|
MinFinder result;
|
|
|
|
Kokkos::parallel_reduce(
|
|
"CalcHydroConstraintForElems", length,
|
|
KOKKOS_LAMBDA(const int i, MinFinder& minf) {
|
|
Index_t indx = domain.regElemlist(r, i);
|
|
|
|
if(domain.vdov(indx) != Real_t(0.))
|
|
{
|
|
Real_t dtdvov = dvovmax / (FABS(domain.vdov(indx)) + Real_t(1.e-20));
|
|
|
|
MinFinder tmp(dtdvov, i);
|
|
if(domain.vdov(indx) != Real_t(0.))
|
|
{
|
|
minf += tmp;
|
|
}
|
|
}
|
|
},
|
|
result);
|
|
|
|
if(result.val > dthydro)
|
|
{
|
|
result.val = dthydro;
|
|
}
|
|
|
|
if(result.i != -1)
|
|
{
|
|
dthydro = result.val;
|
|
}
|
|
|
|
return;
|
|
}
|
|
|
|
static inline void
|
|
CalcTimeConstraintsForElems(Domain& domain)
|
|
{
|
|
domain.dtcourant() = 1.0e+20;
|
|
domain.dthydro() = 1.0e+20;
|
|
|
|
for(Index_t r = 0; r < domain.numReg(); ++r)
|
|
{
|
|
CalcCourantConstraintForElems(domain, domain.regElemSize(r), r, domain.qqc(),
|
|
domain.dtcourant());
|
|
|
|
CalcHydroConstraintForElems(domain, domain.regElemSize(r), r, domain.dvovmax(),
|
|
domain.dthydro());
|
|
}
|
|
}
|
|
|
|
static inline void
|
|
LagrangeLeapFrog(Domain& domain)
|
|
{
|
|
#ifdef SEDOV_SYNC_POS_VEL_LATE
|
|
Domain_member fieldData[6];
|
|
#endif
|
|
LagrangeNodal(domain);
|
|
|
|
#ifdef SEDOV_SYNC_POS_VEL_LATE
|
|
#endif
|
|
LagrangeElements(domain, domain.numElem());
|
|
|
|
#if USE_MPI
|
|
# ifdef SEDOV_SYNC_POS_VEL_LATE
|
|
CommRecv(domain, MSG_SYNC_POS_VEL, 6, domain.sizeX() + 1, domain.sizeY() + 1,
|
|
domain.sizeZ() + 1, false, false);
|
|
|
|
fieldData[0] = &Domain::x;
|
|
fieldData[1] = &Domain::y;
|
|
fieldData[2] = &Domain::z;
|
|
fieldData[3] = &Domain::xd;
|
|
fieldData[4] = &Domain::yd;
|
|
fieldData[5] = &Domain::zd;
|
|
|
|
CommSend(domain, MSG_SYNC_POS_VEL, 6, fieldData, domain.sizeX() + 1,
|
|
domain.sizeY() + 1, domain.sizeZ() + 1, false, false);
|
|
# endif
|
|
#endif
|
|
|
|
CalcTimeConstraintsForElems(domain);
|
|
|
|
#if USE_MPI
|
|
# ifdef SEDOV_SYNC_POS_VEL_LATE
|
|
CommSyncPosVel(domain);
|
|
# endif
|
|
#endif
|
|
}
|
|
|
|
int
|
|
main(int argc, char* argv[])
|
|
{
|
|
Int_t numRanks;
|
|
Int_t myRank;
|
|
struct cmdLineOpts opts;
|
|
|
|
#if USE_MPI
|
|
Domain_member fieldData;
|
|
|
|
MPI_Init(&argc, &argv);
|
|
MPI_Comm_size(MPI_COMM_WORLD, &numRanks);
|
|
MPI_Comm_rank(MPI_COMM_WORLD, &myRank);
|
|
#else
|
|
numRanks = 1;
|
|
myRank = 0;
|
|
#endif
|
|
|
|
Kokkos::initialize(argc, argv);
|
|
{
|
|
Kokkos::Tools::pushRegion("initialization");
|
|
opts.its = 9999999;
|
|
opts.nx = 30;
|
|
opts.numReg = 11;
|
|
opts.numFiles = (int) (numRanks + 10) / 9;
|
|
opts.showProg = 0;
|
|
opts.quiet = 0;
|
|
opts.viz = 0;
|
|
opts.balance = 1;
|
|
opts.cost = 1;
|
|
opts.do_atomic = 0;
|
|
|
|
ParseCommandLineOptions(argc, argv, myRank, &opts);
|
|
|
|
if(opts.do_atomic == 1)
|
|
do_atomic = 1;
|
|
else
|
|
do_atomic = 0;
|
|
|
|
if((myRank == 0) && (opts.quiet == 0))
|
|
{
|
|
printf("Running problem size %d^3 per domain until completion\n", opts.nx);
|
|
printf("Num processors: %d\n", numRanks);
|
|
printf("Total number of elements: %lld\n\n",
|
|
(long long int) (numRanks * opts.nx * opts.nx * opts.nx));
|
|
printf("To run other sizes, use -s <integer>.\n");
|
|
printf("To run a fixed number of iterations, use -i <integer>.\n");
|
|
printf("To run a more or less balanced region set, use -b <integer>.\n");
|
|
printf("To change the relative costs of regions, use -c <integer>.\n");
|
|
printf("To print out progress, use -p\n");
|
|
printf("To write an output file for VisIt, use -v\n");
|
|
printf("See help (-h) for more options\n\n");
|
|
}
|
|
|
|
Int_t col, row, plane, side;
|
|
InitMeshDecomp(numRanks, myRank, &col, &row, &plane, &side);
|
|
|
|
// Build the main data structure and initialize it
|
|
Domain locDom(numRanks, col, row, plane, opts.nx, side, opts.numReg, opts.balance,
|
|
opts.cost);
|
|
|
|
#if USE_MPI
|
|
fieldData = &Domain::nodalMass;
|
|
|
|
// Initial domain boundary communication
|
|
CommRecv(locDom, MSG_COMM_SBN, 1, locDom.sizeX() + 1, locDom.sizeY() + 1,
|
|
locDom.sizeZ() + 1, true, false);
|
|
CommSend(locDom, MSG_COMM_SBN, 1, &fieldData, locDom.sizeX() + 1,
|
|
locDom.sizeY() + 1, locDom.sizeZ() + 1, true, false);
|
|
CommSBN(locDom, 1, &fieldData);
|
|
|
|
// End initialization
|
|
MPI_Barrier(MPI_COMM_WORLD);
|
|
#endif
|
|
|
|
Kokkos::Tools::popRegion();
|
|
|
|
#if USE_MPI
|
|
double start = MPI_Wtime();
|
|
#else
|
|
timeval start;
|
|
gettimeofday(&start, nullptr);
|
|
#endif
|
|
|
|
uint32_t _time_incrp = 0;
|
|
uint32_t _leap_frogp = 0;
|
|
Kokkos::Tools::createProfileSection("TimeIncr", &_time_incrp);
|
|
Kokkos::Tools::createProfileSection("LeapFrog", &_leap_frogp);
|
|
|
|
while((locDom.time() < locDom.stoptime()) && (locDom.cycle() < opts.its))
|
|
{
|
|
Kokkos::Tools::startSection(_time_incrp);
|
|
// CAUSAL_BEGIN("Iteration")
|
|
TimeIncrement(locDom);
|
|
Kokkos::Tools::stopSection(_time_incrp);
|
|
|
|
Kokkos::Tools::startSection(_leap_frogp);
|
|
LagrangeLeapFrog(locDom);
|
|
Kokkos::Tools::stopSection(_leap_frogp);
|
|
CAUSAL_PROGRESS_NAMED("Iteration")
|
|
// CAUSAL_END("Iteration")
|
|
|
|
if((opts.showProg != 0) && (opts.quiet == 0) && (myRank == 0))
|
|
{
|
|
printf("cycle = %d, time = %e, dt=%e\n", locDom.cycle(),
|
|
double(locDom.time()), double(locDom.deltatime()));
|
|
}
|
|
Kokkos::Tools::markEvent("completed_timestep");
|
|
CAUSAL_PROGRESS
|
|
}
|
|
|
|
Kokkos::Tools::destroyProfileSection(_time_incrp);
|
|
Kokkos::Tools::destroyProfileSection(_leap_frogp);
|
|
|
|
double elapsed_time;
|
|
#if USE_MPI
|
|
elapsed_time = MPI_Wtime() - start;
|
|
#else
|
|
timeval end;
|
|
gettimeofday(&end, NULL);
|
|
elapsed_time = (double) (end.tv_sec - start.tv_sec) +
|
|
((double) (end.tv_usec - start.tv_usec)) / 1000000;
|
|
#endif
|
|
double elapsed_timeG;
|
|
#if USE_MPI
|
|
MPI_Reduce(&elapsed_time, &elapsed_timeG, 1, MPI_DOUBLE, MPI_MAX, 0,
|
|
MPI_COMM_WORLD);
|
|
#else
|
|
elapsed_timeG = elapsed_time;
|
|
#endif
|
|
|
|
Kokkos::Tools::pushRegion("finalization");
|
|
if(opts.viz)
|
|
{
|
|
DumpToVisit(locDom, opts.numFiles, myRank, numRanks);
|
|
}
|
|
|
|
if((myRank == 0) && (opts.quiet == 0))
|
|
{
|
|
VerifyAndWriteFinalOutput(elapsed_timeG, locDom, opts.nx, numRanks);
|
|
}
|
|
Kokkos::Tools::popRegion();
|
|
|
|
buffer = Kokkos::View<Real_t*>();
|
|
}
|
|
Kokkos::finalize();
|
|
#if USE_MPI
|
|
MPI_Finalize();
|
|
#endif
|
|
|
|
return 0;
|
|
}
|