modify derivative of inv jacobian; test left?

jyoo1042 · Aug 28, 2024 · c574681 · c574681
1 parent e038674
commit c574681
Showing 1 changed file with 32 additions and 30 deletions.
diff --git a/src/PointwiseFunctions/AnalyticData/GrMhd/SphericalTorus.cpp b/src/PointwiseFunctions/AnalyticData/GrMhd/SphericalTorus.cpp
@@ -326,8 +326,10 @@ tnsr::Ijk<T, 3, Frame::NoFrame> SphericalTorus::derivative_of_inv_jacobian(
 
   auto& cos_theta = get<2, 2, 1>(result);
   auto& sin_theta = get<2, 2, 2>(result);
+  get<0, 0, 0>(result) = compression_level_ * pow<3>(get<1>(source_coords)) +
+                         (1.0 - compression_level_) * (get<1>(source_coords));
   get<0, 0, 0>(result) =
-      M_PI_2 - (pi_over_2_minus_theta_min_ * get<1>(source_coords));
+      M_PI_2 - (pi_over_2_minus_theta_min_ * get<0, 0, 0>(result));
   cos_theta = cos(get<0, 0, 0>(result));
   sin_theta = sin(get<0, 0, 0>(result));
 
@@ -344,45 +346,45 @@ tnsr::Ijk<T, 3, Frame::NoFrame> SphericalTorus::derivative_of_inv_jacobian(
   radius(make_not_null(&r), get<0>(source_coords));
   r_factor = 0.5 * (r_max_ - r_min_);
 
-  get<0, 0, 1>(result) = -theta_factor / r_factor * cos_theta * cos_phi;
-  // *= eta_prime
+  get<0, 0, 1>(result) =
+      -theta_factor / r_factor * cos_theta * cos_phi * eta_prime;
   get<0, 0, 2>(result) = -phi_factor / r_factor * sin_theta * sin_phi;
-  get<0, 1, 1>(result) = -theta_factor / r_factor * cos_theta * sin_phi;
-  // *= eta_prime
+  get<0, 1, 1>(result) =
+      -theta_factor / r_factor * cos_theta * sin_phi * eta_prime;
   get<0, 1, 2>(result) = phi_factor / r_factor * sin_theta * cos_phi;
-  get<0, 2, 1>(result) = theta_factor / r_factor * sin_theta;
-  // *= eta_prime
+  get<0, 2, 1>(result) = theta_factor / r_factor * sin_theta * eta_prime;
   get<1, 0, 0>(result) =
-      r_factor / theta_factor * cos_theta * cos_phi / square(r);
-  // *= (eta_prime)^-1
-  get<1, 0, 1>(result) = -sin_theta * cos_phi / r;
-  // copy over notes
-
-  get<1, 0, 2>(result) = phi_factor / theta_factor * cos_theta * sin_phi / r;
-  // *= (eta_prime)^-1
+      r_factor / theta_factor * cos_theta * cos_phi / (square(r) * eta_prime);
+  get<1, 0, 1>(result) =
+      (-cos_phi / (r * theta_factor)) *
+      (sin_theta * theta_factor * eta_prime + cos_theta * eta_double_prime) /
+      (square(eta_prime));
+  get<1, 0, 2>(result) =
+      phi_factor / theta_factor * cos_theta * sin_phi / (r * eta_prime);
   get<1, 1, 0>(result) =
-      r_factor / theta_factor * cos_theta * sin_phi / square(r);
-  // *= (eta_prime)^-1
-  get<1, 1, 1>(result) = -sin_theta * sin_phi / r;
-  // copy over notes
-
-  get<1, 1, 2>(result) = -phi_factor / theta_factor * cos_theta * cos_phi / r;
-  // *= (eta_prime)^-1
-  get<1, 2, 0>(result) = -r_factor / theta_factor * sin_theta / square(r);
-  // *= (eta_prime)^-1
-  get<1, 2, 1>(result) = -cos_theta / r;
+      r_factor / theta_factor * cos_theta * sin_phi / (square(r) * eta_prime);
+  get<1, 1, 1>(result) =
+      (-sin_phi / (r * theta_factor)) *
+      (sin_theta * theta_factor * eta_prime + cos_theta * eta_double_prime) /
+      (square(eta_prime));
+  get<1, 1, 2>(result) =
+      -phi_factor / theta_factor * cos_theta * cos_phi / (r * eta_prime);
+  get<1, 2, 0>(result) =
+      -r_factor / theta_factor * sin_theta / (square(r) * eta_prime);
+  get<1, 2, 1>(result) =
+      (-1.0 / (theta_factor * r)) *
+      (cos_theta * theta_factor + sin_theta * eta_double_prime) /
+      (square(eta_prime));
   // copy over notes
   get<2, 0, 0>(result) =
       r_factor / phi_factor * sin_phi / (square(r) * sin_theta);
-  get<2, 0, 1>(result) = -theta_factor / phi_factor * cos_theta * sin_phi /
-                         (r * square(sin_theta));
-  // *= eta_prime
+  get<2, 0, 1>(result) = -theta_factor / phi_factor * cos_theta * sin_phi *
+                         eta_prime / (r * square(sin_theta));
   get<2, 0, 2>(result) = -cos_phi / (r * sin_theta);
   get<2, 1, 0>(result) =
       -r_factor / phi_factor * cos_phi / (square(r) * sin_theta);
-  get<2, 1, 1>(result) =
-      theta_factor / phi_factor * cos_theta * cos_phi / (r * square(sin_theta));
-  // *= eta_prime
+  get<2, 1, 1>(result) = theta_factor / phi_factor * cos_theta * cos_phi *
+                         eta_prime / (r * square(sin_theta));
   get<2, 1, 2>(result) = -sin_phi / (r * sin_theta);
 
   // remove temp vars and restore to zero