noe/pipeline
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

[postprocessor] [bug] take dx and dy into account in gradients

Browse Source
This commit is contained in:
Noe Brucy
2021-02-11 09:06:25 +01:00
parent 67bff23457
commit 4d533ab464
1 changed files with 18 additions and 10 deletions
Download Patch File Download Diff File Expand all files Collapse all files
postprocessor.py
+18 -10
View File
@@ -902,17 +902,25 @@ class PostProcessor(HDF5Container):
dy = (im_extent[3] - im_extent[2]) / map_size dy = (im_extent[3] - im_extent[2]) / map_size
# Physical coordinates of the center of the cells # Physical coordinates of the center of the cells
x = np.linspace(im_extent[0], im_extent[1], map_size) + 0.5 * dx x = (
y = np.linspace(im_extent[2], im_extent[3], map_size) + 0.5 * dy np.linspace(im_extent[0], im_extent[1], map_size)
+ 0.5 * dx
- pos_star[0]
)
y = (
np.linspace(im_extent[2], im_extent[3], map_size)
+ 0.5 * dy
- pos_star[1]
)
xx, yy = np.meshgrid(x, y) xx, yy = np.meshgrid(x, y)
rr = np.sqrt((xx - pos_star[0]) ** 2 + (yy - pos_star[1]) ** 2) rr = np.sqrt(xx ** 2 + yy ** 2)
# Compute kappa**2 = (2omega/R)*d(R**2*omega)/dR # Compute kappa**2 = (2omega/R)*d(R**2*omega)/dR
Gx, Gy = np.gradient(rr ** 2 * omega) Gy, Gx = np.gradient(rr ** 2 * omega, dy, dx, edge_order=2)
Gr = (xx * Gx + yy * Gy) / rr Gr = (xx * Gx + yy * Gy) / rr
kappa_square = (2 * omega / rr) * Gr kappa_square = (2 * omega / rr) * Gr
kappa = np.square(kappa_square) kappa = np.sqrt(kappa_square)
map_Q = (np.sqrt(cs2) * kappa) / (np.pi * self.G * coldens) map_Q = (np.sqrt(cs2) * kappa) / (np.pi * self.G * coldens)
@@ -1303,22 +1311,22 @@ class PostProcessor(HDF5Container):
dy = (im_extent[3] - im_extent[2]) / map_size dy = (im_extent[3] - im_extent[2]) / map_size
# Physical coordinates of the center of the cells # Physical coordinates of the center of the cells
x = np.linspace(im_extent[0], im_extent[1], map_size) + 0.5 * dx x = np.linspace(im_extent[0], im_extent[1], map_size) + 0.5 * dx - pos_star[0]
y = np.linspace(im_extent[2], im_extent[3], map_size) + 0.5 * dy y = np.linspace(im_extent[2], im_extent[3], map_size) + 0.5 * dy - pos_star[1]
xx, yy = np.meshgrid(x, y) xx, yy = np.meshgrid(x, y)
rr = np.sqrt((xx - pos_star[0]) ** 2 + (yy - pos_star[1]) ** 2) rr = np.sqrt(xx ** 2 + yy ** 2)
# Rotational support # Rotational support
R = vphi ** 2 / rr R = vphi ** 2 / rr
# Equilibrium # Equilibrium
Gvrx, Gvry = np.gradient(vr) Gvry, Gvrx = np.gradient(vr, dy, dx, edge_order=2)
gradvr = (xx * Gvrx + yy * Gvry) / rr gradvr = (xx * Gvrx + yy * Gvry) / rr
dvr = gradvr + vr * gradvr # Complete derivative dvr = gradvr + vr * gradvr # Complete derivative
# Thermal support # Thermal support
GPx, GPy = np.gradient(Pz) GPy, GPx = np.gradient(Pz, dy, dx, edge_order=2)
gradPr = (xx * GPx + yy * GPy) / rr gradPr = (xx * GPx + yy * GPy) / rr
fP = gradPr / rho fP = gradPr / rho