Skip to content

Commit

Permalink
Psi4 Fix
Browse files Browse the repository at this point in the history
  • Loading branch information
@AlexCarpenter46
AlexCarpenter46 committed Jan 26, 2023
1 parent 19712dd commit 5469717
Show file tree
Hide file tree
Showing 3 changed files with 66 additions and 23 deletions.
63 changes: 50 additions & 13 deletions src/PointwiseFunctions/GeneralRelativity/Psi4.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,11 +28,13 @@ void psi_4(const gsl::not_null<Scalar<ComplexDataVector>*> psi_4_result,
const tnsr::ii<DataVector, 3, Frame>& spatial_metric,
const tnsr::II<DataVector, 3, Frame>& inverse_spatial_metric,
const tnsr::I<DataVector, 3, Frame>& inertial_coords) {
Variables<tmpl::list<::Tags::TempScalar<0>, ::Tags::TempI<0, 3, Frame>,
::Tags::TempI<1, 3, Frame>, ::Tags::Tempi<0, 3, Frame>,
::Tags::Tempij<0, 3, Frame>, ::Tags::Tempii<0, 3, Frame>,
::Tags::Tempii<1, 3, Frame>, ::Tags::TempIj<0, 3, Frame>,
::Tags::TempII<0, 3, Frame>>>
Variables<
tmpl::list<::Tags::TempScalar<0>, ::Tags::TempScalar<1>,
::Tags::TempI<0, 3, Frame>, ::Tags::TempI<1, 3, Frame>,
::Tags::TempI<2, 3, Frame>, ::Tags::TempI<3, 3, Frame>,
::Tags::Tempi<0, 3, Frame>, ::Tags::Tempij<0, 3, Frame>,
::Tags::Tempii<0, 3, Frame>, ::Tags::Tempii<1, 3, Frame>,
::Tags::TempIj<0, 3, Frame>, ::Tags::TempII<0, 3, Frame>>>
temp_buffer{get<0>(inertial_coords).size()};
auto& magnitude_cartesian = get<::Tags::TempScalar<0>>(temp_buffer);
magnitude(make_not_null(&magnitude_cartesian), inertial_coords,
Expand Down Expand Up @@ -69,27 +71,62 @@ void psi_4(const gsl::not_null<Scalar<ComplexDataVector>*> psi_4_result,
inverse_spatial_metric, cov_deriv_extrinsic_curvature, r_hat,
inverse_projection_tensor, projection_tensor, projection_up_lo, 1.0);

auto& x_coord = get<::Tags::TempI<0, 3, Frame>>(temp_buffer);
// Gram-Schmidt x_hat
auto& x_coord = get<::Tags::TempI<1, 3, Frame>>(temp_buffer);
x_coord.get(0) = 1.0;
x_coord.get(1) = x_coord.get(2) = 0.0;
auto& magnitude_x = get<::Tags::TempScalar<0>>(temp_buffer);
magnitude(make_not_null(&magnitude_x), x_coord, spatial_metric);
auto& x_hat = get<::Tags::TempI<0, 3, Frame>>(temp_buffer);
tenex::evaluate<ti::I>(make_not_null(&x_hat), x_coord(ti::I) / magnitude_x());
auto& x_component = get<::Tags::TempScalar<1>>(temp_buffer);
dot_product(make_not_null(&x_component), x_coord, r_hat, spatial_metric);
auto& x_hat = get<::Tags::TempI<2, 3, Frame>>(temp_buffer);
tenex::evaluate<ti::I>(make_not_null(&x_hat),
x_coord(ti::I) - (x_component() * r_hat(ti::I)));
auto& magnitude_x = get<::Tags::TempScalar<1>>(temp_buffer);
magnitude(make_not_null(&magnitude_x), x_hat, spatial_metric);
for (size_t j = 0; j < 3; j++) {
for (size_t i = 0; i < get(magnitude_x).size(); i++) {
if (magnitude_x.get()[i] != 0.0) {
x_hat.get(j)[i] = x_hat.get(j)[i] / magnitude_x.get()[i];
} else {
x_hat.get(j)[i] = 0.0;
}
}
}

Variables<tmpl::list<::Tags::TempI<0, 3, Frame, ComplexDataVector>,
::Tags::TempI<1, 3, Frame, ComplexDataVector>>>
y_hat_buffer{get<0>(inertial_coords).size()};

// Grad-Schmidt y_hat
auto& y_coord = get<::Tags::TempI<1, 3, Frame>>(temp_buffer);
y_coord.get(1) = 1.0;
y_coord.get(0) = y_coord.get(2) = 0.0;
auto& magnitude_y = get<::Tags::TempScalar<0>>(temp_buffer);
magnitude(make_not_null(&magnitude_y), y_coord, spatial_metric);
auto& y_component = get<::Tags::TempScalar<1>>(temp_buffer);
dot_product(make_not_null(&y_component), y_coord, r_hat, spatial_metric);
auto& y_hat_not_complex = get<::Tags::TempI<3, 3, Frame>>(temp_buffer);
tenex::evaluate<ti::I>(make_not_null(&y_hat_not_complex),
y_coord(ti::I) - (y_component() * r_hat(ti::I)));
dot_product(make_not_null(&y_component), y_coord, x_hat, spatial_metric);
tenex::evaluate<ti::I>(
make_not_null(&y_hat_not_complex),
y_hat_not_complex(ti::I) - (y_component() * x_hat(ti::I)));
auto& magnitude_y = get<::Tags::TempScalar<1>>(temp_buffer);
magnitude(make_not_null(&magnitude_y), y_hat_not_complex, spatial_metric);
for (size_t j = 0; j < 3; j++) {
for (size_t i = 0; i < get(magnitude_y).size(); i++) {
if (magnitude_y.get()[i] != 0.0) {
y_hat_not_complex.get(j)[i] =
y_hat_not_complex.get(j)[i] / magnitude_y.get()[i];
} else {
y_hat_not_complex.get(j)[i] = 0.0;
}
}
}
const std::complex<double> imag = std::complex<double>(0.0, 1.0);
auto& y_hat =
get<::Tags::TempI<0, 3, Frame, ComplexDataVector>>(y_hat_buffer);
tenex::evaluate<ti::I>(make_not_null(&y_hat),
imag * y_coord(ti::I) / magnitude_y());
imag * y_hat_not_complex(ti::I));

auto& m_bar =
get<::Tags::TempI<1, 3, Frame, ComplexDataVector>>(y_hat_buffer);
tenex::evaluate<ti::I>(make_not_null(&m_bar), x_hat(ti::I) - y_hat(ti::I));
Expand Down
23 changes: 16 additions & 7 deletions tests/Unit/PointwiseFunctions/GeneralRelativity/Psi4.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,15 +35,24 @@ def psi_4(spatial_ricci, extrinsic_curvature, cov_deriv_extrinsic_curvature,

x_coord = np.zeros((3))
x_coord[0] = 1
magnitude_x = math.sqrt(
np.einsum("a,b,ab", x_coord, x_coord, spatial_metric))
x_hat = x_coord / magnitude_x
x_component = np.einsum("a,b,ab", x_coord, r_hat, spatial_metric)
x_hat = x_coord - (x_component * r_hat)
magnitude_x = math.sqrt(np.einsum("a,b,ab", x_hat, x_hat, spatial_metric))
if (magnitude_x != 0.0):
x_hat = np.einsum("a", x_hat / magnitude_x)
else:
x_hat = np.einsum("a", x_hat * 0.0)
y_coord = np.zeros((3))
y_coord[1] = 1
magnitude_y = math.sqrt(
np.einsum("a,b,ab", y_coord, y_coord, spatial_metric))
y_hat = y_coord / magnitude_y

y_component = np.einsum("a,b,ab", y_coord, r_hat, spatial_metric)
y_hat = y_coord - (y_component * r_hat)
y_component = np.einsum("a,b,ab", y_coord, x_hat, spatial_metric)
y_hat = y_hat - (y_component * x_hat)
magnitude_y = math.sqrt(np.einsum("a,b,ab", y_hat, y_hat, spatial_metric))
if (magnitude_y != 0.0):
y_hat = np.einsum("a", y_hat / magnitude_y)
else:
y_hat = np.einsum("a", y_hat * 0.0)
m_bar = x_hat - (y_hat * complex(0.0, 1.0))

return (-0.5 * np.einsum("ab,a,b", u8_plus, m_bar, m_bar))
3 changes: 0 additions & 3 deletions tests/Unit/PointwiseFunctions/GeneralRelativity/Test_Psi4.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -96,9 +96,6 @@ void test_psi_4(const RealDataType& used_for_size_real,
auto inertial_coords =
make_with_random_values<tnsr::I<RealDataType, 3, Frame::Inertial>>(
nn_generator, nn_distribution, used_for_size_real);
for (size_t i = 0; i < 3; ++i) {
inertial_coords.get(i)[0] = 0.0;
}

const auto python_psi_4 = pypp::call<Scalar<ComplexDataType>>(
"GeneralRelativity.Psi4", "psi_4", spatial_ricci, extrinsic_curvature,
Expand Down

0 comments on commit 5469717

Please sign in to comment.