GADGET-4
artificial_viscosity.cc
Go to the documentation of this file.
1/*******************************************************************************
2 * \copyright This file is part of the GADGET4 N-body/SPH code developed
3 * \copyright by Volker Springel. Copyright (C) 2014-2020 by Volker Springel
4 * \copyright (vspringel@mpa-garching.mpg.de) and all contributing authors.
5 *******************************************************************************/
6
12#include "gadgetconfig.h"
13
14#include <gsl/gsl_linalg.h>
15
16#include "../data/sph_particle_data.h"
17
21#ifdef IMPROVED_VELOCITY_GRADIENTS
22void sph_particle_data::set_velocity_gradients(void)
23{
24#ifdef ONEDIMS
25 if(fabs(dpos.dx_dx) > 0)
26 {
27 dvel[0][0] = dvel[0][0] / dpos.dx_dx;
28 DivVel = dvel[0][0];
29 }
30 else
31 {
32 DivVel = dvel[0][0] / Density;
33 }
34#elif defined(TWODIMS)
35 double det = dpos.dx_dx * dpos.dy_dy - dpos.dx_dy * dpos.dx_dy;
36 if(fabs(det) > 0)
37 {
38 double m11_inv = dpos.dy_dy / det;
39 double m12_inv = -dpos.dx_dy / det;
40 double m21_inv = -dpos.dx_dy / det;
41 double m22_inv = dpos.dx_dx / det;
42
43 double y11 = dvel[0][0];
44 double y21 = dvel[0][1];
45 double y12 = dvel[1][0];
46 double y22 = dvel[1][1];
47
48 dvel[0][0] = m11_inv * y11 + m12_inv * y21;
49 dvel[0][1] = m21_inv * y11 + m22_inv * y21;
50 dvel[1][0] = m11_inv * y12 + m12_inv * y22;
51 dvel[1][1] = m21_inv * y12 + m22_inv * y22;
52 DivVel = dvel[0][0] + dvel[1][1];
53 CurlVel = fabs(dvel[0][1] - dvel[1][0]);
54 }
55 else
56 {
57 DivVel = (dvel[0][0] + dvel[1][1]) / Density;
58 CurlVel = fabs((dvel[1][0] - dvel[0][1]) / Density);
59 }
60#else
61 gsl_matrix* distance_matrix = gsl_matrix_alloc(3, 3);
62 gsl_matrix_set(distance_matrix, 0, 0, dpos.dx_dx);
63 gsl_matrix_set(distance_matrix, 0, 1, dpos.dx_dy);
64 gsl_matrix_set(distance_matrix, 0, 2, dpos.dx_dz);
65 gsl_matrix_set(distance_matrix, 1, 0, dpos.dx_dy);
66 gsl_matrix_set(distance_matrix, 1, 1, dpos.dy_dy);
67 gsl_matrix_set(distance_matrix, 1, 2, dpos.dy_dz);
68 gsl_matrix_set(distance_matrix, 2, 0, dpos.dx_dz);
69 gsl_matrix_set(distance_matrix, 2, 1, dpos.dy_dz);
70 gsl_matrix_set(distance_matrix, 2, 2, dpos.dz_dz);
71
72 int sign = 1;
73 gsl_permutation* permut = gsl_permutation_alloc(3);
74 gsl_linalg_LU_decomp(distance_matrix, permut, &sign);
75
76 double det = gsl_linalg_LU_det(distance_matrix, sign);
77
78 if(fabs(det) > 0)
79 {
80 gsl_matrix* inv_distance_matrix = gsl_matrix_alloc(3, 3);
81 gsl_linalg_LU_invert(distance_matrix, permut, inv_distance_matrix);
82
83 double m_inv[3][3];
84 for(int i = 0; i < 3; i++)
85 {
86 for(int j = 0; j < 3; j++)
87 {
88 m_inv[i][j] = gsl_matrix_get(inv_distance_matrix, i, j);
89 }
90 }
91
92 double y[3][3];
93 for(int i = 0; i < 3; i++)
94 {
95 for(int j = 0; j < 3; j++)
96 {
97 y[i][j] = dvel[i][j];
98 }
99 }
100
101 for(int i = 0; i < 3; i++)
102 {
103 for(int j = 0; j < 3; j++)
104 {
105 dvel[i][j] = 0;
106 for(int k = 0; k < 3; k++)
107 {
108 dvel[i][j] += m_inv[j][k] * y[k][i];
109 }
110 }
111 }
112
113 DivVel = dvel[0][0] + dvel[1][1] + dvel[2][2];
114 Rot[0] = dvel[2][1] - dvel[1][2];
115 Rot[1] = dvel[0][2] - dvel[2][0];
116 Rot[2] = dvel[1][0] - dvel[0][1];
117 CurlVel = sqrt(Rot[0] * Rot[0] + Rot[1] * Rot[1] + Rot[2] * Rot[2]);
118
119 gsl_permutation_free(permut);
120 gsl_matrix_free(distance_matrix);
121 gsl_matrix_free(inv_distance_matrix);
122 }
123 else
124 {
125 DivVel = (dvel[0][0] + dvel[1][1] + dvel[2][2]) / Density;
126 Rot[0] = (dvel[2][1] - dvel[1][2]) / Density;
127 Rot[1] = (dvel[0][2] - dvel[2][0]) / Density;
128 Rot[2] = (dvel[1][0] - dvel[0][1]) / Density;
129 CurlVel = sqrt(Rot[0] * Rot[0] + Rot[1] * Rot[1] + Rot[2] * Rot[2]);
130 }
131
132#endif
133}
134#endif
135
136#ifdef TIMEDEP_ART_VISC
137void sph_particle_data::set_viscosity_coefficient(double dt)
138{
139 double dDivVel_dt = dt > 0 ? (DivVel - DivVelOld) / (dt) : 0;
140 dDivVel_dt *= All.cf_a2inv; //now in physical coordinates
141 double shockIndicator = -dDivVel_dt > 0 ? -dDivVel_dt : 0;
142 double hsml_p = Hsml * All.cf_atime;
143 double hsml2_p = hsml_p * hsml_p;
144 double csnd_p = Csnd * All.cf_afac3;
145 double alpha_tar = (hsml2_p * shockIndicator) / (hsml2_p * shockIndicator + csnd_p * csnd_p) * All.ArtBulkViscConst;
146
147 double DivVel2 = (DivVel * All.cf_a2inv) * (DivVel * All.cf_a2inv);
148 double CurlVel2 = (CurlVel * All.cf_a2inv) * (CurlVel * All.cf_a2inv);
149 double CsndOverHsml2 = (csnd_p / hsml_p) * (csnd_p / hsml_p);
150 double limiter = DivVel2 / (DivVel2 + CurlVel2 + 0.00001 * CsndOverHsml2);
151#ifdef NO_SHEAR_VISCOSITY_LIMITER
152 limiter = 1.;
153#endif
154
155 if(Alpha < alpha_tar)
156 {
157 Alpha = alpha_tar * limiter;
158 return;
159 }
160
161 double devayVel_p = decayVel * All.cf_afac3; //has the same a factor as sound speed
162
163 double DecayTime = 10. * hsml_p / devayVel_p;
164 Alpha = limiter * (alpha_tar + (Alpha - alpha_tar) * exp(-dt / DecayTime));
165 if(Alpha < All.AlphaMin)
166 Alpha = All.AlphaMin;
167}
168
169#endif
global_data_all_processes All
Definition: main.cc:40
expr exp(half arg)
Definition: half.hpp:2724
half fabs(half arg)
Definition: half.hpp:2627
expr sqrt(half arg)
Definition: half.hpp:2777