Skip to content

Commit 0f774e9

Browse files
lsawadelsawade
authored
Merge pull request #1128 from PrincetonUniversity/issue-990-assembly-sources
Issue 990 assembly sources - Implement assembly sources
2 parents 94b5e97 + 5d2bc4f commit 0f774e9

26 files changed

+1847
-182
lines changed

‎core/specfem/assembly/compute_source_array.hpp‎

Lines changed: 27 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@
77
#include "specfem/assembly/mesh.hpp"
88
#include "specfem/source.hpp"
99
#include "specfem_setup.hpp"
10+
#include <Kokkos_Core.hpp>
1011

1112
namespace specfem::assembly {
1213

@@ -25,12 +26,37 @@ namespace specfem::assembly {
2526
* @param source_array The output source array to be filled.
2627
*
2728
*/
29+
template <typename SourceArrayViewType>
2830
void compute_source_array(
2931
const std::shared_ptr<
3032
specfem::sources::source<specfem::dimension::type::dim2> > &source,
3133
const specfem::assembly::mesh<specfem::dimension::type::dim2> &mesh,
3234
const specfem::assembly::jacobian_matrix<specfem::dimension::type::dim2>
3335
&jacobian_matrix,
34-
specfem::kokkos::HostView3d<type_real> source_array);
36+
SourceArrayViewType &source_array);
37+
38+
/**
39+
* @brief Compute the lagrange interpolants for a specific 3d source location
40+
* in a 3d element.
41+
*
42+
* This is a helper function that dispatches the computation to the appropriate
43+
* implementation based on the source type.
44+
*
45+
* @param source The source for which the source array is computed.
46+
* @param mesh The mesh on which the source is defined.
47+
* @param jacobian_matrix The Jacobian matrix for the mesh.
48+
* @param source_array The output source array to be filled.
49+
*/
50+
template <typename SourceArrayViewType>
51+
void compute_source_array(
52+
const std::shared_ptr<
53+
specfem::sources::source<specfem::dimension::type::dim3> > &source,
54+
const specfem::assembly::mesh<specfem::dimension::type::dim3> &mesh,
55+
const specfem::assembly::jacobian_matrix<specfem::dimension::type::dim3>
56+
&jacobian_matrix,
57+
SourceArrayViewType &source_array);
3558

3659
} // namespace specfem::assembly
60+
61+
#include "specfem/assembly/compute_source_array/dim2/compute_source_array.tpp"
62+
#include "specfem/assembly/compute_source_array/dim3/compute_source_array.tpp"

‎core/specfem/assembly/compute_source_array/dim2/compute_source_array.cpp‎ renamed to ‎core/specfem/assembly/compute_source_array/dim2/compute_source_array.tpp‎

Lines changed: 10 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,5 @@
1+
#pragma once
2+
13
#include "specfem/assembly/compute_source_array.hpp"
24
#include "enumerations/dimension.hpp"
35
#include "impl/compute_source_array_from_tensor.hpp"
@@ -6,14 +8,19 @@
68
#include "specfem/assembly/element_types.hpp"
79
#include "specfem/source.hpp"
810
#include "specfem_setup.hpp"
11+
#include <Kokkos_Core.hpp>
912

13+
template<typename SourceArrayViewType>
1014
void specfem::assembly::compute_source_array(
1115
const std::shared_ptr<
1216
specfem::sources::source<specfem::dimension::type::dim2> > &source,
1317
const specfem::assembly::mesh<specfem::dimension::type::dim2> &mesh,
1418
const specfem::assembly::jacobian_matrix<specfem::dimension::type::dim2>
1519
&jacobian_matrix,
16-
specfem::kokkos::HostView3d<type_real> source_array) {
20+
SourceArrayViewType &source_array) {
21+
22+
// Ensure source_array is a 3D view
23+
static_assert(SourceArrayViewType::rank == 3, "Source array must be in rank 3.");
1724

1825
switch (source->get_source_type()) {
1926
case specfem::sources::source_type::vector_source: {
@@ -27,9 +34,10 @@ void specfem::assembly::compute_source_array(
2734
"Source is not of vector type. Cannot compute vector source "
2835
"array.");
2936
}
30-
37+
std::cout << "computing for vector" << std::endl;
3138
specfem::assembly::compute_source_array_impl::from_vector(*vector_source,
3239
source_array);
40+
std::cout << "computed for vector: " << std::endl;
3341
break;
3442
}
3543
case specfem::sources::source_type::tensor_source: {

‎core/specfem/assembly/compute_source_array/dim2/impl/compute_source_array_from_tensor.cpp‎

Lines changed: 5 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,7 @@
99
#include "specfem/point.hpp"
1010
#include "specfem/source.hpp"
1111
#include "specfem_setup.hpp"
12+
#include <Kokkos_Core.hpp>
1213

1314
// Local namespace for implementation details
1415
namespace specfem::assembly::compute_source_array_impl {
@@ -19,7 +20,8 @@ void compute_source_array_from_tensor_and_element_jacobian(
1920
const JacobianViewType &element_jacobian_matrix,
2021
const specfem::assembly::mesh_impl::quadrature<
2122
specfem::dimension::type::dim2> &quadrature,
22-
specfem::kokkos::HostView3d<type_real> source_array) {
23+
Kokkos::View<type_real ***, Kokkos::LayoutRight, Kokkos::HostSpace>
24+
source_array) {
2325

2426
const int ngllx = quadrature.N;
2527
const int ngllz = quadrature.N;
@@ -90,7 +92,8 @@ void specfem::assembly::compute_source_array_impl::from_tensor(
9092
const specfem::assembly::mesh<specfem::dimension::type::dim2> &mesh,
9193
const specfem::assembly::jacobian_matrix<specfem::dimension::type::dim2>
9294
&jacobian_matrix,
93-
specfem::kokkos::HostView3d<type_real> source_array) {
95+
Kokkos::View<type_real ***, Kokkos::LayoutRight, Kokkos::HostSpace>
96+
source_array) {
9497

9598
const int ngllx = source_array.extent(2);
9699
const int ngllz = source_array.extent(1);

‎core/specfem/assembly/compute_source_array/dim2/impl/compute_source_array_from_tensor.hpp‎

Lines changed: 5 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@
33
#include "specfem/assembly/jacobian_matrix.hpp"
44
#include "specfem/assembly/mesh.hpp"
55
#include "specfem/source.hpp"
6+
#include <Kokkos_Core.hpp>
67

78
namespace specfem::assembly::compute_source_array_impl {
89

@@ -17,7 +18,8 @@ void from_tensor(
1718
const specfem::assembly::mesh<specfem::dimension::type::dim2> &mesh,
1819
const specfem::assembly::jacobian_matrix<specfem::dimension::type::dim2>
1920
&jacobian_matrix,
20-
specfem::kokkos::HostView3d<type_real> source_array);
21+
Kokkos::View<type_real ***, Kokkos::LayoutRight, Kokkos::HostSpace>
22+
source_array);
2123

2224
// Implementation details for computing source array from tensor and element
2325
// jacobian
@@ -27,6 +29,7 @@ void compute_source_array_from_tensor_and_element_jacobian(
2729
const JacobianViewType &element_jacobian_matrix,
2830
const specfem::assembly::mesh_impl::quadrature<
2931
specfem::dimension::type::dim2> &quadrature,
30-
specfem::kokkos::HostView3d<type_real> source_array);
32+
Kokkos::View<type_real ***, Kokkos::LayoutRight, Kokkos::HostSpace>
33+
source_array);
3134

3235
} // namespace specfem::assembly::compute_source_array_impl

‎core/specfem/assembly/compute_source_array/dim2/impl/compute_source_array_from_vector.cpp‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -9,11 +9,13 @@
99
#include "specfem/point.hpp"
1010
#include "specfem/source.hpp"
1111
#include "specfem_setup.hpp"
12+
#include <Kokkos_Core.hpp>
1213

1314
void specfem::assembly::compute_source_array_impl::from_vector(
1415
const specfem::sources::vector_source<specfem::dimension::type::dim2>
1516
&vector_source,
16-
specfem::kokkos::HostView3d<type_real> source_array) {
17+
Kokkos::View<type_real ***, Kokkos::LayoutRight, Kokkos::HostSpace>
18+
source_array) {
1719

1820
const int ngllx = source_array.extent(2);
1921
const int ngllz = source_array.extent(1);

‎core/specfem/assembly/compute_source_array/dim2/impl/compute_source_array_from_vector.hpp‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2,12 +2,14 @@
22
#include "kokkos_abstractions.h"
33
#include "specfem/assembly/mesh.hpp"
44
#include "specfem/source.hpp"
5+
#include <Kokkos_Core.hpp>
56

67
namespace specfem::assembly::compute_source_array_impl {
78

89
void from_vector(
910
const specfem::sources::vector_source<specfem::dimension::type::dim2>
1011
&source,
11-
specfem::kokkos::HostView3d<type_real> source_array);
12+
Kokkos::View<type_real ***, Kokkos::LayoutRight, Kokkos::HostSpace>
13+
source_array);
1214

1315
} // namespace specfem::assembly::compute_source_array_impl

‎core/specfem/assembly/compute_source_array/dim3/compute_source_array.tpp‎

Lines changed: 65 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,65 @@
1+
#pragma once
2+
3+
#include "specfem/assembly/compute_source_array.hpp"
4+
#include "enumerations/dimension.hpp"
5+
#include "impl/compute_source_array_from_tensor.hpp"
6+
#include "impl/compute_source_array_from_vector.hpp"
7+
#include "kokkos_abstractions.h"
8+
#include "specfem/assembly/element_types.hpp"
9+
#include "specfem/source.hpp"
10+
#include "specfem_setup.hpp"
11+
#include <Kokkos_Core.hpp>
12+
13+
template<typename SourceArrayViewType>
14+
void specfem::assembly::compute_source_array(
15+
const std::shared_ptr<
16+
specfem::sources::source<specfem::dimension::type::dim3> > &source,
17+
const specfem::assembly::mesh<specfem::dimension::type::dim3> &mesh,
18+
const specfem::assembly::jacobian_matrix<specfem::dimension::type::dim3>
19+
&jacobian_matrix,
20+
SourceArrayViewType &source_array) {
21+
22+
// Ensure source_array is a 4D view
23+
static_assert(SourceArrayViewType::rank == 4, "Source array must be in rank 4.");
24+
25+
26+
switch (source->get_source_type()) {
27+
case specfem::sources::source_type::vector_source: {
28+
29+
// Cast to derived class to access specific methods
30+
auto vector_source = static_cast<const specfem::sources::vector_source<
31+
specfem::dimension::type::dim3> *>(source.get());
32+
33+
if (!vector_source) {
34+
KOKKOS_ABORT_WITH_LOCATION(
35+
"Source is not of vector type. Cannot compute vector source "
36+
"array.");
37+
}
38+
39+
specfem::assembly::compute_source_array_impl::from_vector(*vector_source,
40+
source_array);
41+
break;
42+
}
43+
// case specfem::sources::source_type::tensor_source: {
44+
45+
// // Cast to derived class to access specific methods
46+
// auto tensor_source = static_cast<const specfem::sources::tensor_source<
47+
// specfem::dimension::type::dim3> *>(source.get());
48+
49+
// if (!tensor_source) {
50+
// KOKKOS_ABORT_WITH_LOCATION(
51+
// "Source is not of tensor type. Cannot compute tensor source "
52+
// "array.");
53+
// }
54+
55+
// specfem::assembly::compute_source_array_impl::from_tensor(
56+
// *tensor_source, mesh, jacobian_matrix, source_array);
57+
// break;
58+
// }
59+
default:
60+
// Handle unsupported source types
61+
KOKKOS_ABORT_WITH_LOCATION(
62+
"Unsupported source type for compute_source_array.");
63+
}
64+
return;
65+
}

‎core/specfem/assembly/compute_source_array/dim3/impl/compute_source_array_from_tensor.cpp‎

Lines changed: 152 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,152 @@
1+
#include "compute_source_array_from_tensor.hpp"
2+
#include "algorithms/interface.hpp"
3+
#include "enumerations/macros.hpp"
4+
#include "kokkos_abstractions.h"
5+
#include "quadrature/interface.hpp"
6+
#include "specfem/assembly/element_types.hpp"
7+
#include "specfem/assembly/jacobian_matrix.hpp"
8+
#include "specfem/assembly/mesh.hpp"
9+
#include "specfem/point.hpp"
10+
#include "specfem/source.hpp"
11+
#include "specfem_setup.hpp"
12+
13+
/*
14+
// Local namespace for implementation details
15+
namespace specfem::assembly::compute_source_array_impl {
16+
17+
void compute_source_array_from_tensor_and_element_jacobian(
18+
const specfem::sources::tensor_source<specfem::dimension::type::dim3>
19+
&tensor_source,
20+
const JacobianViewType &element_jacobian_matrix,
21+
const specfem::assembly::mesh_impl::quadrature<
22+
specfem::dimension::type::dim3> &quadrature,
23+
specfem::kokkos::HostView3d<type_real> source_array) {
24+
25+
const int ngllx = quadrature.N;
26+
const int nglly = quadrature.N;
27+
const int ngllz = quadrature.N;
28+
29+
// Create quadrature and compute xi/gamma arrays
30+
auto xi = quadrature.h_xi;
31+
auto eta = quadrature.h_xi;
32+
auto gamma = quadrature.h_xi;
33+
34+
// Compute lagrange interpolants at the local source location
35+
auto [hxi_source, hpxi_source] =
36+
specfem::quadrature::gll::Lagrange::compute_lagrange_interpolants(
37+
tensor_source.get_local_coordinates().xi, ngllx, xi);
38+
auto [heta_source, hpeta_source] =
39+
specfem::quadrature::gll::Lagrange::compute_lagrange_interpolants(
40+
tensor_source.get_local_coordinates().eta, nglly, eta);
41+
auto [hgamma_source, hpgamma_source] =
42+
specfem::quadrature::gll::Lagrange::compute_lagrange_interpolants(
43+
tensor_source.get_local_coordinates().gamma, ngllz, gamma);
44+
45+
specfem::kokkos::HostView3d<type_real> source_polynomial("source_polynomial",
46+
ngllz, nglly, ngllx);
47+
48+
// Use pre-computed jacobian data instead of loading from jacobian_matrix
49+
for (int iz = 0; iz < ngllz; ++iz) {
50+
for (int iy = 0; iy < nglly; ++iy) {
51+
for (int ix = 0; ix < ngllx; ++ix) {
52+
type_real hlagrange =
53+
hxi_source(ix) * heta_source(iy) * hgamma_source(iz);
54+
source_polynomial(iz, iy, ix) = hlagrange;
55+
}
56+
}
57+
}
58+
59+
// Store the derivatives in a function object for interpolation
60+
auto derivatives_source = specfem::algorithms::interpolate_function(
61+
source_polynomial, element_jacobian_matrix);
62+
63+
const auto source_tensor = tensor_source.get_source_tensor();
64+
65+
int ncomponents = source_array.extent(0);
66+
67+
// Sanity check
68+
if (ncomponents != source_tensor.extent(0)) {
69+
KOKKOS_ABORT_WITH_LOCATION(
70+
"source_array_components and tensor components do not match")
71+
}
72+
73+
for (int iz = 0; iz < ngllz; ++iz) {
74+
for (int iy = 0; iy < nglly; ++iy) {
75+
for (int ix = 0; ix < ngllx; ++ix) {
76+
77+
// Compute the derivatives at the source location
78+
type_real dsrc_dx =
79+
(hpxi_source(ix) * derivatives_source.xix) * heta_source(iy) *
80+
hgamma_source(iz) +
81+
hxi_source(ix) * (hpeta_source * derivatives_source.etax) *
82+
hgamma_source(iz) +
83+
hxi_source(ix) * heta_source(iy) *
84+
(hpgamma_source(iz) * derivatives_source.gammax);
85+
type_real dsrc_dy =
86+
(hpxi_source(ix) * derivatives_source.xiy) * heta_source(iy) *
87+
hgamma_source(iz) +
88+
hxi_source(ix) * (hpeta_source * derivatives_source.etay) *
89+
hgamma_source(iz) +
90+
hxi_source(ix) * heta_source(iy) *
91+
(hpgamma_source(iz) * derivatives_source.gammay);
92+
type_real dsrc_dz =
93+
(hpxi_source(ix) * derivatives_source.xiz) * heta_source(iy) *
94+
hgamma_source(iz) +
95+
hxi_source(ix) * (hpeta_source * derivatives_source.etaz) *
96+
hgamma_source(iz) +
97+
hxi_source(ix) * heta_source(iy) *
98+
(hpgamma_source(iz) * derivatives_source.gammaz);
99+
100+
for (int i = 0; i < ncomponents; ++i) {
101+
source_array(i, iz, ix) = source_tensor(i, 0) * dsrc_dx +
102+
source_tensor(i, 1) * dsrc_dy +
103+
source_tensor(i, 2) * dsrc_dz;
104+
}
105+
}
106+
}
107+
108+
return;
109+
}
110+
111+
} // namespace specfem::assembly::compute_source_array_impl
112+
113+
void specfem::assembly::compute_source_array_impl::from_tensor(
114+
const specfem::sources::tensor_source<specfem::dimension::type::dim3>
115+
&tensor_source,
116+
const specfem::assembly::mesh<specfem::dimension::type::dim3> &mesh,
117+
const specfem::assembly::jacobian_matrix<specfem::dimension::type::dim3>
118+
&jacobian_matrix,
119+
specfem::kokkos::HostView4d<type_real> source_array) {
120+
121+
const int ngllz = source_array.extent(1);
122+
const int nglly = source_array.extent(2);
123+
const int ngllx = source_array.extent(3);
124+
125+
using PointJacobianMatrix =
126+
specfem::point::jacobian_matrix<specfem::dimension::type::dim3, false,
127+
false>;
128+
Kokkos::View<PointJacobianMatrix ***, Kokkos::LayoutRight, Kokkos::HostSpace>
129+
element_jacobian( "element_jacobian", ngllz, nglly, ngllx);
130+
131+
// Extract jacobian data from jacobian_matrix
132+
for (int iz = 0; iz < ngllz; ++iz) {
133+
for (int iy = 0; iy < nglly; ++iy) {
134+
for (int ix = 0; ix < ngllx; ++ix) {
135+
const specfem::point::index<specfem::dimension::type::dim3> index(
136+
tensor_source.get_local_coordinates().ispec, iz, iy, ix);
137+
PointJacobianMatrix derivatives;
138+
specfem::assembly::load_on_host(index, jacobian_matrix, derivatives);
139+
element_jacobian(iz, iy, ix) = derivatives;
140+
}
141+
}
142+
}
143+
const auto &quadrature =
144+
static_cast<const specfem::assembly::mesh_impl::quadrature<
145+
specfem::dimension::type::dim3> &>(mesh);
146+
147+
specfem::assembly::compute_source_array_impl::
148+
compute_source_array_from_tensor_and_element_jacobian(
149+
tensor_source, element_jacobian, quadrature, source_array);
150+
return;
151+
}
152+
*/

0 commit comments

Comments
 (0)