artificial_viscosity.cc 4.99 KB
Newer Older
Volker Springel's avatar
Volker Springel committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
/*******************************************************************************
 * \copyright   This file is part of the GADGET4 N-body/SPH code developed
 * \copyright   by Volker Springel. Copyright (C) 2014-2020 by Volker Springel
 * \copyright   (vspringel@mpa-garching.mpg.de) and all contributing authors.
 *******************************************************************************/

/*! \file  artificial_viscosity.cc
 *
 *  \brief Calculates time-dependent artificial viscosity parameter
 */

#include "gadgetconfig.h"

#include <gsl/gsl_linalg.h>

#include "../data/sph_particle_data.h"

/*! This file contains the function for the time-dependent artificial viscosity
 */

#ifdef IMPROVED_VELOCITY_GRADIENTS
void sph_particle_data::set_velocity_gradients(void)
{
#ifdef ONEDIMS
  if(fabs(dpos.dx_dx) > 0)
    {
      dvel[0][0] = dvel[0][0] / dpos.dx_dx;
      DivVel     = dvel[0][0];
    }
  else
    {
      DivVel = dvel[0][0] / Density;
    }
#elif defined(TWODIMS)
  double det = dpos.dx_dx * dpos.dy_dy - dpos.dx_dy * dpos.dx_dy;
  if(fabs(det) > 0)
    {
      double m11_inv = dpos.dy_dy / det;
      double m12_inv = -dpos.dx_dy / det;
      double m21_inv = -dpos.dx_dy / det;
      double m22_inv = dpos.dx_dx / det;

      double y11 = dvel[0][0];
      double y21 = dvel[0][1];
      double y12 = dvel[1][0];
      double y22 = dvel[1][1];

      dvel[0][0] = m11_inv * y11 + m12_inv * y21;
      dvel[0][1] = m21_inv * y11 + m22_inv * y21;
      dvel[1][0] = m11_inv * y12 + m12_inv * y22;
      dvel[1][1] = m21_inv * y12 + m22_inv * y22;
      DivVel     = dvel[0][0] + dvel[1][1];
      CurlVel    = fabs(dvel[0][1] - dvel[1][0]);
    }
  else
    {
      DivVel  = (dvel[0][0] + dvel[1][1]) / Density;
      CurlVel = fabs((dvel[1][0] - dvel[0][1]) / Density);
    }
#else
  gsl_matrix* distance_matrix = gsl_matrix_alloc(3, 3);
  gsl_matrix_set(distance_matrix, 0, 0, dpos.dx_dx);
  gsl_matrix_set(distance_matrix, 0, 1, dpos.dx_dy);
  gsl_matrix_set(distance_matrix, 0, 2, dpos.dx_dz);
  gsl_matrix_set(distance_matrix, 1, 0, dpos.dx_dy);
  gsl_matrix_set(distance_matrix, 1, 1, dpos.dy_dy);
  gsl_matrix_set(distance_matrix, 1, 2, dpos.dy_dz);
  gsl_matrix_set(distance_matrix, 2, 0, dpos.dx_dz);
  gsl_matrix_set(distance_matrix, 2, 1, dpos.dy_dz);
  gsl_matrix_set(distance_matrix, 2, 2, dpos.dz_dz);

  int sign                = 1;
  gsl_permutation* permut = gsl_permutation_alloc(3);
  gsl_linalg_LU_decomp(distance_matrix, permut, &sign);

  double det = gsl_linalg_LU_det(distance_matrix, sign);

  if(fabs(det) > 0)
    {
      gsl_matrix* inv_distance_matrix = gsl_matrix_alloc(3, 3);
      gsl_linalg_LU_invert(distance_matrix, permut, inv_distance_matrix);

      double m_inv[3][3];
      for(int i = 0; i < 3; i++)
        {
          for(int j = 0; j < 3; j++)
            {
              m_inv[i][j] = gsl_matrix_get(inv_distance_matrix, i, j);
            }
        }

      double y[3][3];
      for(int i = 0; i < 3; i++)
        {
          for(int j = 0; j < 3; j++)
            {
              y[i][j] = dvel[i][j];
            }
        }

      for(int i = 0; i < 3; i++)
        {
          for(int j = 0; j < 3; j++)
            {
              dvel[i][j] = 0;
              for(int k = 0; k < 3; k++)
                {
                  dvel[i][j] += m_inv[j][k] * y[k][i];
                }
            }
        }

      DivVel  = dvel[0][0] + dvel[1][1] + dvel[2][2];
      Rot[0]  = dvel[2][1] - dvel[1][2];
      Rot[1]  = dvel[0][2] - dvel[2][0];
      Rot[2]  = dvel[1][0] - dvel[0][1];
      CurlVel = sqrt(Rot[0] * Rot[0] + Rot[1] * Rot[1] + Rot[2] * Rot[2]);

      gsl_permutation_free(permut);
      gsl_matrix_free(distance_matrix);
      gsl_matrix_free(inv_distance_matrix);
    }
  else
    {
      DivVel  = (dvel[0][0] + dvel[1][1] + dvel[2][2]) / Density;
      Rot[0]  = (dvel[2][1] - dvel[1][2]) / Density;
      Rot[1]  = (dvel[0][2] - dvel[2][0]) / Density;
      Rot[2]  = (dvel[1][0] - dvel[0][1]) / Density;
      CurlVel = sqrt(Rot[0] * Rot[0] + Rot[1] * Rot[1] + Rot[2] * Rot[2]);
    }

#endif
}
#endif

#ifdef TIMEDEP_ART_VISC
void sph_particle_data::set_viscosity_coefficient(double dt)
{
  double dDivVel_dt = dt > 0 ? (DivVel - DivVelOld) / (dt) : 0;
  dDivVel_dt *= All.cf_a2inv;
  double shockIndicator = -dDivVel_dt > 0 ? -dDivVel_dt : 0;
  double hsml           = Hsml * All.cf_atime;
  double Hsml2          = hsml * hsml;
Oliver Zier's avatar
Oliver Zier committed
144
  double alpha_tar      = (Hsml2 * shockIndicator) / (Hsml2 * shockIndicator + Csnd * Csnd) * All.ArtBulkViscConst;
Volker Springel's avatar
Volker Springel committed
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165

  double DivVel2       = DivVel * DivVel;
  double CurlVel2      = CurlVel * CurlVel;
  double CsndOverHsml2 = (Csnd / Hsml) * (Csnd / Hsml);
  double limiter       = DivVel2 / (DivVel2 + CurlVel2 + 0.00001 * CsndOverHsml2);
#ifdef NO_SHEAR_VISCOSITY_LIMITER
  limiter = 1.;
#endif

  if(Alpha < alpha_tar)
    {
      Alpha = alpha_tar * limiter;
      return;
    }
  double DecayTime = 10. * Hsml / decayVel;
  Alpha            = limiter * (alpha_tar + (Alpha - alpha_tar) * exp(-dt / DecayTime));
  if(Alpha < All.AlphaMin)
    Alpha = All.AlphaMin;
}

#endif