Skip to content

Commit 3f770dc

Browse files
EmilyBourneEmilyBourne
authoredFeb 24, 2025
Use new vector type in existing code (#85)
Use the new vector type that arises from #81 in the existing code.
1 parent 5399bc4 commit 3f770dc

19 files changed

+318
-98
lines changed
 

‎.github/actions/test_filter/action.yml

+1-1
Original file line numberDiff line numberDiff line change
@@ -32,7 +32,7 @@ runs:
3232
changed_files=$(awk '{printf("%s ",$0)} END { printf "\n" }' changed_files.txt)
3333
echo "CHANGED_FILES=${changed_files}" >> $GITHUB_ENV
3434
FOUND_POLAR_SPLINES=$(grep "src/interpolation/polar_splines/.*" -x changed_files.txt || true)
35-
if [ -z "$FOUND_POLAR_SPLINES" ]; then echo "POLAR_SPLINES_TEST_DEGREE_MIN=1" >> build.env; echo "POLAR_SPLINES_TEST_DEGREE_MIN=6" >> $GITHUB_ENV; fi
35+
if [ -n "$FOUND_POLAR_SPLINES" ]; then echo "POLAR_SPLINES_TEST_DEGREE_MIN=1" >> build.env; echo "POLAR_SPLINES_TEST_DEGREE_MIN=6" >> $GITHUB_ENV; fi
3636
FOUND_POISSON_2D=$(grep "src/geometryRTheta/poisson/.*" -x changed_files.txt || true)
3737
if [ -z "$FOUND_POISSON_2D" ]; then echo "POISSON_2D_BUILD_TESTING=OFF" >> $GITHUB_ENV; else echo "POISSON_2D_BUILD_TESTING=ON" >> $GITHUB_ENV; fi
3838
FOUND_DDC=$(grep "vendor/ddc" changed_files.txt || true)

‎src/data_types/tensor.hpp

+34-5
Original file line numberDiff line numberDiff line change
@@ -109,6 +109,22 @@ class Tensor
109109
"Filling the tensor on initialisation is only permitted for 1D vector objects");
110110
}
111111

112+
/**
113+
* @brief Construct a 1D tensor object from a coordinate.
114+
*
115+
* @param coord The coordinate.
116+
*/
117+
template <class... Dims>
118+
explicit KOKKOS_FUNCTION Tensor(Coord<Dims...> coord) : m_data(coord.array())
119+
{
120+
static_assert(
121+
rank() == 1,
122+
"Filling the tensor on initialisation is only permitted for 1D vector objects");
123+
static_assert(
124+
std::is_same_v<VectorIndexSet<Dims...>, ddc::type_seq_element_t<0, index_set>>,
125+
"The coordinate must have the same memory layout to make a clean conversion.");
126+
}
127+
112128
/**
113129
* @brief Construct a tensor object by copying an existing tensor.
114130
*
@@ -151,6 +167,24 @@ class Tensor
151167
*/
152168
KOKKOS_DEFAULTED_FUNCTION Tensor& operator=(Tensor const& other) = default;
153169

170+
/**
171+
* @brief A copy operator.
172+
* @param coord The coordinate to be copied into the vector.
173+
* @return A reference to the current tensor.
174+
*/
175+
template <class... Dims>
176+
KOKKOS_FUNCTION Tensor& operator=(Coord<Dims...> coord)
177+
{
178+
static_assert(
179+
rank() == 1,
180+
"Copying a coordinate into a tensor is only possible for 1D tensor objects");
181+
static_assert(
182+
std::is_same_v<VectorIndexSet<Dims...>, ddc::type_seq_element_t<0, index_set>>,
183+
"The coordinate must have the same memory layout to make a clean conversion.");
184+
m_data = coord.array();
185+
return *this;
186+
}
187+
154188
/**
155189
* @brief An operator to multiply all the element of the current tensor by
156190
* a value.
@@ -295,9 +329,6 @@ using to_tensor_t = typename detail::ToTensor<ElementType, TypeSeqValidIndexSet>
295329
template <class... ValidIndexSet>
296330
using DTensor = Tensor<double, ValidIndexSet...>;
297331

298-
// These objects should be extricated from the namespace in a later PR
299-
namespace tensor_tools {
300-
301332
/**
302333
* @brief A helper type alias to get a 1D tensor (a vector).
303334
* @tparam ElementType The type of the elements of the tensor (usually double/complex).
@@ -313,8 +344,6 @@ using Vector = Tensor<ElementType, VectorIndexSet<Dims...>>;
313344
template <class... Dims>
314345
using DVector = Vector<double, Dims...>;
315346

316-
} // namespace tensor_tools
317-
318347
//////////////////////////////////////////////////////////////////////////
319348
// Operators
320349
//////////////////////////////////////////////////////////////////////////

‎src/data_types/vector_field_common.hpp

+2
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@
44
#include <ddc/ddc.hpp>
55

66
#include "ddc_aliases.hpp"
7+
#include "tensor.hpp"
78

89
template <class FieldType, class NDTypeSeq>
910
class VectorFieldCommon;
@@ -53,6 +54,7 @@ inline constexpr auto get(VectorFieldType const& field) noexcept
5354
template <class QueryTag, class VectorFieldType>
5455
inline constexpr auto get(VectorFieldType& field) noexcept
5556
{
57+
static_assert(is_vector_field_v<VectorFieldType>);
5658
return field.template get<QueryTag>();
5759
}
5860

‎src/geometryRTheta/geometry/geometry.hpp

+18
Original file line numberDiff line numberDiff line change
@@ -381,6 +381,15 @@ struct X
381381
* Here, not periodic.
382382
*/
383383
static bool constexpr PERIODIC = false;
384+
385+
/// A boolean indicating if dimension describes a covariant coordinate.
386+
static bool constexpr IS_COVARIANT = true;
387+
388+
/// A boolean indicating if dimension describes a contravariant coordinate.
389+
static bool constexpr IS_CONTRAVARIANT = true;
390+
391+
/// A type-alias mapping to the covariant counterpart.
392+
using Dual = X;
384393
};
385394
/**
386395
* @brief Define non periodic real Y dimension.
@@ -392,6 +401,15 @@ struct Y
392401
* Here, not periodic.
393402
*/
394403
static bool constexpr PERIODIC = false;
404+
405+
/// A boolean indicating if dimension describes a covariant coordinate.
406+
static bool constexpr IS_COVARIANT = true;
407+
408+
/// A boolean indicating if dimension describes a contravariant coordinate.
409+
static bool constexpr IS_CONTRAVARIANT = true;
410+
411+
/// A type-alias mapping to the covariant counterpart.
412+
using Dual = Y;
395413
};
396414

397415
/**

‎src/mapping/geometry_pseudo_cartesian.hpp

+19
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,17 @@ struct X_pC
1616
* Here, not periodic.
1717
*/
1818
static bool constexpr PERIODIC = false;
19+
20+
/// A boolean indicating if dimension describes a covariant coordinate.
21+
static bool constexpr IS_COVARIANT = true;
22+
23+
/// A boolean indicating if dimension describes a contravariant coordinate.
24+
static bool constexpr IS_CONTRAVARIANT = true;
25+
26+
/// A type-alias mapping to the co/contra-variant counterpart.
27+
using Dual = X_pC;
1928
};
29+
2030
/**
2131
* @brief Tag the second non periodic dimension
2232
* in the pseudo_Cartesian index range.
@@ -28,6 +38,15 @@ struct Y_pC
2838
* Here, not periodic.
2939
*/
3040
static bool constexpr PERIODIC = false;
41+
42+
/// A boolean indicating if dimension describes a covariant coordinate.
43+
static bool constexpr IS_COVARIANT = true;
44+
45+
/// A boolean indicating if dimension describes a contravariant coordinate.
46+
static bool constexpr IS_CONTRAVARIANT = true;
47+
48+
/// A type-alias mapping to the co/contra-variant counterpart.
49+
using Dual = Y_pC;
3150
};
3251

3352
using CoordX_pC = Coord<X_pC>;

‎src/mapping/vector_mapper.hpp

+7-8
Original file line numberDiff line numberDiff line change
@@ -37,11 +37,6 @@ class VectorMapper<NDTag<XIn, YIn>, NDTag<XOut, YOut>, Mapping, ExecSpace>
3737
/// @brief The type of the memory space where the field is saved (CPU vs GPU).
3838
using memory_space = typename ExecSpace::memory_space;
3939

40-
/// The vector type in the coordinate system taken as input.
41-
using vector_element_type_in = DVector<XIn, YIn>;
42-
/// The vector type in the coordinate system returned as output.
43-
using vector_element_type_out = DVector<XOut, YOut>;
44-
4540
private:
4641
Mapping m_mapping;
4742

@@ -84,7 +79,7 @@ class VectorMapper<NDTag<XIn, YIn>, NDTag<XOut, YOut>, Mapping, ExecSpace>
8479
Matrix_2x2 map_J;
8580
mapping_proxy.jacobian_matrix(ddc::coordinate(idx), map_J);
8681

87-
vector_element_type_out vector_out;
82+
Coord<XOut, YOut> vector_out;
8883
vector_out.array() = mat_vec_mul(map_J, vector_field_input(idx).array());
8984
ddcHelper::get<XOut>(vector_field_output)(idx) = ddc::get<XOut>(vector_out);
9085
ddcHelper::get<YOut>(vector_field_output)(idx) = ddc::get<YOut>(vector_out);
@@ -97,8 +92,12 @@ class VectorMapper<NDTag<XIn, YIn>, NDTag<XOut, YOut>, Mapping, ExecSpace>
9792
KOKKOS_LAMBDA(IdxType idx) {
9893
Matrix_2x2 map_J = inv_mapping(ddc::coordinate(idx));
9994

100-
vector_element_type_out vector_out;
101-
vector_out.array() = mat_vec_mul(map_J, vector_field_input(idx).array());
95+
Coord<XOut, YOut> vector_out;
96+
// mat_vec_mul should be replaced with a tensor calculus function
97+
// when map_J is stored in a Tensor
98+
vector_out.array() = mat_vec_mul(
99+
map_J,
100+
ddcHelper::to_coord(vector_field_input(idx)).array());
102101
ddcHelper::get<XOut>(vector_field_output)(idx) = ddc::get<XOut>(vector_out);
103102
ddcHelper::get<YOut>(vector_field_output)(idx) = ddc::get<YOut>(vector_out);
104103
});

‎src/math_tools/l_norm_tools.hpp

+22
Original file line numberDiff line numberDiff line change
@@ -30,6 +30,28 @@ KOKKOS_FUNCTION double norm_inf(ddc::Coordinate<Tags...> coord)
3030
return result;
3131
}
3232

33+
/**
34+
* @brief Compute the infinity norm of a vector on an
35+
* orthonormal coordinate system.
36+
*
37+
* For a given vector @f$ x @f$ , compute
38+
* @f$|Vert x |Vert_{\infty} = \sup_n |x_n| @f$.
39+
*
40+
* @param[in] vec The given vector.
41+
*
42+
* @return A double containing the value of the infinty norm.
43+
*/
44+
template <class... Tags>
45+
KOKKOS_FUNCTION double norm_inf(DVector<Tags...> vec)
46+
{
47+
using index_set = typename DVector<Tags...>::vector_index_set_t<0>;
48+
static_assert(
49+
std::is_same_v<index_set, vector_index_set_dual_t<index_set>>,
50+
"Mapping is needed to calculate norm_inf on a non-orthonormal coordinate system");
51+
double result = 0.0;
52+
((result = Kokkos::max(result, Kokkos::fabs(ddcHelper::get<Tags>(vec)))), ...);
53+
return result;
54+
}
3355

3456
/**
3557
* @brief Compute the infinity norm.

‎src/timestepper/itimestepper.hpp

+5-2
Original file line numberDiff line numberDiff line change
@@ -161,12 +161,15 @@ class ITimeStepper
161161
* @param[in] new_val The coordinate that should be saved to the vector field.
162162
*/
163163
template <class DerivFieldType, class Idx, class... DDims>
164-
KOKKOS_FUNCTION static void fill_k_total(DerivFieldType k_total, Idx i, Coord<DDims...> new_val)
164+
KOKKOS_FUNCTION static void fill_k_total(
165+
DerivFieldType k_total,
166+
Idx i,
167+
DVector<DDims...> new_val)
165168
{
166169
static_assert(
167170
(std::is_same_v<DerivField, DerivFieldType>)
168171
|| (is_multipatch_field_v<DerivField>));
169-
((ddcHelper::get<DDims>(k_total)(i) = ddc::get<DDims>(new_val)), ...);
172+
((ddcHelper::get<DDims>(k_total)(i) = ddcHelper::get<DDims>(new_val)), ...);
170173
}
171174

172175
/**

‎src/utils/ddc_aliases.hpp

-24
Original file line numberDiff line numberDiff line change
@@ -77,30 +77,6 @@ using UniformGridBase = ddc::UniformPointSampling<GridType>;
7777
template <class GridType>
7878
using NonUniformGridBase = ddc::NonUniformPointSampling<GridType>;
7979

80-
/// A type describing a vector e.g. (E_x, E_y)
81-
template <class ElementType, class... Dims>
82-
class Vector : public ddc::detail::TaggedVector<ElementType, Dims...>
83-
{
84-
public:
85-
/**
86-
* @brief A constructor for a Vector to build from individual elements:
87-
* E.g. Vector<double, X, Y>(1.0, 2.0)
88-
*
89-
* or to build from component vectors
90-
* E.g. Vector<double, X, Y>(vec_x, vec_y)
91-
*
92-
* @param[in] p The element of the Vector.
93-
*/
94-
template <class... Params>
95-
KOKKOS_FUNCTION Vector(Params... p) : ddc::detail::TaggedVector<ElementType, Dims...>(p...)
96-
{
97-
}
98-
};
99-
100-
/// A type describing a vector whose elements are doubles e.g. (E_x, E_y)
101-
template <class... Dims>
102-
using DVector = Vector<double, Dims...>;
103-
10480

10581
//----------------------------------------------
10682
// Class for batched 1D spline builder

0 commit comments

Comments
 (0)