Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Added artificial conductivity #256

Open
wants to merge 1 commit into
base: develop
Choose a base branch
from
+68 −11
Open

Added artificial conductivity #256

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
37 changes: 28 additions & 9 deletions sph/include/sph/kernel_ve/momentum_energy_kern.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,7 @@ template<typename T>
CUDA_DEVICE_HOST_FUN inline void
momentumAndEnergyJLoop(int i, T sincIndex, T K, const cstone::Box<T>& box, const int* neighbors, int neighborsCount,
const T* x, const T* y, const T* z, const T* vx, const T* vy, const T* vz, const T* h,
const T* m, const T* p, const T* c, const T* c11, const T* c12, const T* c13, const T* c22,
const T* m, const T* p, const T* c, const T* u, const T* c11, const T* c12, const T* c13, const T* c22,
const T* c23, const T* c33, const T Atmin, const T Atmax, const T ramp, const T* wh,
const T* whd, const T* kx, const T* xm, const T* alpha, const T* gradh, T* grad_P_x, T* grad_P_y,
T* grad_P_z, T* du, T* maxvsignal)
Expand All @@ -61,6 +61,7 @@ momentumAndEnergyJLoop(int i, T sincIndex, T K, const cstone::Box<T>& box, const
T hi = h[i];
T mi = m[i];
T ci = c[i];
T ui = u[i];
T kxi = kx[i];

T alpha_i = alpha[i];
Expand All @@ -79,7 +80,8 @@ momentumAndEnergyJLoop(int i, T sincIndex, T K, const cstone::Box<T>& box, const

T maxvsignali = 0.0;
T momentum_x = 0.0, momentum_y = 0.0, momentum_z = 0.0, energy = 0.0;
T a_visc_energy = 0.0;
T a_visc_gradp = 0.0;
T a_heat_cond = 0.0;

T c11i = c11[i];
T c12i = c12[i];
Expand All @@ -101,6 +103,10 @@ momentumAndEnergyJLoop(int i, T sincIndex, T K, const cstone::Box<T>& box, const
T r2 = rx * rx + ry * ry + rz * rz;
T dist = std::sqrt(r2);

T ux_ij = rx / dist;
T uy_ij = ry / dist;
T uz_ij = rz / dist;

T vx_ij = vxi - vx[j];
T vy_ij = vyi - vy[j];
T vz_ij = vzi - vz[j];
Expand Down Expand Up @@ -134,21 +140,24 @@ momentumAndEnergyJLoop(int i, T sincIndex, T K, const cstone::Box<T>& box, const

T mj = m[j];
T cj = c[j];
T uj = u[j];
T kxj = kx[j];
T xmassj = xm[j];
T rhoj = kxj * mj / xmassj;

T alpha_j = alpha[j];

T wij = rv / dist;
T viscosity_ij = artificial_viscosity(alpha_i, alpha_j, ci, cj, wij);
T viscosity_jj = viscosity_ij;
T wij = rv / dist;
T delta_u = ui - uj;
T proj = p[j] / (kxj * mj * mj * gradh[j]);
T viscosity_ij = artificial_viscosity(alpha_i, alpha_j, ci, cj, wij);
T heat_conduction = AV_heat_conduction(wij, rhoi, rhoj, proi, proj, delta_u);
T viscosity_jj = viscosity_ij;

// For time-step calculations
T vijsignal = ci + cj - T(3) * wij;
maxvsignali = (vijsignal > maxvsignali) ? vijsignal : maxvsignali;

T proj = p[j] / (kxj * mj * mj * gradh[j]);

T Atwood = (std::abs(rhoi - rhoj)) / (rhoi + rhoj);
if (Atwood < Atmin)
Expand Down Expand Up @@ -179,14 +188,23 @@ momentumAndEnergyJLoop(int i, T sincIndex, T K, const cstone::Box<T>& box, const
T a_visc_y = T(0.5) * (a_visc * termA2_i + b_visc * termA2_j);
T a_visc_z = T(0.5) * (a_visc * termA3_i + b_visc * termA3_j);

T a_heat = voli * mj / mi * heat_conduction;
T b_heat = volj * heat_conduction;
T a_heat_x = T(0.5) * (a_heat * termA1_i + b_heat * termA1_j);
T a_heat_y = T(0.5) * (a_heat * termA2_i + b_heat * termA2_j);
T a_heat_z = T(0.5) * (a_heat * termA3_i + b_heat * termA3_j);

{
// T a = Wi * (mj_pro_i + viscosity_ij * mi_roi);
// T b = mj_roj_Wj * (p[j] / (roj * gradh_j) + viscosity_ij);

momentum_x += momentum_i * termA1_i + momentum_j * termA1_j + a_visc_x;
momentum_y += momentum_i * termA2_i + momentum_j * termA2_j + a_visc_y;
momentum_z += momentum_i * termA3_i + momentum_j * termA3_j + a_visc_z;
a_visc_energy += a_visc_x * vx_ij + a_visc_y * vy_ij + a_visc_z * vz_ij;

a_visc_gradp += a_visc_x * vx_ij + a_visc_y * vy_ij + a_visc_z * vz_ij;
a_heat_cond += a_heat_x * ux_ij + a_heat_y * uy_ij + a_heat_z * uz_ij;

}
{
// T a = Wi * (2.0 * mj_pro_i + viscosity_ij * mi_roi);
Expand All @@ -200,8 +218,9 @@ momentumAndEnergyJLoop(int i, T sincIndex, T K, const cstone::Box<T>& box, const

// with the choice of calculating coordinate (r) and velocity (v_ij) differences as i - j,
// we add the negative sign only here at the end instead of to termA123_ij in each iteration
a_visc_energy = std::max(T(0), a_visc_energy);
du[i] = K * (energy + T(0.5) * a_visc_energy); // factor of 2 already removed from 2P/rho
a_visc_gradp = std::max(T(0), a_visc_gradp);
du[i] = K * (energy + T(0.5) * a_visc_gradp + a_heat_cond); // factor of 2 already removed from 2P/rho

// grad_P_xyz is stored as the acceleration, accel = -grad_P / rho
grad_P_x[i] = -K * momentum_x;
grad_P_y[i] = -K * momentum_y;
Expand Down
34 changes: 34 additions & 0 deletions sph/include/sph/kernels.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -109,5 +109,39 @@ CUDA_DEVICE_FUN inline T artificial_viscosity(T alpha_i, T alpha_j, T c_i, T c_j

return viscosity_ij;
}
/*! @brief calculate the AV heat conduction between a pair of two particles
*
* @tparam T float or double
* @param p_i baryonic pressure of particle i
* @param p_j baryonic pressure of particle j
* @param rho_i baryonic density of particle i
* @param rho_j baryonic density of particle j
* @param delta_u internal energy difference between particles i and j (u_i - u_j)
* @param w_ij relative velocity (v_i - v_j), projected onto the connecting axis (r_i - r_j)
* @return the heat conduction AV term
*/
template<typename T>
CUDA_DEVICE_FUN inline T AV_heat_conduction(T w_ij, T rho_i, T rho_j, T p_i, T p_j, T delta_u)
{
constexpr T alfa_u = T(0.05);

T heat_conduction = T(0.0);
if (w_ij < T(0.0))
{
T vij_signal_u = T(0.0);

// with gravity
vij_signal_u = abs(w_ij);

// without gravity
//T rho_ij = T(0.5) * (rho_i + rho_j);
//vij_signal_u = std::sqrt(abs(p_i - p_j) / rho_ij);

heat_conduction = alfa_u * vij_signal_u * delta_u;
}

return heat_conduction;
}


} // namespace sphexa
2 changes: 2 additions & 0 deletions sph/include/sph/momentum_energy_ve.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@ void computeGradPVeImpl(size_t startIndex, size_t endIndex, size_t ngmax, Datase
const T* vy = d.vy.data();
const T* vz = d.vz.data();
const T* c = d.c.data();
const T* u = d.u.data();
const T* p = d.p.data();
const T* alpha = d.alpha.data();
const T* gradh = d.gradh.data();
Expand Down Expand Up @@ -77,6 +78,7 @@ void computeGradPVeImpl(size_t startIndex, size_t endIndex, size_t ngmax, Datase
m,
p,
c,
u,
c11,
c12,
c13,
Expand Down
6 changes: 4 additions & 2 deletions sph/test/kernel_ve/momentum_energy_kern.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -70,6 +70,7 @@ TEST(MomentumEnergy, JLoop)
std::vector<T> vz{0.091, -0.081, 0.071, -0.061, 0.055};

std::vector<T> c{0.4, 0.5, 0.6, 0.7, 0.8};
std::vector<T> u{0.4, 0.5, 0.6, 0.7, 0.8};
std::vector<T> p{0.2, 0.3, 0.4, 0.5, 0.6};

std::vector<T> alpha{1.0, 0.05, 0.3, 0.5, 0.3};
Expand Down Expand Up @@ -119,6 +120,7 @@ TEST(MomentumEnergy, JLoop)
m.data(),
p.data(),
c.data(),
u.data(),
c11.data(),
c12.data(),
c13.data(),
Expand All @@ -143,6 +145,6 @@ TEST(MomentumEnergy, JLoop)
EXPECT_NEAR(grad_Px, 4.6852624676440924e-1, 1e-10);
EXPECT_NEAR(grad_Py, -8.2810161944474575e-2, 1e-10);
EXPECT_NEAR(grad_Pz, 5.209843022360216e-1, 1e-10);
EXPECT_NEAR(du, -3.8445778269613888e-3, 1e-10);
EXPECT_NEAR(maxvsignal, 1.4112466829, 1e-10);
EXPECT_NEAR(du, -3.6517141514543563e-3, 1e-10);
EXPECT_NEAR(maxvsignal, 1.4112466828564341, 1e-10);
}