From a0b9361a8a9f5fd4a8193c6b9915954c57894227 Mon Sep 17 00:00:00 2001 From: Philipp Windischhofer Date: Thu, 25 Jan 2024 17:28:44 -0600 Subject: [PATCH] moving WeightingFieldUtils to use DistributedWeightingField --- include/Eisvogel/WeightingFieldUtils.hh | 5 +- src/WeightingFieldUtils.cxx | 218 +++++++++++------------- 2 files changed, 104 insertions(+), 119 deletions(-) diff --git a/include/Eisvogel/WeightingFieldUtils.hh b/include/Eisvogel/WeightingFieldUtils.hh index 532d6f03f..8e5e06b6b 100644 --- a/include/Eisvogel/WeightingFieldUtils.hh +++ b/include/Eisvogel/WeightingFieldUtils.hh @@ -10,8 +10,9 @@ namespace WeightingFieldUtils { const CoordVector& start_coords, const CoordVector& end_coords, scalar_t tp, unsigned int N, scalar_t r_min, scalar_t os_factor = 1.5, scalar_t n = 1); - WeightingField SampleElectricDipoleWeightingField(const CoordVector& start_coords, const CoordVector& end_coords, - scalar_t tp, unsigned int N, scalar_t r_min, scalar_t os_factor, scalar_t n = 1); + void SampleElectricDipoleWeightingFieldChunk(ScalarField3D& E_r_buffer, ScalarField3D& E_z_buffer, ScalarField3D& E_phi_buffer, + const CoordVector& start_coords, const CoordVector& end_coords, + scalar_t tp, unsigned int N, scalar_t r_min, scalar_t os_factor, scalar_t n); } diff --git a/src/WeightingFieldUtils.cxx b/src/WeightingFieldUtils.cxx index 776fe9482..48aa3746e 100644 --- a/src/WeightingFieldUtils.cxx +++ b/src/WeightingFieldUtils.cxx @@ -4,6 +4,7 @@ #include "Eisvogel/CoordUtils.hh" #include "Eisvogel/MathUtils.hh" #include "Eisvogel/Serialization.hh" +#include "Eisvogel/DistributedWeightingField.hh" #include #include #include @@ -15,103 +16,11 @@ namespace WeightingFieldUtils { void CreateElectricDipoleWeightingField(std::string wf_path, const CoordVector& start_coords, const CoordVector& end_coords, scalar_t tp, unsigned int N, scalar_t r_min, scalar_t os_factor, scalar_t n) { - - std::fstream ofs; - ofs.open(wf_path, std::ios::out | std::ios::binary); - stor::Serializer oser(ofs); - - // TODO: later, this will also break up the weighting field into multiple chunks just like for meep - std::cout << "Building weighting field ..." << std::endl; - WeightingField wf_out = SampleElectricDipoleWeightingField(start_coords, end_coords, tp, N, r_min, os_factor, n); - std::cout << "Saving weighting field ..." << std::endl; - oser.serialize(wf_out); - ofs.close(); - } - - scalar_t filtered_theta(scalar_t t, scalar_t tp, unsigned int N) { - if(t <= 0) { - return 0.0; - } - return 1.0 - MathUtils::incomplete_gamma(1 + N, N * t / tp) / std::exp(std::lgamma(N + 1)); - }; - - scalar_t filtered_delta(scalar_t t, scalar_t tp, unsigned int N) { - if(t <= 0) { - return 0.0; - } - return std::pow(t / tp * N, N) * std::exp(-t / tp * N) / (tp * std::exp(std::lgamma(N))); - }; - - scalar_t filtered_delta_prime(scalar_t t, scalar_t tp, unsigned int N) { - if(t <= 0) { - return 0.0; - } - return filtered_delta(t, tp, N) * (tp - t) * N / (tp * t); - }; - - // Weighting field in spherical coordinates - scalar_t E_r(scalar_t t, scalar_t r_xy, scalar_t z, scalar_t ior, scalar_t c, scalar_t r_min, scalar_t Qw, scalar_t ds, scalar_t eps0, scalar_t N, scalar_t tp) { - - r_xy = std::fabs(r_xy); - scalar_t r = std::sqrt(std::pow(r_xy, 2) + std::pow(z, 2)); - if(r < r_min) { - return std::nan(""); - } - scalar_t t_prop = r * ior / c, t_del = t - t_prop; - scalar_t cos_theta = z / r; - - return -2.0 * Qw * ds / (eps0 * 4 * M_PI) * cos_theta / std::pow(r, 3) * (filtered_theta(t_del, tp, N) + - t_prop * filtered_delta(t_del, tp, N)); - }; - - scalar_t E_theta(scalar_t t, scalar_t r_xy, scalar_t z, scalar_t ior, scalar_t c, scalar_t r_min, scalar_t Qw, scalar_t ds, scalar_t eps0, scalar_t N, scalar_t tp) { - - r_xy = std::fabs(r_xy); - scalar_t r = std::sqrt(std::pow(r_xy, 2) + std::pow(z, 2)); - if(r < r_min) { - return std::nan(""); - } - scalar_t t_prop = r * ior / c, t_del = t - t_prop; - scalar_t sin_theta = r_xy / r; - return -Qw * ds / (eps0 * 4 * M_PI) * sin_theta / std::pow(r, 3) * (filtered_theta(t_del, tp, N) + t_prop * filtered_delta(t_del, tp, N) - + std::pow(t_prop, 2) * filtered_delta_prime(t_del, tp, N)); - }; - - // Weighting field in cylindrical coordinates - scalar_t E_rxy(scalar_t t, scalar_t r_xy, scalar_t z, scalar_t ior, scalar_t c, scalar_t r_min, scalar_t Qw, scalar_t ds, scalar_t eps0, scalar_t N, scalar_t tp) { - r_xy = std::fabs(r_xy); - scalar_t r = std::sqrt(std::pow(r_xy, 2) + std::pow(z, 2)); - scalar_t cos_theta = z / r, sin_theta = r_xy / r; - return E_r(t, r_xy, z, ior, c, r_min, Qw, ds, eps0, N, tp) * sin_theta + E_theta(t, r_xy, z, ior, c, r_min, Qw, ds, eps0, N, tp) * cos_theta; - }; - - scalar_t E_z(scalar_t t, scalar_t r_xy, scalar_t z, scalar_t ior, scalar_t c, scalar_t r_min, scalar_t Qw, scalar_t ds, scalar_t eps0, scalar_t N, scalar_t tp) { - r_xy = std::fabs(r_xy); - scalar_t r = std::sqrt(std::pow(r_xy, 2) + std::pow(z, 2)); - scalar_t cos_theta = z / r, sin_theta = r_xy / r; - return E_r(t, r_xy, z, ior, c, r_min, Qw, ds, eps0, N, tp) * cos_theta - E_theta(t, r_xy, z, ior, c, r_min, Qw, ds, eps0, N, tp) * sin_theta; - }; - - scalar_t E_phi(scalar_t t, scalar_t r_xy, scalar_t z) { - return 0.0; - }; - - ScalarField3D SampleElectricDipoleWeightingField_E_r(const CoordVector& start_coords, const CoordVector& end_coords, - scalar_t tp, unsigned int N, scalar_t delta_t, scalar_t delta_pos, - scalar_t os_factor, scalar_t ior) { - - } - - WeightingField SampleElectricDipoleWeightingField(const CoordVector& start_coords, const CoordVector& end_coords, - scalar_t tp, unsigned int N, scalar_t r_min, scalar_t os_factor, scalar_t n) { - - scalar_t Qw = 1.0; - scalar_t eps0 = 1.0; // vacuum dielectric constant - scalar_t ds = 1.0; - scalar_t c = 1.0; // speed of light in vacuum + DistributedWeightingField dwf(wf_path, start_coords, end_coords); // compute required step size for sampling of weighting field + scalar_t c = 1.0; // speed of light in vacuum scalar_t fmax = (scalar_t)N / (2 * M_PI * tp) * std::sqrt(std::pow(2.0, 1.0 / (N + 1)) - 1); scalar_t lambda_min = c / (fmax * n); scalar_t delta_t = 1.0 / (2 * fmax * os_factor); @@ -122,45 +31,120 @@ namespace WeightingFieldUtils { std::cout << "delta_t = " << delta_t << std::endl; std::cout << "delta_r = delta_z = " << delta_pos << std::endl; std::cout << "start_coords: t = " << C::getT(start_coords) << ", r = " << C::getR(start_coords) << ", z = " << C::getZ(start_coords) << std::endl; - std::cout << "end_coords: t = " << C::getT(end_coords) << ", r = " << C::getR(end_coords) << ", z = " << C::getZ(end_coords) << std::endl; std::cout << "---------------------------" << std::endl; - DeltaVector step_requested = C::MakeCoordVectorTRZ(delta_t, delta_pos, delta_pos); + DeltaVector stepsize_requested = C::MakeCoordVectorTRZ(delta_t, delta_pos, delta_pos); - // ============== - - CoordVector number_pts = (end_coords - start_coords) / step_requested; + CoordVector number_pts = (end_coords - start_coords) / stepsize_requested; std::size_t pts_t = C::getT(number_pts), pts_r = C::getR(number_pts), pts_z = C::getZ(number_pts); - // TODO: find a better way to make sure the ordering t, z, r etc is not messed up - ScalarField3D E_r_sampled({pts_t, pts_z, pts_r}, 0.0); - ScalarField3D E_z_sampled({pts_t, pts_z, pts_r}, 0.0); - ScalarField3D E_phi_sampled({pts_t, pts_z, pts_r}, 0.0); - - DeltaVector step = (end_coords - start_coords) / E_r_sampled.shape(); + // TODO: generalize this to produce more than one chunk + + ScalarField3D chunk_buffer_E_r({pts_t, pts_z, pts_r}, 0.0); + ScalarField3D chunk_buffer_E_z({pts_t, pts_z, pts_r}, 0.0); + ScalarField3D chunk_buffer_E_phi({pts_t, pts_z, pts_r}, 0.0); + + // Fill chunk buffers + SampleElectricDipoleWeightingFieldChunk(chunk_buffer_E_r, chunk_buffer_E_z, chunk_buffer_E_phi, start_coords, end_coords, tp, N, r_min, os_factor, n); + + // Register chunk buffers + dwf.RegisterChunk(chunk_buffer_E_r, chunk_buffer_E_z, chunk_buffer_E_phi, IndexVector({0, 0, 0})); + + dwf.Flush(); + } + + // TODO: three separate `ScalarField3D`s to be replaced with single vector field + void SampleElectricDipoleWeightingFieldChunk(ScalarField3D& E_r_buffer, ScalarField3D& E_z_buffer, ScalarField3D& E_phi_buffer, + const CoordVector& start_coords, const CoordVector& end_coords, + scalar_t tp, unsigned int N, scalar_t r_min, scalar_t os_factor, scalar_t n) { + + scalar_t Qw = 1.0; + scalar_t eps0 = 1.0; // vacuum dielectric constant + scalar_t ds = 1.0; + scalar_t c = 1.0; // speed of light in vacuum + + auto filtered_theta = [&](scalar_t t) -> scalar_t { + if(t <= 0) { + return 0.0; + } + return 1.0 - MathUtils::incomplete_gamma(1 + N, N * t / tp) / std::exp(std::lgamma(N + 1)); + }; + + auto filtered_delta = [&](scalar_t t) -> scalar_t { + if(t <= 0) { + return 0.0; + } + return std::pow(t / tp * N, N) * std::exp(-t / tp * N) / (tp * std::exp(std::lgamma(N))); + }; + + auto filtered_delta_prime = [&](scalar_t t) -> scalar_t { + if(t <= 0) { + return 0.0; + } + return filtered_delta(t) * (tp - t) * N / (tp * t); + }; + + // Weighting field in spherical coordinates + auto E_r = [&](scalar_t t, scalar_t r_xy, scalar_t z) -> scalar_t { + r_xy = std::fabs(r_xy); + scalar_t r = std::sqrt(std::pow(r_xy, 2) + std::pow(z, 2)); + if(r < r_min) { + return std::nan(""); + } + scalar_t t_prop = r * n / c, t_del = t - t_prop; + scalar_t cos_theta = z / r; + + return -2.0 * Qw * ds / (eps0 * 4 * M_PI) * cos_theta / std::pow(r, 3) * (filtered_theta(t_del) + + t_prop * filtered_delta(t_del)); + }; + + auto E_theta = [&](scalar_t t, scalar_t r_xy, scalar_t z) -> scalar_t { + r_xy = std::fabs(r_xy); + scalar_t r = std::sqrt(std::pow(r_xy, 2) + std::pow(z, 2)); + if(r < r_min) { + return std::nan(""); + } + scalar_t t_prop = r * n / c, t_del = t - t_prop; + scalar_t sin_theta = r_xy / r; + return -Qw * ds / (eps0 * 4 * M_PI) * sin_theta / std::pow(r, 3) * (filtered_theta(t_del) + t_prop * filtered_delta(t_del) + + std::pow(t_prop, 2) * filtered_delta_prime(t_del)); + }; + + // Weighting field in cylindrical coordinates + auto E_rxy = [&](scalar_t t, scalar_t r_xy, scalar_t z) -> scalar_t { + r_xy = std::fabs(r_xy); + scalar_t r = std::sqrt(std::pow(r_xy, 2) + std::pow(z, 2)); + scalar_t cos_theta = z / r, sin_theta = r_xy / r; + return E_r(t, r_xy, z) * sin_theta + E_theta(t, r_xy, z) * cos_theta; + }; + + auto E_z = [&](scalar_t t, scalar_t r_xy, scalar_t z) -> scalar_t { + r_xy = std::fabs(r_xy); + scalar_t r = std::sqrt(std::pow(r_xy, 2) + std::pow(z, 2)); + scalar_t cos_theta = z / r, sin_theta = r_xy / r; + return E_r(t, r_xy, z) * cos_theta - E_theta(t, r_xy, z) * sin_theta; + }; + + auto E_phi = [&](scalar_t t, scalar_t r_xy, scalar_t z) -> scalar_t { + return 0.0; + }; IndexVector start_inds({0, 0, 0}); - IndexVector end_inds({pts_t, pts_z, pts_r}); + IndexVector end_inds(E_r_buffer.shape()); for(IndexCounter cnt(start_inds, end_inds); cnt.running(); ++cnt) { IndexVector ind = cnt.index(); - CoordVector coords = WeightingField::FracIndsToCoord(ind, start_coords, end_coords, E_r_sampled.shape()); + CoordVector coords = WeightingField::FracIndsToCoord(ind, start_coords, end_coords, E_r_buffer.shape()); scalar_t t = C::getT(coords); scalar_t r = C::getR(coords); scalar_t z = C::getZ(coords); - scalar_t cur_E_rxy = E_rxy(t, r, z, n, c, r_min, Qw, ds, eps0, N, tp); - scalar_t cur_E_z = E_z(t, r, z, n, c, r_min, Qw, ds, eps0, N, tp); - scalar_t cur_E_phi = E_phi(t, r, z); - - E_r_sampled(ind) = cur_E_rxy; - E_z_sampled(ind) = cur_E_z; - E_phi_sampled(ind) = cur_E_phi; + E_r_buffer(ind) = E_rxy(t, r, z); + E_z_buffer(ind) = E_z(t, r, z); + E_phi_buffer(ind) = E_phi(t, r, z); + } - - return WeightingField(std::move(E_r_sampled), std::move(E_z_sampled), std::move(E_phi_sampled), - std::move(start_coords), std::move(end_coords)); } }