Cholla 3.0.1-dev
Cholla - Massively parallel hydro on GPUs
Loading...
Searching...
No Matches
static_grav.h
1
6#pragma once
7
8#include <math.h> // provides sqrt log cos sin atan etc.
9#include <stdio.h>
10
11#include "../global/global.h" // provides GN etc.
12
13// Work around lack of pow(Real,int) in Hip Clang for Rocm 3.5
14static inline __device__ Real pow2(const Real x) { return x * x; }
15
16inline __device__ void calc_g_1D(int xid, int x_off, int n_ghost, int custom_grav, Real dx, Real xbound, Real *gx)
17{
18 Real x_pos, r_disk, r_halo;
19 x_pos = (x_off + xid - n_ghost + 0.5) * dx + xbound;
20 // set gravity field according to parameter file input
21 switch (custom_grav) {
22 case 1:
23 // 1D NFW halo & Miyamoto-Nagai disk
24 // for disk components, calculate polar r
25 // r_disk = 0.220970869121;
26 // r_disk = 6.85009694274;
27 r_disk = 13.9211647546;
28 // r_disk = 20.9922325665;
29 // for halo, calculate spherical r
30 r_halo = sqrt(x_pos * x_pos + r_disk * r_disk);
31
32 // set properties of halo and disk (these must match initial conditions)
33 Real a_disk_z, a_halo, M_vir, M_d, R_vir, R_d, z_d, R_h, M_h, c_vir, phi_0_h, x;
34 M_vir = 1.0e12; // viral mass of MW in M_sun
35 M_d = 6.5e10; // mass of disk in M_sun
36 M_h = M_vir - M_d; // halo mass in M_sun
37 R_vir = 261; // viral radius in kpc
38 c_vir = 20.0; // halo concentration
39 R_h = R_vir / c_vir; // halo scale length in kpc
40 R_d = 3.5; // disk scale length in kpc
41 z_d = 3.5 / 5.0; // disk scale height in kpc
42 phi_0_h = GN * M_h / (log(1.0 + c_vir) - c_vir / (1.0 + c_vir));
43 x = r_halo / R_h;
44
45 // calculate acceleration due to NFW halo & Miyamoto-Nagai disk
46 a_halo = -phi_0_h * (log(1 + x) - x / (1 + x)) / (r_halo * r_halo);
47 a_disk_z =
48 -GN * M_d * x_pos * (R_d + sqrt(x_pos * x_pos + z_d * z_d)) /
49 (pow(r_disk * r_disk + pow2(R_d + sqrt(x_pos * x_pos + z_d * z_d)), 1.5) * sqrt(x_pos * x_pos + z_d * z_d));
50
51 // total acceleration is the sum of the halo + disk components
52 *gx = (x_pos / r_halo) * a_halo + a_disk_z;
53 break;
54 default:
55 *gx = 0;
56 }
57 return;
58}
59
60inline __device__ void calc_g_2D(int xid, int yid, int x_off, int y_off, int n_ghost, int custom_grav, Real dx, Real dy,
61 Real xbound, Real ybound, Real *gx, Real *gy)
62{
63 Real x_pos, y_pos, r, phi;
64 // use the subgrid offset and global boundaries to calculate absolute
65 // positions on the grid
66 x_pos = (x_off + xid - n_ghost + 0.5) * dx + xbound;
67 y_pos = (y_off + yid - n_ghost + 0.5) * dy + ybound;
68 // for Gresho and disks, also need r & phi
69 r = sqrt(x_pos * x_pos + y_pos * y_pos);
70 phi = atan2(y_pos, x_pos);
71 switch (custom_grav) {
72 case 1:
73 // Gresho vortex
74 // set acceleration to balance v_phi in Gresho problem
75 if (r < 0.2) {
76 *gx = -cos(phi) * 25.0 * r;
77 *gy = -sin(phi) * 25.0 * r;
78 } else if (r >= 0.2 && r < 0.4) {
79 *gx = -cos(phi) * (4.0 - 20.0 * r + 25.0 * r * r) / r;
80 *gy = -sin(phi) * (4.0 - 20.0 * r + 25.0 * r * r) / r;
81 } else {
82 *gx = 0.0;
83 *gy = 0.0;
84 }
85 break;
86 case 2:
87 // Rayleigh-Taylor instability
88 *gx = 0;
89 *gy = -1;
90 break;
91 case 3:
92 // 2D disk in keplerian rotation
93 Real M;
94 M = 1 * MSUN_CGS;
95 *gx = -cos(phi) * GN * M / (r * r);
96 *gy = -sin(phi) * GN * M / (r * r);
97 break;
98 case 4:
99 // set gravitational acceleration for Kuzmin disk + NFW halo
100 Real a_d, a_h, a, M_vir, M_d, R_vir, R_d, R_s, M_h, c_vir, x;
101 M_vir = 1.0e12; // viral mass of MW in M_sun
102 M_d = 6.5e10; // mass of disk in M_sun (assume all gas)
103 M_h = M_vir - M_d; // halo mass in M_sun
104 R_vir = 261; // viral radius in kpc
105 c_vir = 20; // halo concentration
106 R_s = R_vir / c_vir; // halo scale length in kpc
107 R_d = 3.5; // disk scale length in kpc
108
109 // calculate acceleration
110 x = r / R_s;
111 a_d = GN * M_d * r * pow(r * r + R_d * R_d, -1.5);
112 a_h = GN * M_h * (log(1 + x) - x / (1 + x)) / ((log(1 + c_vir) - c_vir / (1 + c_vir)) * r * r);
113 a = a_d + a_h;
114
115 *gx = -cos(phi) * a;
116 *gy = -sin(phi) * a;
117 break;
118 default:
119 *gx = 0;
120 *gy = 0;
121 }
122
123 return;
124}
125
126inline __device__ void calc_g_3D(int xid, int yid, int zid, int x_off, int y_off, int z_off, int n_ghost,
127 int custom_grav, Real dx, Real dy, Real dz, Real xbound, Real ybound, Real zbound,
128 Real *gx, Real *gy, Real *gz)
129{
130 Real x_pos, y_pos, z_pos, r_disk, r_halo;
131 // use the subgrid offset and global boundaries to calculate absolute
132 // positions on the grid
133 x_pos = (x_off + xid - n_ghost + 0.5) * dx + xbound;
134 y_pos = (y_off + yid - n_ghost + 0.5) * dy + ybound;
135 z_pos = (z_off + zid - n_ghost + 0.5) * dz + zbound;
136
137 // for disk components, calculate polar r
138 r_disk = sqrt(x_pos * x_pos + y_pos * y_pos);
139 // for halo, calculate spherical r
140 r_halo = sqrt(x_pos * x_pos + y_pos * y_pos + z_pos * z_pos);
141 Real a_disk_r, a_disk_z, a_halo, a_halo_r, a_halo_z;
142 Real M_vir, M_d, R_vir, R_d, z_d, R_h, M_h, c_vir, phi_0_h, x;
143 switch (custom_grav) {
144 case 1:
145 // Milky way disk model
146 // set properties of halo and disk (these must match initial conditions)
147
148 M_vir = 1.0e12; // viral mass of in M_sun
149 M_d = 6.5e10; // viral mass of in M_sun
150 R_d = 3.5; // disk scale length in kpc
151 z_d = 3.5 / 5.0; // disk scale height in kpc
152 R_vir = 261.; // virial radius in kpc
153 c_vir = 20.0; // halo concentration
154
155 M_h = M_vir - M_d; // halo mass in M_sun
156 R_h = R_vir / c_vir; // halo scale length in kpc
157 phi_0_h = GN * M_h / (log(1.0 + c_vir) - c_vir / (1.0 + c_vir));
158 x = r_halo / R_h;
159
160 // calculate acceleration due to NFW halo & Miyamoto-Nagai disk
161 a_halo = -phi_0_h * (log(1 + x) - x / (1 + x)) / (r_halo * r_halo);
162 a_halo_r = a_halo * (r_disk / r_halo);
163 a_halo_z = a_halo * (z_pos / r_halo);
164 a_disk_r = -GN * M_d * r_disk * pow(r_disk * r_disk + pow2(R_d + sqrt(z_pos * z_pos + z_d * z_d)), -1.5);
165 a_disk_z =
166 -GN * M_d * z_pos * (R_d + sqrt(z_pos * z_pos + z_d * z_d)) /
167 (pow(r_disk * r_disk + pow2(R_d + sqrt(z_pos * z_pos + z_d * z_d)), 1.5) * sqrt(z_pos * z_pos + z_d * z_d));
168
169 // total acceleration is the sum of the halo + disk components
170 *gx = (x_pos / r_disk) * (a_disk_r + a_halo_r);
171 *gy = (y_pos / r_disk) * (a_disk_r + a_halo_r);
172 *gz = a_disk_z + a_halo_z;
173 break;
174 case 2:
175 // M82 model
176 // set properties of halo and disk (these must match initial conditions)
177
178 M_vir = 5.0e10; // viral mass of in M_sun
179 M_d = 1.0e10; // mass of disk in M_sun
180 R_d = 0.8; // disk scale length in kpc
181 z_d = 0.15; // disk scale height in kpc
182 R_vir = R_d / 0.015; // viral radius in kpc
183 c_vir = 10.0; // halo concentration
184
185 M_h = M_vir - M_d; // halo mass in M_sun
186 R_h = R_vir / c_vir; // halo scale length in kpc
187 phi_0_h = GN * M_h / (log(1.0 + c_vir) - c_vir / (1.0 + c_vir));
188 x = r_halo / R_h;
189
190 // calculate acceleration due to NFW halo & Miyamoto-Nagai disk
191 a_halo = -phi_0_h * (log(1 + x) - x / (1 + x)) / (r_halo * r_halo);
192 a_halo_r = a_halo * (r_disk / r_halo);
193 a_halo_z = a_halo * (z_pos / r_halo);
194 a_disk_r = -GN * M_d * r_disk * pow(r_disk * r_disk + pow2(R_d + sqrt(z_pos * z_pos + z_d * z_d)), -1.5);
195 a_disk_z =
196 -GN * M_d * z_pos * (R_d + sqrt(z_pos * z_pos + z_d * z_d)) /
197 (pow(r_disk * r_disk + pow2(R_d + sqrt(z_pos * z_pos + z_d * z_d)), 1.5) * sqrt(z_pos * z_pos + z_d * z_d));
198
199 // total acceleration is the sum of the halo + disk components
200 *gx = (x_pos / r_disk) * (a_disk_r + a_halo_r);
201 *gy = (y_pos / r_disk) * (a_disk_r + a_halo_r);
202 *gz = a_disk_z + a_halo_z;
203 break;
204 default:
205 *gx = 0;
206 *gy = 0;
207 *gz = 0;
208 }
209 return;
210}
Generated by doxygen 1.9.8