From 75aff0086abefda47160b295b8a474812a7680cc Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Wed, 27 Sep 2023 16:32:19 +0200 Subject: [PATCH] Add single precision support to fvm::Nabla --- src/atlas/numerics/fvm/Nabla.cc | 665 +++++++++++++++------------ src/tests/numerics/test_fvm_nabla.cc | 435 ++++++++++-------- 2 files changed, 611 insertions(+), 489 deletions(-) diff --git a/src/atlas/numerics/fvm/Nabla.cc b/src/atlas/numerics/fvm/Nabla.cc index 91a3395fb..cd432e14e 100644 --- a/src/atlas/numerics/fvm/Nabla.cc +++ b/src/atlas/numerics/fvm/Nabla.cc @@ -95,76 +95,95 @@ void Nabla::gradient(const Field& field, Field& grad_field) const { void Nabla::gradient_of_scalar(const Field& scalar_field, Field& grad_field) const { Log::debug() << "Compute gradient of scalar field " << scalar_field.name() << " with fvm method" << std::endl; - const double radius = fvm_->radius(); - const double deg2rad = M_PI / 180.; - const mesh::Edges& edges = fvm_->mesh().edges(); - const mesh::Nodes& nodes = fvm_->mesh().nodes(); + auto dispatch = [&](auto value) { + using Value = std::decay_t; - const idx_t nnodes = fvm_->node_columns().nb_nodes(); - const idx_t nedges = fvm_->edge_columns().nb_edges(); + const Value radius = fvm_->radius(); + const Value deg2rad = M_PI / 180.; - const auto scalar = scalar_field.levels() - ? array::make_view(scalar_field).slice(Range::all(), Range::all()) - : array::make_view(scalar_field).slice(Range::all(), Range::dummy()); - auto grad = grad_field.levels() - ? array::make_view(grad_field).slice(Range::all(), Range::all(), Range::all()) - : array::make_view(grad_field).slice(Range::all(), Range::dummy(), Range::all()); + const mesh::Edges& edges = fvm_->mesh().edges(); + const mesh::Nodes& nodes = fvm_->mesh().nodes(); - const idx_t nlev = scalar.shape(1); - if (grad.shape(1) != nlev) { - throw_AssertionFailed("gradient field should have same number of levels", Here()); - } + const idx_t nnodes = fvm_->node_columns().nb_nodes(); + const idx_t nedges = fvm_->edge_columns().nb_edges(); - const auto lonlat_deg = array::make_view(nodes.lonlat()); - const auto dual_volumes = array::make_view(nodes.field("dual_volumes")); - const auto dual_normals = array::make_view(edges.field("dual_normals")); - const auto node2edge_sign = array::make_view(nodes.field("node2edge_sign")); + const auto scalar = scalar_field.levels() + ? array::make_view(scalar_field).slice(Range::all(), Range::all()) + : array::make_view(scalar_field).slice(Range::all(), Range::dummy()); + auto grad = grad_field.levels() + ? array::make_view(grad_field).slice(Range::all(), Range::all(), Range::all()) + : array::make_view(grad_field).slice(Range::all(), Range::dummy(), Range::all()); - const mesh::Connectivity& node2edge = nodes.edge_connectivity(); - const mesh::MultiBlockConnectivity& edge2node = edges.node_connectivity(); + const idx_t nlev = scalar.shape(1); + if (grad.shape(1) != nlev) { + throw_AssertionFailed("gradient field should have same number of levels", Here()); + } - array::ArrayT avgS_arr(nedges, nlev, 2ul); - auto avgS = array::make_view(avgS_arr); + const auto lonlat_deg = array::make_view(nodes.lonlat()); + const auto dual_volumes = array::make_view(nodes.field("dual_volumes")); + const auto dual_normals = array::make_view(edges.field("dual_normals")); + const auto node2edge_sign = array::make_view(nodes.field("node2edge_sign")); - const double scale = deg2rad * deg2rad * radius; + const mesh::Connectivity& node2edge = nodes.edge_connectivity(); + const mesh::MultiBlockConnectivity& edge2node = edges.node_connectivity(); - atlas_omp_parallel { - atlas_omp_for(idx_t jedge = 0; jedge < nedges; ++jedge) { - idx_t ip1 = edge2node(jedge, 0); - idx_t ip2 = edge2node(jedge, 1); + array::ArrayT avgS_arr(nedges, nlev, 2ul); + auto avgS = array::make_view(avgS_arr); - for (idx_t jlev = 0; jlev < nlev; ++jlev) { - double avg = (scalar(ip1, jlev) + scalar(ip2, jlev)) * 0.5; - avgS(jedge, jlev, LON) = dual_normals(jedge, LON) * deg2rad * avg; - avgS(jedge, jlev, LAT) = dual_normals(jedge, LAT) * deg2rad * avg; - } - } + const Value scale = deg2rad * deg2rad * radius; - atlas_omp_for(idx_t jnode = 0; jnode < nnodes; ++jnode) { - for (idx_t jlev = 0; jlev < nlev; ++jlev) { - grad(jnode, jlev, LON) = 0.; - grad(jnode, jlev, LAT) = 0.; + atlas_omp_parallel { + atlas_omp_for(idx_t jedge = 0; jedge < nedges; ++jedge) { + idx_t ip1 = edge2node(jedge, 0); + idx_t ip2 = edge2node(jedge, 1); + Value S[2] = { static_cast(dual_normals(jedge, LON)), static_cast(dual_normals(jedge, LAT))}; + + for (idx_t jlev = 0; jlev < nlev; ++jlev) { + Value avg = (scalar(ip1, jlev) + scalar(ip2, jlev)) * Value{0.5}; + avgS(jedge, jlev, LON) = S[LON] * deg2rad * avg; + avgS(jedge, jlev, LAT) = S[LAT] * deg2rad * avg; + } } - for (idx_t jedge = 0; jedge < node2edge.cols(jnode); ++jedge) { - const idx_t iedge = node2edge(jnode, jedge); - if (iedge < nedges) { - const double add = node2edge_sign(jnode, jedge); - for (idx_t jlev = 0; jlev < nlev; ++jlev) { - grad(jnode, jlev, LON) += add * avgS(iedge, jlev, LON); - grad(jnode, jlev, LAT) += add * avgS(iedge, jlev, LAT); + + atlas_omp_for(idx_t jnode = 0; jnode < nnodes; ++jnode) { + for (idx_t jlev = 0; jlev < nlev; ++jlev) { + grad(jnode, jlev, LON) = 0.; + grad(jnode, jlev, LAT) = 0.; + } + for (idx_t jedge = 0; jedge < node2edge.cols(jnode); ++jedge) { + const idx_t iedge = node2edge(jnode, jedge); + if (iedge < nedges) { + const Value add = node2edge_sign(jnode, jedge); + for (idx_t jlev = 0; jlev < nlev; ++jlev) { + grad(jnode, jlev, LON) += add * avgS(iedge, jlev, LON); + grad(jnode, jlev, LAT) += add * avgS(iedge, jlev, LAT); + } } } - } - const double y = lonlat_deg(jnode, LAT) * deg2rad; - const double metric_y = 1. / (dual_volumes(jnode) * scale); - const double metric_x = metric_y / std::cos(y); - for (idx_t jlev = 0; jlev < nlev; ++jlev) { - grad(jnode, jlev, LON) *= metric_x; - grad(jnode, jlev, LAT) *= metric_y; + const Value y = lonlat_deg(jnode, LAT) * deg2rad; + const Value metric_y = Value{1.} / (static_cast(dual_volumes(jnode)) * scale); + const Value metric_x = metric_y / std::cos(y); + for (idx_t jlev = 0; jlev < nlev; ++jlev) { + grad(jnode, jlev, LON) *= metric_x; + grad(jnode, jlev, LAT) *= metric_y; + } } } + }; + ATLAS_ASSERT( scalar_field.datatype() == grad_field.datatype() ); + switch (scalar_field.datatype().kind()) { + case (DataType::KIND_REAL32): { + dispatch(float{}); + break; + } + case (DataType::KIND_REAL64): { + dispatch(double{}); + break; + } + default: + ATLAS_NOTIMPLEMENTED; } } @@ -172,317 +191,375 @@ void Nabla::gradient_of_scalar(const Field& scalar_field, Field& grad_field) con void Nabla::gradient_of_vector(const Field& vector_field, Field& grad_field) const { Log::debug() << "Compute gradient of vector field " << vector_field.name() << " with fvm method" << std::endl; - const double radius = fvm_->radius(); - const double deg2rad = M_PI / 180.; - - const mesh::Edges& edges = fvm_->mesh().edges(); - const mesh::Nodes& nodes = fvm_->mesh().nodes(); - - const idx_t nnodes = fvm_->node_columns().nb_nodes(); - const idx_t nedges = fvm_->edge_columns().nb_edges(); - - const auto vector = - vector_field.levels() - ? array::make_view(vector_field).slice(Range::all(), Range::all(), Range::all()) - : array::make_view(vector_field).slice(Range::all(), Range::dummy(), Range::all()); - auto grad = grad_field.levels() - ? array::make_view(grad_field).slice(Range::all(), Range::all(), Range::all()) - : array::make_view(grad_field).slice(Range::all(), Range::dummy(), Range::all()); - - const idx_t nlev = vector.shape(1); - if (grad.shape(1) != nlev) { - throw_AssertionFailed("gradient field should have same number of levels", Here()); - } - const auto lonlat_deg = array::make_view(nodes.lonlat()); - const auto dual_volumes = array::make_view(nodes.field("dual_volumes")); - const auto dual_normals = array::make_view(edges.field("dual_normals")); - const auto node2edge_sign = array::make_view(nodes.field("node2edge_sign")); - const auto edge_flags = array::make_view(edges.flags()); - auto is_pole_edge = [&](idx_t e) { return Topology::check(edge_flags(e), Topology::POLE); }; + auto dispatch = [&](auto value) { + using Value = std::decay_t; - const mesh::Connectivity& node2edge = nodes.edge_connectivity(); - const mesh::MultiBlockConnectivity& edge2node = edges.node_connectivity(); + const Value radius = fvm_->radius(); + const Value deg2rad = M_PI / 180.; - array::ArrayT avgS_arr(nedges, nlev, 4ul); - array::ArrayView avgS = array::make_view(avgS_arr); + const mesh::Edges& edges = fvm_->mesh().edges(); + const mesh::Nodes& nodes = fvm_->mesh().nodes(); - const double scale = deg2rad * deg2rad * radius; + const idx_t nnodes = fvm_->node_columns().nb_nodes(); + const idx_t nedges = fvm_->edge_columns().nb_edges(); - enum - { - LONdLON = 0, - LONdLAT = 1, - LATdLON = 2, - LATdLAT = 3 - }; - - atlas_omp_parallel { - atlas_omp_for(idx_t jedge = 0; jedge < nedges; ++jedge) { - idx_t ip1 = edge2node(jedge, 0); - idx_t ip2 = edge2node(jedge, 1); - double pbc = 1. - 2. * is_pole_edge(jedge); + const auto vector = + vector_field.levels() + ? array::make_view(vector_field).slice(Range::all(), Range::all(), Range::all()) + : array::make_view(vector_field).slice(Range::all(), Range::dummy(), Range::all()); + auto grad = grad_field.levels() + ? array::make_view(grad_field).slice(Range::all(), Range::all(), Range::all()) + : array::make_view(grad_field).slice(Range::all(), Range::dummy(), Range::all()); - for (idx_t jlev = 0; jlev < nlev; ++jlev) { - double avg[2] = {(vector(ip1, jlev, LON) + pbc * vector(ip2, jlev, LON)) * 0.5, - (vector(ip1, jlev, LAT) + pbc * vector(ip2, jlev, LAT)) * 0.5}; - avgS(jedge, jlev, LONdLON) = dual_normals(jedge, LON) * deg2rad * avg[LON]; - // above = 0 at pole because of dual_normals - avgS(jedge, jlev, LONdLAT) = dual_normals(jedge, LAT) * deg2rad * avg[LON]; - avgS(jedge, jlev, LATdLON) = dual_normals(jedge, LON) * deg2rad * avg[LAT]; - // above = 0 at pole because of dual_normals - avgS(jedge, jlev, LATdLAT) = dual_normals(jedge, LAT) * deg2rad * avg[LAT]; - } + const idx_t nlev = vector.shape(1); + if (grad.shape(1) != nlev) { + throw_AssertionFailed("gradient field should have same number of levels", Here()); } - atlas_omp_for(idx_t jnode = 0; jnode < nnodes; ++jnode) { - for (idx_t jlev = 0; jlev < nlev; ++jlev) { - grad(jnode, jlev, LONdLON) = 0.; - grad(jnode, jlev, LONdLAT) = 0.; - grad(jnode, jlev, LATdLON) = 0.; - grad(jnode, jlev, LATdLAT) = 0.; + const auto lonlat_deg = array::make_view(nodes.lonlat()); + const auto dual_volumes = array::make_view(nodes.field("dual_volumes")); + const auto dual_normals = array::make_view(edges.field("dual_normals")); + const auto node2edge_sign = array::make_view(nodes.field("node2edge_sign")); + const auto edge_flags = array::make_view(edges.flags()); + auto is_pole_edge = [&](idx_t e) { return Topology::check(edge_flags(e), Topology::POLE); }; + + const mesh::Connectivity& node2edge = nodes.edge_connectivity(); + const mesh::MultiBlockConnectivity& edge2node = edges.node_connectivity(); + + array::ArrayT avgS_arr(nedges, nlev, 4ul); + array::ArrayView avgS = array::make_view(avgS_arr); + + const Value scale = deg2rad * deg2rad * radius; + + enum + { + LONdLON = 0, + LONdLAT = 1, + LATdLON = 2, + LATdLAT = 3 + }; + + atlas_omp_parallel { + atlas_omp_for(idx_t jedge = 0; jedge < nedges; ++jedge) { + idx_t ip1 = edge2node(jedge, 0); + idx_t ip2 = edge2node(jedge, 1); + Value pbc = 1. - 2. * is_pole_edge(jedge); + + Value S[2] = {static_cast(dual_normals(jedge,LON)), static_cast(dual_normals(jedge,LAT))}; + Value avg[2]; + + for (idx_t jlev = 0; jlev < nlev; ++jlev) { + avg[LON] = (vector(ip1, jlev, LON) + pbc * vector(ip2, jlev, LON)) * Value{0.5}; + avg[LAT] = (vector(ip1, jlev, LAT) + pbc * vector(ip2, jlev, LAT)) * Value{0.5}; + avgS(jedge, jlev, LONdLON) = S[LON] * deg2rad * avg[LON]; + // above = 0 at pole because of dual_normals + avgS(jedge, jlev, LONdLAT) = S[LAT] * deg2rad * avg[LON]; + avgS(jedge, jlev, LATdLON) = S[LON] * deg2rad * avg[LAT]; + // above = 0 at pole because of dual_normals + avgS(jedge, jlev, LATdLAT) = S[LAT] * deg2rad * avg[LAT]; + } } - for (idx_t jedge = 0; jedge < node2edge.cols(jnode); ++jedge) { - const idx_t iedge = node2edge(jnode, jedge); - if (iedge < nedges) { - double add = node2edge_sign(jnode, jedge); - for (idx_t jlev = 0; jlev < nlev; ++jlev) { - grad(jnode, jlev, LONdLON) += add * avgS(iedge, jlev, LONdLON); - grad(jnode, jlev, LONdLAT) += add * avgS(iedge, jlev, LONdLAT); - grad(jnode, jlev, LATdLON) += add * avgS(iedge, jlev, LATdLON); - grad(jnode, jlev, LATdLAT) += add * avgS(iedge, jlev, LATdLAT); + + atlas_omp_for(idx_t jnode = 0; jnode < nnodes; ++jnode) { + for (idx_t jlev = 0; jlev < nlev; ++jlev) { + grad(jnode, jlev, LONdLON) = 0.; + grad(jnode, jlev, LONdLAT) = 0.; + grad(jnode, jlev, LATdLON) = 0.; + grad(jnode, jlev, LATdLAT) = 0.; + } + for (idx_t jedge = 0; jedge < node2edge.cols(jnode); ++jedge) { + const idx_t iedge = node2edge(jnode, jedge); + if (iedge < nedges) { + Value add = node2edge_sign(jnode, jedge); + for (idx_t jlev = 0; jlev < nlev; ++jlev) { + grad(jnode, jlev, LONdLON) += add * avgS(iedge, jlev, LONdLON); + grad(jnode, jlev, LONdLAT) += add * avgS(iedge, jlev, LONdLAT); + grad(jnode, jlev, LATdLON) += add * avgS(iedge, jlev, LATdLON); + grad(jnode, jlev, LATdLAT) += add * avgS(iedge, jlev, LATdLAT); + } } } + const Value y = lonlat_deg(jnode, LAT) * deg2rad; + const Value metric_y = Value{1.} / (static_cast(dual_volumes(jnode)) * scale); + const Value metric_x = metric_y / std::cos(y); + for (idx_t jlev = 0; jlev < nlev; ++jlev) { + grad(jnode, jlev, LONdLON) *= metric_x; + grad(jnode, jlev, LATdLON) *= metric_x; + grad(jnode, jlev, LONdLAT) *= metric_y; + grad(jnode, jlev, LATdLAT) *= metric_y; + } } - const double y = lonlat_deg(jnode, LAT) * deg2rad; - const double metric_y = 1. / (dual_volumes(jnode) * scale); - const double metric_x = metric_y / std::cos(y); + } + // Fix wrong node2edge_sign for vector quantities + for (size_t jedge = 0; jedge < pole_edges_.size(); ++jedge) { + const idx_t iedge = pole_edges_[jedge]; + const idx_t jnode = edge2node(iedge, 1); + const Value metric_y = Value{1.} / (static_cast(dual_volumes(jnode)) * scale); for (idx_t jlev = 0; jlev < nlev; ++jlev) { - grad(jnode, jlev, LONdLON) *= metric_x; - grad(jnode, jlev, LATdLON) *= metric_x; - grad(jnode, jlev, LONdLAT) *= metric_y; - grad(jnode, jlev, LATdLAT) *= metric_y; + grad(jnode, jlev, LONdLAT) -= 2. * avgS(iedge, jlev, LONdLAT) * metric_y; + grad(jnode, jlev, LATdLAT) -= 2. * avgS(iedge, jlev, LATdLAT) * metric_y; } } - } - // Fix wrong node2edge_sign for vector quantities - for (size_t jedge = 0; jedge < pole_edges_.size(); ++jedge) { - const idx_t iedge = pole_edges_[jedge]; - const idx_t jnode = edge2node(iedge, 1); - const double metric_y = 1. / (dual_volumes(jnode) * scale); - for (idx_t jlev = 0; jlev < nlev; ++jlev) { - grad(jnode, jlev, LONdLAT) -= 2. * avgS(iedge, jlev, LONdLAT) * metric_y; - grad(jnode, jlev, LATdLAT) -= 2. * avgS(iedge, jlev, LATdLAT) * metric_y; + }; + ATLAS_ASSERT( vector_field.datatype() == grad_field.datatype() ); + switch (vector_field.datatype().kind()) { + case (DataType::KIND_REAL32): { + dispatch(float{}); + break; + } + case (DataType::KIND_REAL64): { + dispatch(double{}); + break; } + default: + ATLAS_NOTIMPLEMENTED; } } // ================================================================================ void Nabla::divergence(const Field& vector_field, Field& div_field) const { - const double radius = fvm_->radius(); - const double deg2rad = M_PI / 180.; - const mesh::Edges& edges = fvm_->mesh().edges(); - const mesh::Nodes& nodes = fvm_->mesh().nodes(); - - const idx_t nnodes = fvm_->node_columns().nb_nodes(); - const idx_t nedges = fvm_->edge_columns().nb_edges(); + auto dispatch = [&](auto value) { + using Value = std::decay_t; - const auto vector = - vector_field.levels() - ? array::make_view(vector_field).slice(Range::all(), Range::all(), Range::all()) - : array::make_view(vector_field).slice(Range::all(), Range::dummy(), Range::all()); - auto div = div_field.levels() ? array::make_view(div_field).slice(Range::all(), Range::all()) - : array::make_view(div_field).slice(Range::all(), Range::dummy()); - - const idx_t nlev = vector.shape(1); - if (div.shape(1) != nlev) { - throw_AssertionFailed("div_field should have same number of levels", Here()); - } + const Value radius = fvm_->radius(); + const Value deg2rad = M_PI / 180.; - const auto lonlat_deg = array::make_view(nodes.lonlat()); - const auto dual_volumes = array::make_view(nodes.field("dual_volumes")); - const auto dual_normals = array::make_view(edges.field("dual_normals")); - const auto node2edge_sign = array::make_view(nodes.field("node2edge_sign")); - const auto edge_flags = array::make_view(edges.flags()); - auto is_pole_edge = [&](idx_t e) { return Topology::check(edge_flags(e), Topology::POLE); }; + const mesh::Edges& edges = fvm_->mesh().edges(); + const mesh::Nodes& nodes = fvm_->mesh().nodes(); - const mesh::Connectivity& node2edge = nodes.edge_connectivity(); - const mesh::MultiBlockConnectivity& edge2node = edges.node_connectivity(); + const idx_t nnodes = fvm_->node_columns().nb_nodes(); + const idx_t nedges = fvm_->edge_columns().nb_edges(); - array::ArrayT avgS_arr(nedges, nlev, 2ul); - array::ArrayView avgS = array::make_view(avgS_arr); + const auto vector = + vector_field.levels() + ? array::make_view(vector_field).slice(Range::all(), Range::all(), Range::all()) + : array::make_view(vector_field).slice(Range::all(), Range::dummy(), Range::all()); + auto div = div_field.levels() ? array::make_view(div_field).slice(Range::all(), Range::all()) + : array::make_view(div_field).slice(Range::all(), Range::dummy()); - const double scale = deg2rad * deg2rad * radius; + const idx_t nlev = vector.shape(1); + if (div.shape(1) != nlev) { + throw_AssertionFailed("div_field should have same number of levels", Here()); + } - enum - { - LONdLON = 0, - LATdLAT = 1 - }; - atlas_omp_parallel { - atlas_omp_for(idx_t jedge = 0; jedge < nedges; ++jedge) { - double pbc = 1 - is_pole_edge(jedge); - - idx_t ip1 = edge2node(jedge, 0); - idx_t ip2 = edge2node(jedge, 1); - double y1 = lonlat_deg(ip1, LAT) * deg2rad; - double y2 = lonlat_deg(ip2, LAT) * deg2rad; - - double cosy1, cosy2; - if (metric_approach_ == 0) { - cosy1 = std::cos(y1) * pbc; - cosy2 = std::cos(y2) * pbc; - } - else { - cosy1 = cosy2 = std::cos(0.5 * (y1 + y2)) * pbc; - } + const auto lonlat_deg = array::make_view(nodes.lonlat()); + const auto dual_volumes = array::make_view(nodes.field("dual_volumes")); + const auto dual_normals = array::make_view(edges.field("dual_normals")); + const auto node2edge_sign = array::make_view(nodes.field("node2edge_sign")); + const auto edge_flags = array::make_view(edges.flags()); + auto is_pole_edge = [&](idx_t e) { return Topology::check(edge_flags(e), Topology::POLE); }; + + const mesh::Connectivity& node2edge = nodes.edge_connectivity(); + const mesh::MultiBlockConnectivity& edge2node = edges.node_connectivity(); + + array::ArrayT avgS_arr(nedges, nlev, 2ul); + array::ArrayView avgS = array::make_view(avgS_arr); + + const Value scale = deg2rad * deg2rad * radius; + + enum + { + LONdLON = 0, + LATdLAT = 1 + }; + atlas_omp_parallel { + atlas_omp_for(idx_t jedge = 0; jedge < nedges; ++jedge) { + Value pbc = 1 - is_pole_edge(jedge); + + idx_t ip1 = edge2node(jedge, 0); + idx_t ip2 = edge2node(jedge, 1); + Value y1 = lonlat_deg(ip1, LAT) * deg2rad; + Value y2 = lonlat_deg(ip2, LAT) * deg2rad; + + Value cosy1, cosy2; + if (metric_approach_ == 0) { + cosy1 = std::cos(y1) * pbc; + cosy2 = std::cos(y2) * pbc; + } + else { + cosy1 = cosy2 = std::cos(0.5 * (y1 + y2)) * pbc; + } - double S[2] = {dual_normals(jedge, LON) * deg2rad, dual_normals(jedge, LAT) * deg2rad}; - double avg[2]; + Value S[2] = {static_cast(dual_normals(jedge, LON)) * deg2rad, static_cast(dual_normals(jedge, LAT)) * deg2rad}; + Value avg[2]; - for (idx_t jlev = 0; jlev < nlev; ++jlev) { - double u1 = vector(ip1, jlev, LON); - double u2 = vector(ip2, jlev, LON); - double v1 = vector(ip1, jlev, LAT) * cosy1; - double v2 = vector(ip2, jlev, LAT) * cosy2; + for (idx_t jlev = 0; jlev < nlev; ++jlev) { + Value u1 = vector(ip1, jlev, LON); + Value u2 = vector(ip2, jlev, LON); + Value v1 = vector(ip1, jlev, LAT) * cosy1; + Value v2 = vector(ip2, jlev, LAT) * cosy2; - avg[LON] = (u1 + u2) * 0.5; - avg[LAT] = (v1 + v2) * 0.5; + avg[LON] = (u1 + u2) * 0.5; + avg[LAT] = (v1 + v2) * 0.5; - avgS(jedge, jlev, LONdLON) = avg[LON] * S[LON]; - avgS(jedge, jlev, LATdLAT) = avg[LAT] * S[LAT]; + avgS(jedge, jlev, LONdLON) = avg[LON] * S[LON]; + avgS(jedge, jlev, LATdLAT) = avg[LAT] * S[LAT]; + } } - } - atlas_omp_for(idx_t jnode = 0; jnode < nnodes; ++jnode) { - for (idx_t jlev = 0; jlev < nlev; ++jlev) { - div(jnode, jlev) = 0.; - } - for (idx_t jedge = 0; jedge < node2edge.cols(jnode); ++jedge) { - idx_t iedge = node2edge(jnode, jedge); - if (iedge < nedges) { - double add = node2edge_sign(jnode, jedge); - for (idx_t jlev = 0; jlev < nlev; ++jlev) { - div(jnode, jlev) += add * (avgS(iedge, jlev, LONdLON) + avgS(iedge, jlev, LATdLAT)); + atlas_omp_for(idx_t jnode = 0; jnode < nnodes; ++jnode) { + for (idx_t jlev = 0; jlev < nlev; ++jlev) { + div(jnode, jlev) = 0.; + } + for (idx_t jedge = 0; jedge < node2edge.cols(jnode); ++jedge) { + idx_t iedge = node2edge(jnode, jedge); + if (iedge < nedges) { + Value add = node2edge_sign(jnode, jedge); + for (idx_t jlev = 0; jlev < nlev; ++jlev) { + div(jnode, jlev) += add * (avgS(iedge, jlev, LONdLON) + avgS(iedge, jlev, LATdLAT)); + } } } + const Value y = lonlat_deg(jnode, LAT) * deg2rad; + Value metric = Value{1.} / (static_cast(dual_volumes(jnode)) * scale * std::cos(y)); + for (idx_t jlev = 0; jlev < nlev; ++jlev) { + div(jnode, jlev) *= metric; + } } - const double y = lonlat_deg(jnode, LAT) * deg2rad; - double metric = 1. / (dual_volumes(jnode) * scale * std::cos(y)); - for (idx_t jlev = 0; jlev < nlev; ++jlev) { - div(jnode, jlev) *= metric; - } } + + }; + ATLAS_ASSERT( vector_field.datatype() == div_field.datatype() ); + switch (vector_field.datatype().kind()) { + case (DataType::KIND_REAL32): { + dispatch(float{}); + break; + } + case (DataType::KIND_REAL64): { + dispatch(double{}); + break; + } + default: + ATLAS_NOTIMPLEMENTED; } } void Nabla::curl(const Field& vector_field, Field& curl_field) const { - const double radius = fvm_->radius(); - const double deg2rad = M_PI / 180.; + auto dispatch = [&](auto value) { + using Value = std::decay_t; - const mesh::Edges& edges = fvm_->mesh().edges(); - const mesh::Nodes& nodes = fvm_->mesh().nodes(); - const idx_t nnodes = fvm_->node_columns().nb_nodes(); - const idx_t nedges = fvm_->edge_columns().nb_edges(); + const Value radius = fvm_->radius(); + const Value deg2rad = M_PI / 180.; - const auto vector = - vector_field.levels() - ? array::make_view(vector_field).slice(Range::all(), Range::all(), Range::all()) - : array::make_view(vector_field).slice(Range::all(), Range::dummy(), Range::all()); - auto curl = curl_field.levels() ? array::make_view(curl_field).slice(Range::all(), Range::all()) - : array::make_view(curl_field).slice(Range::all(), Range::dummy()); + const mesh::Edges& edges = fvm_->mesh().edges(); + const mesh::Nodes& nodes = fvm_->mesh().nodes(); - const idx_t nlev = vector.shape(1); - if (curl.shape(1) != nlev) { - throw_AssertionFailed("curl field should have same number of levels", Here()); - } + const idx_t nnodes = fvm_->node_columns().nb_nodes(); + const idx_t nedges = fvm_->edge_columns().nb_edges(); - const auto lonlat_deg = array::make_view(nodes.lonlat()); - const auto dual_volumes = array::make_view(nodes.field("dual_volumes")); - const auto dual_normals = array::make_view(edges.field("dual_normals")); - const auto node2edge_sign = array::make_view(nodes.field("node2edge_sign")); - const auto edge_flags = array::make_view(edges.flags()); - auto is_pole_edge = [&](idx_t e) { return Topology::check(edge_flags(e), Topology::POLE); }; + const auto vector = + vector_field.levels() + ? array::make_view(vector_field).slice(Range::all(), Range::all(), Range::all()) + : array::make_view(vector_field).slice(Range::all(), Range::dummy(), Range::all()); + auto curl = curl_field.levels() ? array::make_view(curl_field).slice(Range::all(), Range::all()) + : array::make_view(curl_field).slice(Range::all(), Range::dummy()); + const idx_t nlev = vector.shape(1); + if (curl.shape(1) != nlev) { + throw_AssertionFailed("curl field should have same number of levels", Here()); + } - const mesh::Connectivity& node2edge = nodes.edge_connectivity(); - const mesh::MultiBlockConnectivity& edge2node = edges.node_connectivity(); + const auto lonlat_deg = array::make_view(nodes.lonlat()); + const auto dual_volumes = array::make_view(nodes.field("dual_volumes")); + const auto dual_normals = array::make_view(edges.field("dual_normals")); + const auto node2edge_sign = array::make_view(nodes.field("node2edge_sign")); + const auto edge_flags = array::make_view(edges.flags()); + auto is_pole_edge = [&](idx_t e) { return Topology::check(edge_flags(e), Topology::POLE); }; - array::ArrayT avgS_arr(nedges, nlev, 2ul); - array::ArrayView avgS = array::make_view(avgS_arr); - const double scale = deg2rad * deg2rad * radius; + const mesh::Connectivity& node2edge = nodes.edge_connectivity(); + const mesh::MultiBlockConnectivity& edge2node = edges.node_connectivity(); - enum - { - LONdLAT = 0, - LATdLON = 1 - }; + array::ArrayT avgS_arr(nedges, nlev, 2ul); + array::ArrayView avgS = array::make_view(avgS_arr); + const Value scale = deg2rad * deg2rad * radius; - atlas_omp_parallel { - atlas_omp_for(idx_t jedge = 0; jedge < nedges; ++jedge) { - idx_t ip1 = edge2node(jedge, 0); - idx_t ip2 = edge2node(jedge, 1); - double y1 = lonlat_deg(ip1, LAT) * deg2rad; - double y2 = lonlat_deg(ip2, LAT) * deg2rad; + enum + { + LONdLAT = 0, + LATdLON = 1 + }; - double pbc = 1 - is_pole_edge(jedge); - double cosy1; - double cosy2; - if (metric_approach_ == 0) { - cosy1 = std::cos(y1) * pbc; - cosy2 = std::cos(y2) * pbc; - } - else { - cosy1 = cosy2 = std::cos(0.5 * (y1 + y2)) * pbc; - } - double S[2] = {dual_normals(jedge, LON) * deg2rad, dual_normals(jedge, LAT) * deg2rad}; - double avg[2]; + atlas_omp_parallel { + atlas_omp_for(idx_t jedge = 0; jedge < nedges; ++jedge) { + idx_t ip1 = edge2node(jedge, 0); + idx_t ip2 = edge2node(jedge, 1); + Value y1 = lonlat_deg(ip1, LAT) * deg2rad; + Value y2 = lonlat_deg(ip2, LAT) * deg2rad; - for (idx_t jlev = 0; jlev < nlev; ++jlev) { - double u1 = vector(ip1, jlev, LON) * cosy1; - double u2 = vector(ip2, jlev, LON) * cosy2; - double v1 = vector(ip1, jlev, LAT); - double v2 = vector(ip2, jlev, LAT); + Value pbc = 1 - is_pole_edge(jedge); + Value cosy1; + Value cosy2; + if (metric_approach_ == 0) { + cosy1 = std::cos(y1) * pbc; + cosy2 = std::cos(y2) * pbc; + } + else { + cosy1 = cosy2 = std::cos(Value{0.5} * (y1 + y2)) * pbc; + } - avg[LON] = (u1 + u2) * 0.5; - avg[LAT] = (v1 + v2) * 0.5; + Value S[2] = {static_cast(dual_normals(jedge, LON)) * deg2rad, static_cast(dual_normals(jedge, LAT)) * deg2rad}; + Value avg[2]; - avgS(jedge, jlev, LONdLAT) = avg[LON] * S[LAT]; - avgS(jedge, jlev, LATdLON) = avg[LAT] * S[LON]; - } - } + for (idx_t jlev = 0; jlev < nlev; ++jlev) { + Value u1 = vector(ip1, jlev, LON) * cosy1; + Value u2 = vector(ip2, jlev, LON) * cosy2; + Value v1 = vector(ip1, jlev, LAT); + Value v2 = vector(ip2, jlev, LAT); - atlas_omp_for(idx_t jnode = 0; jnode < nnodes; ++jnode) { - for (idx_t jlev = 0; jlev < nlev; ++jlev) { - curl(jnode, jlev) = 0.; + avg[LON] = (u1 + u2) * Value{0.5}; + avg[LAT] = (v1 + v2) * Value{0.5}; + + avgS(jedge, jlev, LONdLAT) = avg[LON] * S[LAT]; + avgS(jedge, jlev, LATdLON) = avg[LAT] * S[LON]; + } } - for (idx_t jedge = 0; jedge < node2edge.cols(jnode); ++jedge) { - idx_t iedge = node2edge(jnode, jedge); - if (iedge < nedges) { - double add = node2edge_sign(jnode, jedge); - for (idx_t jlev = 0; jlev < nlev; ++jlev) { - curl(jnode, jlev) += add * (avgS(iedge, jlev, LATdLON) - avgS(iedge, jlev, LONdLAT)); + + atlas_omp_for(idx_t jnode = 0; jnode < nnodes; ++jnode) { + for (idx_t jlev = 0; jlev < nlev; ++jlev) { + curl(jnode, jlev) = 0.; + } + for (idx_t jedge = 0; jedge < node2edge.cols(jnode); ++jedge) { + idx_t iedge = node2edge(jnode, jedge); + if (iedge < nedges) { + double add = node2edge_sign(jnode, jedge); + for (idx_t jlev = 0; jlev < nlev; ++jlev) { + curl(jnode, jlev) += add * (avgS(iedge, jlev, LATdLON) - avgS(iedge, jlev, LONdLAT)); + } } } + Value y = static_cast(lonlat_deg(jnode, LAT)) * deg2rad; + Value metric = Value{1.} / (static_cast(dual_volumes(jnode)) * scale * std::cos(y)); + for (idx_t jlev = 0; jlev < nlev; ++jlev) { + curl(jnode, jlev) *= metric; + } } - double y = lonlat_deg(jnode, LAT) * deg2rad; - double metric = 1. / (dual_volumes(jnode) * scale * std::cos(y)); - for (idx_t jlev = 0; jlev < nlev; ++jlev) { - curl(jnode, jlev) *= metric; - } } + }; + ATLAS_ASSERT( vector_field.datatype() == curl_field.datatype() ); + switch (vector_field.datatype().kind()) { + case (DataType::KIND_REAL32): { + dispatch(float{}); + break; + } + case (DataType::KIND_REAL64): { + dispatch(double{}); + break; + } + default: + ATLAS_NOTIMPLEMENTED; } } void Nabla::laplacian(const Field& scalar, Field& lapl) const { - Field grad(fvm_->node_columns().createField(option::name("grad") | option::levels(scalar.levels()) | - option::variables(2))); + Field grad(fvm_->node_columns().createField(option::name("grad") | option::levels(scalar.levels()) | + option::variables(2) | option::datatype(scalar.datatype()))); gradient(scalar, grad); if (fvm_->node_columns().halo().size() < 2) { fvm_->node_columns().haloExchange(grad); diff --git a/src/tests/numerics/test_fvm_nabla.cc b/src/tests/numerics/test_fvm_nabla.cc index 167cb698a..82ae1d2b7 100644 --- a/src/tests/numerics/test_fvm_nabla.cc +++ b/src/tests/numerics/test_fvm_nabla.cc @@ -57,16 +57,18 @@ static std::string griduid() { return "Slat20"; } -array::LocalView make_vectorview(Field& field) { +template +array::LocalView make_vectorview(Field& field) { using array::Range; - return field.levels() ? array::make_view(field).slice(Range::all(), Range::all(), Range::all()) - : array::make_view(field).slice(Range::all(), Range::dummy(), Range::all()); + return field.levels() ? array::make_view(field).slice(Range::all(), Range::all(), Range::all()) + : array::make_view(field).slice(Range::all(), Range::dummy(), Range::all()); } -array::LocalView make_scalarview(Field& field) { +template +array::LocalView make_scalarview(Field& field) { using array::Range; - return field.levels() ? array::make_view(field).slice(Range::all(), Range::all()) - : array::make_view(field).slice(Range::all(), Range::dummy()); + return field.levels() ? array::make_view(field).slice(Range::all(), Range::all()) + : array::make_view(field).slice(Range::all(), Range::dummy()); } @@ -100,6 +102,7 @@ double dual_volume(const Mesh& mesh) { /// @brief Compute magnitude of flow with rotation-angle beta /// (beta=0 --> zonal, beta=pi/2 --> meridional) +template void rotated_flow(const fvm::Method& fvm, Field& field, const double& beta) { const double radius = fvm.radius(); const double USCAL = 20.; @@ -107,7 +110,7 @@ void rotated_flow(const fvm::Method& fvm, Field& field, const double& beta) { const double deg2rad = M_PI / 180.; auto lonlat_deg = array::make_view(fvm.mesh().nodes().lonlat()); - auto var = make_vectorview(field); + auto var = make_vectorview(field); const idx_t nnodes = fvm.mesh().nodes().size(); const idx_t nlev = var.shape(1); @@ -125,6 +128,7 @@ void rotated_flow(const fvm::Method& fvm, Field& field, const double& beta) { /// @brief Compute magnitude of flow with rotation-angle beta /// (beta=0 --> zonal, beta=pi/2 --> meridional) +template void rotated_flow_magnitude(const fvm::Method& fvm, Field& field, const double& beta) { const double radius = fvm.radius(); const double USCAL = 20.; @@ -132,7 +136,7 @@ void rotated_flow_magnitude(const fvm::Method& fvm, Field& field, const double& const double deg2rad = M_PI / 180.; auto lonlat_deg = array::make_view(fvm.mesh().nodes().lonlat()); - auto var = make_scalarview(field); + auto var = make_scalarview(field); const idx_t nnodes = fvm.mesh().nodes().size(); const idx_t nlev = var.shape(1); @@ -166,7 +170,6 @@ CASE("test_build") { } CASE("test_grad") { - Log::info() << "test_grad" << std::endl; auto radius = option::radius("Earth"); Grid grid(griduid()); MeshGenerator meshgenerator("structured"); @@ -177,93 +180,102 @@ CASE("test_grad") { idx_t nnodes = mesh.nodes().size(); idx_t nlev = std::max(1, test_levels()); - FieldSet fields; - fields.add(fvm.node_columns().createField(option::name("scalar"))); - fields.add(fvm.node_columns().createField(option::name("rscalar"))); - fields.add(fvm.node_columns().createField(option::name("grad") | option::variables(2))); - fields.add(fvm.node_columns().createField(option::name("rgrad") | option::variables(2))); - fields.add(fvm.node_columns().createField(option::name("xder"))); - fields.add(fvm.node_columns().createField(option::name("yder"))); - fields.add(fvm.node_columns().createField(option::name("rxder"))); - fields.add(fvm.node_columns().createField(option::name("ryder"))); - - EXPECT(fields["scalar"].rank() == (test_levels() ? 2 : 1)); - EXPECT(fields["grad"].rank() == fields["scalar"].rank() + 1); - EXPECT(fields["scalar"].levels() == test_levels()); - EXPECT(fields["grad"].levels() == test_levels()); - - rotated_flow_magnitude(fvm, fields["scalar"], 0.); - rotated_flow_magnitude(fvm, fields["rscalar"], M_PI_2 * 0.75); - - nabla.gradient(fields["scalar"], fields["grad"]); - nabla.gradient(fields["rscalar"], fields["rgrad"]); - auto xder = make_scalarview(fields["xder"]); - auto yder = make_scalarview(fields["yder"]); - auto rxder = make_scalarview(fields["rxder"]); - auto ryder = make_scalarview(fields["ryder"]); - const auto grad = make_vectorview(fields["grad"]); - const auto rgrad = make_vectorview(fields["rgrad"]); - for (idx_t jnode = 0; jnode < nnodes; ++jnode) { - for (idx_t jlev = 0; jlev < nlev; ++jlev) { - xder(jnode, jlev) = grad(jnode, jlev, LON); - yder(jnode, jlev) = grad(jnode, jlev, LAT); - rxder(jnode, jlev) = rgrad(jnode, jlev, LON); - ryder(jnode, jlev) = rgrad(jnode, jlev, LAT); + auto do_test = [&](auto value, double tolerance) { + using Value = std::decay_t; + FieldSet fields; + fields.add(fvm.node_columns().createField(option::name("scalar"))); + fields.add(fvm.node_columns().createField(option::name("rscalar"))); + fields.add(fvm.node_columns().createField(option::name("grad") | option::variables(2))); + fields.add(fvm.node_columns().createField(option::name("rgrad") | option::variables(2))); + fields.add(fvm.node_columns().createField(option::name("xder"))); + fields.add(fvm.node_columns().createField(option::name("yder"))); + fields.add(fvm.node_columns().createField(option::name("rxder"))); + fields.add(fvm.node_columns().createField(option::name("ryder"))); + + EXPECT(fields["scalar"].rank() == (test_levels() ? 2 : 1)); + EXPECT(fields["grad"].rank() == fields["scalar"].rank() + 1); + EXPECT(fields["scalar"].levels() == test_levels()); + EXPECT(fields["grad"].levels() == test_levels()); + + rotated_flow_magnitude(fvm, fields["scalar"], 0.); + rotated_flow_magnitude(fvm, fields["rscalar"], M_PI_2 * 0.75); + + nabla.gradient(fields["scalar"], fields["grad"]); + nabla.gradient(fields["rscalar"], fields["rgrad"]); + auto xder = make_scalarview(fields["xder"]); + auto yder = make_scalarview(fields["yder"]); + auto rxder = make_scalarview(fields["rxder"]); + auto ryder = make_scalarview(fields["ryder"]); + const auto grad = make_vectorview(fields["grad"]); + const auto rgrad = make_vectorview(fields["rgrad"]); + for (idx_t jnode = 0; jnode < nnodes; ++jnode) { + for (idx_t jlev = 0; jlev < nlev; ++jlev) { + xder(jnode, jlev) = grad(jnode, jlev, LON); + yder(jnode, jlev) = grad(jnode, jlev, LAT); + rxder(jnode, jlev) = rgrad(jnode, jlev, LON); + ryder(jnode, jlev) = rgrad(jnode, jlev, LAT); + } } - } - // output to gmsh - { - fvm.node_columns().haloExchange(fields); - output::Gmsh(grid.name() + ".msh").write(mesh); - output::Gmsh gmsh_fields(grid.name() + "_fields.msh"); - gmsh_fields.write(fields["scalar"]); - gmsh_fields.write(fields["xder"]); - gmsh_fields.write(fields["yder"]); - gmsh_fields.write(fields["rscalar"]); - gmsh_fields.write(fields["rxder"]); - gmsh_fields.write(fields["ryder"]); - } + // output to gmsh + { + fvm.node_columns().haloExchange(fields); + output::Gmsh(grid.name() + ".msh").write(mesh); + output::Gmsh gmsh_fields(grid.name() + "_fields.msh"); + gmsh_fields.write(fields["scalar"]); + gmsh_fields.write(fields["xder"]); + gmsh_fields.write(fields["yder"]); + gmsh_fields.write(fields["rscalar"]); + gmsh_fields.write(fields["rxder"]); + gmsh_fields.write(fields["ryder"]); + } - double min, max, mean; - idx_t N; - - fvm.node_columns().minimum(fields["xder"], min); - fvm.node_columns().maximum(fields["xder"], max); - fvm.node_columns().mean(fields["xder"], mean, N); - print_min_max_mean("xder"); - EXPECT_APPROX_EQ(min, 0., 1.e-20); - EXPECT_APPROX_EQ(max, 0., 1.e-20); - EXPECT_APPROX_EQ(mean, 0., 1.e-20); - - fvm.node_columns().minimum(fields["yder"], min); - fvm.node_columns().maximum(fields["yder"], max); - fvm.node_columns().mean(fields["yder"], mean, N); - print_min_max_mean("yder"); - EXPECT_APPROX_EQ(min, -3.1141489788326316614e-06); - EXPECT_APPROX_EQ(max, 3.1141489788326316614e-06); - EXPECT_APPROX_EQ(mean, 0., 1.e-20); - - - fvm.node_columns().minimum(fields["rxder"], min); - fvm.node_columns().maximum(fields["rxder"], max); - fvm.node_columns().mean(fields["rxder"], mean, N); - print_min_max_mean("rxder"); - EXPECT_APPROX_EQ(min, -3.02863817262107e-06); - EXPECT_APPROX_EQ(max, +3.02863817262107e-06); - EXPECT_APPROX_EQ(mean, 0., 1.e-20); - - fvm.node_columns().minimum(fields["ryder"], min); - fvm.node_columns().maximum(fields["ryder"], max); - fvm.node_columns().mean(fields["ryder"], mean, N); - print_min_max_mean("ryder"); - EXPECT_APPROX_EQ(min, -3.114148978832633e-06); - EXPECT_APPROX_EQ(max, +3.114148978832633e-06); - EXPECT_APPROX_EQ(mean, 0., 1.e-20); + double min, max, mean; + idx_t N; + + fvm.node_columns().minimum(fields["xder"], min); + fvm.node_columns().maximum(fields["xder"], max); + fvm.node_columns().mean(fields["xder"], mean, N); + print_min_max_mean("xder"); + EXPECT_APPROX_EQ(min, 0., tolerance); + EXPECT_APPROX_EQ(max, 0., tolerance); + EXPECT_APPROX_EQ(mean, 0., tolerance); + + fvm.node_columns().minimum(fields["yder"], min); + fvm.node_columns().maximum(fields["yder"], max); + fvm.node_columns().mean(fields["yder"], mean, N); + print_min_max_mean("yder"); + EXPECT_APPROX_EQ(min, -3.1141489788326316614e-06, tolerance); + EXPECT_APPROX_EQ(max, 3.1141489788326316614e-06, tolerance); + EXPECT_APPROX_EQ(mean, 0., tolerance); + + + fvm.node_columns().minimum(fields["rxder"], min); + fvm.node_columns().maximum(fields["rxder"], max); + fvm.node_columns().mean(fields["rxder"], mean, N); + print_min_max_mean("rxder"); + EXPECT_APPROX_EQ(min, -3.02863817262107e-06, tolerance); + EXPECT_APPROX_EQ(max, +3.02863817262107e-06, tolerance); + EXPECT_APPROX_EQ(mean, 0., tolerance); + + fvm.node_columns().minimum(fields["ryder"], min); + fvm.node_columns().maximum(fields["ryder"], max); + fvm.node_columns().mean(fields["ryder"], mean, N); + print_min_max_mean("ryder"); + EXPECT_APPROX_EQ(min, -3.114148978832633e-06, tolerance); + EXPECT_APPROX_EQ(max, +3.114148978832633e-06, tolerance); + EXPECT_APPROX_EQ(mean, 0., tolerance); + }; + SECTION("double precision") { + do_test(double{}, 1.e-20); + } + SECTION("single precision") { + do_test(float{}, 1.e-10); + } } + CASE("test_div") { - Log::info() << "test_div" << std::endl; const double radius = util::Earth::radius(); // const double radius = 1.; Grid grid(griduid()); @@ -272,33 +284,44 @@ CASE("test_div") { fvm::Method fvm(mesh, util::Config("radius", radius) | option::levels(test_levels())); Nabla nabla(fvm); - FieldSet fields; - fields.add(fvm.node_columns().createField(option::name("wind") | option::variables(2))); - fields.add(fvm.node_columns().createField(option::name("div"))); + auto do_test = [&](auto value, double tolerance) { + using Value = std::decay_t; - rotated_flow(fvm, fields["wind"], M_PI_2 * 0.75); + FieldSet fields; + fields.add(fvm.node_columns().createField(option::name("wind") | option::variables(2))); + fields.add(fvm.node_columns().createField(option::name("div"))); - nabla.divergence(fields["wind"], fields["div"]); + rotated_flow(fvm, fields["wind"], M_PI_2 * 0.75); - // output to gmsh - { - fvm.node_columns().haloExchange(fields); - output::Gmsh gmsh(grid.name() + "_fields.msh", "a"); - gmsh.write(fields["wind"]); - gmsh.write(fields["div"]); - } + nabla.divergence(fields["wind"], fields["div"]); - double min, max, mean; - idx_t N; - fvm.node_columns().minimum(fields["div"], min); - fvm.node_columns().maximum(fields["div"], max); - fvm.node_columns().mean(fields["div"], mean, N); - - // Divergence free flow! - print_min_max_mean("div"); - EXPECT_APPROX_EQ(min, 0., 1.e-18); - EXPECT_APPROX_EQ(max, 0., 1.e-18); - EXPECT_APPROX_EQ(mean, 0., 1.e-20); + // output to gmsh + { + fvm.node_columns().haloExchange(fields); + output::Gmsh gmsh(grid.name() + "_fields.msh", "a"); + gmsh.write(fields["wind"]); + gmsh.write(fields["div"]); + } + + double min, max, mean; + idx_t N; + fvm.node_columns().minimum(fields["div"], min); + fvm.node_columns().maximum(fields["div"], max); + fvm.node_columns().mean(fields["div"], mean, N); + + // Divergence free flow! + print_min_max_mean("div"); + EXPECT_APPROX_EQ(min, 0.,tolerance); + EXPECT_APPROX_EQ(max, 0., tolerance); + EXPECT_APPROX_EQ(mean, 0., tolerance); + }; + SECTION("double precision") { + do_test(double{}, 1.e-18); + } + SECTION("single precision") { + do_test(float{}, 1.e-10); + } + } CASE("test_curl") { @@ -311,71 +334,81 @@ CASE("test_curl") { fvm::Method fvm(mesh, util::Config("radius", radius) | option::levels(test_levels())); Nabla nabla(fvm); - FieldSet fields; - fields.add(fvm.node_columns().createField(option::name("wind") | option::variables(2))); - fields.add(fvm.node_columns().createField(option::name("vor"))); - - rotated_flow(fvm, fields["wind"], M_PI_2 * 0.75); - - nabla.curl(fields["wind"], fields["vor"]); - - fields.add(fvm.node_columns().createField(option::name("windgrad") | option::variables(2 * 2))); - nabla.gradient(fields["wind"], fields["windgrad"]); - - fields.add(fvm.node_columns().createField(option::name("windX") | option::levels(false))); - fields.add(fvm.node_columns().createField(option::name("windY") | option::levels(false))); - fields.add(fvm.node_columns().createField(option::name("windXgradX"))); - fields.add(fvm.node_columns().createField(option::name("windXgradY"))); - fields.add(fvm.node_columns().createField(option::name("windYgradX"))); - fields.add(fvm.node_columns().createField(option::name("windYgradY"))); - auto wind = make_vectorview(fields["wind"]); - auto windgrad = make_vectorview(fields["windgrad"]); - - auto windX = array::make_view(fields["windX"]); - auto windY = array::make_view(fields["windY"]); - auto windXgradX = make_scalarview(fields["windXgradX"]); - auto windXgradY = make_scalarview(fields["windXgradY"]); - auto windYgradX = make_scalarview(fields["windYgradX"]); - auto windYgradY = make_scalarview(fields["windYgradY"]); - for (idx_t j = 0; j < windX.size(); ++j) { - static const idx_t lev0 = 0; - static const idx_t XdX = XX * 2 + XX; - static const idx_t XdY = XX * 2 + YY; - static const idx_t YdX = YY * 2 + XX; - static const idx_t YdY = YY * 2 + YY; - windX(j) = wind(j, lev0, XX); - windY(j) = wind(j, lev0, YY); - windXgradX(j, lev0) = windgrad(j, lev0, XdX); - windXgradY(j, lev0) = windgrad(j, lev0, XdY); - windYgradX(j, lev0) = windgrad(j, lev0, YdX); - windYgradY(j, lev0) = windgrad(j, lev0, YdY); - } + auto do_test = [&](auto value, double tolerance) { + using Value = std::decay_t; + + FieldSet fields; + fields.add(fvm.node_columns().createField(option::name("wind") | option::variables(2))); + fields.add(fvm.node_columns().createField(option::name("vor"))); + + rotated_flow(fvm, fields["wind"], M_PI_2 * 0.75); + + nabla.curl(fields["wind"], fields["vor"]); + + fields.add(fvm.node_columns().createField(option::name("windgrad") | option::variables(2 * 2))); + nabla.gradient(fields["wind"], fields["windgrad"]); + + fields.add(fvm.node_columns().createField(option::name("windX") | option::levels(false))); + fields.add(fvm.node_columns().createField(option::name("windY") | option::levels(false))); + fields.add(fvm.node_columns().createField(option::name("windXgradX"))); + fields.add(fvm.node_columns().createField(option::name("windXgradY"))); + fields.add(fvm.node_columns().createField(option::name("windYgradX"))); + fields.add(fvm.node_columns().createField(option::name("windYgradY"))); + auto wind = make_vectorview(fields["wind"]); + auto windgrad = make_vectorview(fields["windgrad"]); + + auto windX = array::make_view(fields["windX"]); + auto windY = array::make_view(fields["windY"]); + auto windXgradX = make_scalarview(fields["windXgradX"]); + auto windXgradY = make_scalarview(fields["windXgradY"]); + auto windYgradX = make_scalarview(fields["windYgradX"]); + auto windYgradY = make_scalarview(fields["windYgradY"]); + for (idx_t j = 0; j < windX.size(); ++j) { + static const idx_t lev0 = 0; + static const idx_t XdX = XX * 2 + XX; + static const idx_t XdY = XX * 2 + YY; + static const idx_t YdX = YY * 2 + XX; + static const idx_t YdY = YY * 2 + YY; + windX(j) = wind(j, lev0, XX); + windY(j) = wind(j, lev0, YY); + windXgradX(j, lev0) = windgrad(j, lev0, XdX); + windXgradY(j, lev0) = windgrad(j, lev0, XdY); + windYgradX(j, lev0) = windgrad(j, lev0, YdX); + windYgradY(j, lev0) = windgrad(j, lev0, YdY); + } - // output to gmsh - { - fvm.node_columns().haloExchange(fields); - output::Gmsh gmsh(grid.name() + "_fields.msh", "a"); - gmsh.write(fields["vor"]); - gmsh.write(fields["windX"]); - gmsh.write(fields["windXgradX"]); - gmsh.write(fields["windXgradY"]); - gmsh.write(fields["windY"]); - gmsh.write(fields["windYgradX"]); - gmsh.write(fields["windYgradY"]); - gmsh.write(fields["windgrad"]); - } + // output to gmsh + { + fvm.node_columns().haloExchange(fields); + output::Gmsh gmsh(grid.name() + "_fields.msh", "a"); + gmsh.write(fields["vor"]); + gmsh.write(fields["windX"]); + gmsh.write(fields["windXgradX"]); + gmsh.write(fields["windXgradY"]); + gmsh.write(fields["windY"]); + gmsh.write(fields["windYgradX"]); + gmsh.write(fields["windYgradY"]); + gmsh.write(fields["windgrad"]); + } - double min, max, mean; - idx_t N; - - // Vorticity! - fvm.node_columns().minimum(fields["vor"], min); - fvm.node_columns().maximum(fields["vor"], max); - fvm.node_columns().mean(fields["vor"], mean, N); - print_min_max_mean("vor"); - EXPECT_APPROX_EQ(min, -6.257451225821150e-06); - EXPECT_APPROX_EQ(max, 6.257451225821150e-06); - EXPECT_APPROX_EQ(mean, 0., 1.e-20); + double min, max, mean; + idx_t N; + + // Vorticity! + fvm.node_columns().minimum(fields["vor"], min); + fvm.node_columns().maximum(fields["vor"], max); + fvm.node_columns().mean(fields["vor"], mean, N); + print_min_max_mean("vor"); + EXPECT_APPROX_EQ(min, -6.257451225821150e-06, tolerance); + EXPECT_APPROX_EQ(max, 6.257451225821150e-06, tolerance); + EXPECT_APPROX_EQ(mean, 0., tolerance); + }; + SECTION("double precision") { + do_test(double{}, 1.e-18); + } + SECTION("single precision") { + do_test(float{}, 1.e-10); + } } CASE("test_lapl") { @@ -388,30 +421,42 @@ CASE("test_lapl") { fvm::Method fvm(mesh, util::Config("radius", radius) | option::levels(test_levels())); Nabla nabla(fvm); - FieldSet fields; - fields.add(fvm.node_columns().createField(option::name("scal"))); - fields.add(fvm.node_columns().createField(option::name("lapl"))); - rotated_flow_magnitude(fvm, fields["scal"], M_PI_2 * 0.75); + auto do_test = [&](auto value, double tolerance) { + using Value = std::decay_t; + + FieldSet fields; + fields.add(fvm.node_columns().createField(option::name("scal"))); + fields.add(fvm.node_columns().createField(option::name("lapl"))); + + rotated_flow_magnitude(fvm, fields["scal"], M_PI_2 * 0.75); - nabla.laplacian(fields["scal"], fields["lapl"]); + nabla.laplacian(fields["scal"], fields["lapl"]); - // output to gmsh - { - fvm.node_columns().haloExchange(fields); - output::Gmsh gmsh(grid.name() + "_fields.msh", "a"); - gmsh.write(fields["lapl"]); + // output to gmsh + { + fvm.node_columns().haloExchange(fields); + output::Gmsh gmsh(grid.name() + "_fields.msh", "a"); + gmsh.write(fields["lapl"]); + } + + double min, max, mean; + idx_t N; + fvm.node_columns().minimum(fields["lapl"], min); + fvm.node_columns().maximum(fields["lapl"], max); + fvm.node_columns().mean(fields["lapl"], mean, N); + print_min_max_mean("lapl"); + EXPECT_APPROX_EQ(min, -6.4088005677811607095e-13, tolerance); + EXPECT_APPROX_EQ(max, 9.8984499569639476135e-12, tolerance); + EXPECT_APPROX_EQ(mean, -1.03409e-13, tolerance); + }; + SECTION("double precision") { + do_test(double{}, 1.e-16); + } + SECTION("single precision") { + do_test(float{}, 1.e-10); } - double min, max, mean; - idx_t N; - fvm.node_columns().minimum(fields["lapl"], min); - fvm.node_columns().maximum(fields["lapl"], max); - fvm.node_columns().mean(fields["lapl"], mean, N); - print_min_max_mean("lapl"); - EXPECT_APPROX_EQ(min, -6.4088005677811607095e-13); - EXPECT_APPROX_EQ(max, 9.8984499569639476135e-12); - EXPECT_APPROX_EQ(mean, -1.03409e-13); } //-----------------------------------------------------------------------------