Skip to content

Commit

Permalink
moving WeightingFieldUtils to use DistributedWeightingField
Browse files Browse the repository at this point in the history
  • Loading branch information
@philippwindischhofer
philippwindischhofer committed Jan 25, 2024
1 parent cf2e4b2 commit a0b9361
Show file tree
Hide file tree
Showing 2 changed files with 104 additions and 119 deletions.
5 changes: 3 additions & 2 deletions include/Eisvogel/WeightingFieldUtils.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<scalar_t>& E_r_buffer, ScalarField3D<scalar_t>& E_z_buffer, ScalarField3D<scalar_t>& 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);

}

Expand Down
218 changes: 101 additions & 117 deletions src/WeightingFieldUtils.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@
#include "Eisvogel/CoordUtils.hh"
#include "Eisvogel/MathUtils.hh"
#include "Eisvogel/Serialization.hh"
#include "Eisvogel/DistributedWeightingField.hh"
#include <cmath>
#include <iostream>
#include <fstream>
Expand All @@ -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<scalar_t> 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);
Expand All @@ -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<scalar_t> E_r_sampled({pts_t, pts_z, pts_r}, 0.0);
ScalarField3D<scalar_t> E_z_sampled({pts_t, pts_z, pts_r}, 0.0);
ScalarField3D<scalar_t> 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<scalar_t> chunk_buffer_E_r({pts_t, pts_z, pts_r}, 0.0);
ScalarField3D<scalar_t> chunk_buffer_E_z({pts_t, pts_z, pts_r}, 0.0);
ScalarField3D<scalar_t> 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<scalar_t>& E_r_buffer, ScalarField3D<scalar_t>& E_z_buffer, ScalarField3D<scalar_t>& 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));
}
}

0 comments on commit a0b9361

Please sign in to comment.