Skip to content

Commit

Permalink
.
Browse files Browse the repository at this point in the history
  • Loading branch information
maksimKravchenko1 committed Jun 1, 2023
1 parent beaa522 commit 619e421
Show file tree
Hide file tree
Showing 7 changed files with 171 additions and 91 deletions.
30 changes: 15 additions & 15 deletions Wz.txt

Large diffs are not rendered by default.

17 changes: 11 additions & 6 deletions dist/test-addon.js
Original file line number Diff line number Diff line change
Expand Up @@ -32,10 +32,15 @@ var test1D = function () {
var test2D = function () {
var lambda = 1;
var beamsize = 1;
var materialVector = [1, 0, 2, 0, 1, 1, 1, 1, 1];
var eps = [1.0, 1.2, 1.1];
var mu = [0.51, 0.5, 0.57];
var sigma = [1.0, 0.001, 1.0];
// const materialVector = [1, 0, 2, 0, 1, 1, 1,1 ,1];
// const materialVector = [1, 0, 1, 0, 1, 1, 1, 1 ,1];
var materialVector = [0, 0, 0, 0];
// const eps = [1.0, 1.2, 1.1];
// const mu = [0.51, 0.5, 0.57];
// const sigma = [1.0, 0.001, 1.0];
var eps = [1.2, 1.15];
var mu = [1, 1];
var sigma = [0, 1e+7];
var srcPosition = [0.1, 0.1];
var dataReturnType = 0;
var isReload = true;
Expand Down Expand Up @@ -70,6 +75,6 @@ var test2DTFSF = function () {
console.log(data);
};
// test1D();
test2D();
// test2DTFSF();
// test2D();
test2DTFSF();
// testMemoryUsage();
2 changes: 1 addition & 1 deletion src/fdtd/2d-pml copy/fdtd-pml-2d.h
Original file line number Diff line number Diff line change
Expand Up @@ -97,7 +97,7 @@ class FdtdPml2D {

size_t GetStep();
size_t GetCurrentTimeStep();

sqrt
void SetParams(std::vector<std::vector<double>> &eps,
std::vector<std::vector<double>> &mu,
std::vector<std::vector<double>> &sigma,
Expand Down
131 changes: 87 additions & 44 deletions src/fdtd/2d-pml/fdtd-pml-2d.h
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@
#include <math.h>

#define MEDIACONSTANT (2)
#define NUMBEROFITERATIONCONSTANT (1400) // was 300
#define NUMBEROFITERATIONCONSTANT (1000) // was 300
#define BIGLINESIZE (8192)


Expand All @@ -46,7 +46,7 @@ class FdtdPml2D {
double rmax, orderbc;
static const size_t ie = 220; //number of grid cells in x-direction
static const size_t je = 220; //number of grid cells in y-direction
int ib, jb, is, js, nmax, iebc, jebc, ibbc, jbbc, iefbc, jefbc, ibfbc, jbfbc ;
int ib, jb, is, js, nmax, iebc, jebc, ibbc, jbbc, iefbc, jefbc, ibfbc, jbfbc;
int media;
double rtau, tau, delay;
double **ex; //fields in main grid
Expand Down Expand Up @@ -74,22 +74,22 @@ class FdtdPml2D {
double **cbey;
double **dahz;
double **dbhz;
double **caexbcf ; // pml coefficients
double **cbexbcf ;
double **caexbcb ;
double **cbexbcb ;
double **caexbcl ;
double **cbexbcl ;
double **caexbcr ;
double **cbexbcr ;
double **caeybcf ;
double **cbeybcf ;
double **caeybcb ;
double **cbeybcb ;
double **caeybcl ;
double **cbeybcl ;
double **caeybcr ;
double **cbeybcr ;
double **caexbcf; // pml coefficients
double **cbexbcf;
double **caexbcb;
double **cbexbcb;
double **caexbcl;
double **cbexbcl;
double **caexbcr;
double **cbexbcr;
double **caeybcf;
double **cbeybcf;
double **caeybcb;
double **cbeybcb;
double **caeybcl;
double **cbeybcl;
double **caeybcr;
double **cbeybcr;
double **dahzxbcf;
double **dbhzxbcf;
double **dahzxbcb;
Expand All @@ -110,10 +110,10 @@ class FdtdPml2D {
double ca1,cb1,y1,y2, sigmay,sigmays,da1,db1;
double x1,x2,sigmax,sigmaxs;
double minimumValue, maximumValue;
char filename[BIGLINESIZE] ;
FILE *filePointer ;
char filename[BIGLINESIZE];
FILE *filePointer;
double scaleValue;
int iValue,plottingInterval,centery,centerx ;
int iValue,plottingInterval,centery,centerx;



Expand Down Expand Up @@ -149,10 +149,18 @@ class FdtdPml2D {
output.cols = je;

// size_t flatten_array_size = ie * nyje
double min = -0.000001;
double max = 0.000001;
for(size_t i = 0; i < ie; i += 1) {
for(size_t j = 0; j < je; j += 1) {
// output.Ez[i*je + j] = Ez[i][j];
output.Hz[i*je + j] = hz[i][j];
if(hz[i][j] > max) {
max = hz[i][j];
}
if(hz[i][j] < min) {
min = hz[i][j];
}
output.X[i*je + j] = i;
output.Y[i*je + j] = j;
}
Expand All @@ -161,8 +169,10 @@ class FdtdPml2D {
// output.maxEz = *std::max_element(std::begin(output.Ez), std::end(output.Ez));
// output.minEz = *std::min_element(std::begin(output.Ez), std::end(output.Ez));

output.maxHz = *std::max_element(std::begin(output.Hz), std::end(output.Hz));
output.minHz = *std::min_element(std::begin(output.Hz), std::end(output.Hz));
// output.maxHz = *std::max_element(std::begin(output.Hz), std::end(output.Hz));
// output.minHz = *std::min_element(std::begin(output.Hz), std::end(output.Hz));
output.maxHz = max;
output.minHz = min;

return output;
}
Expand Down Expand Up @@ -342,8 +352,8 @@ class FdtdPml2D {
const size_t tau = 20;
double src = -2.0 * ((n - t0) / tau) * std::exp(-1.0 * std::pow((n - t0) / tau, 2));

hz[is][js] = src;
// hz[is][js] = source[n];
// hz[is][js] = src;
hz[is][js] = source[n];


//***********************************************************************
Expand Down Expand Up @@ -469,11 +479,18 @@ double **AllocateMemory (int imax, int jmax, double initialValue)
pointer[i][j] = initialValue;
} /* jForLoop */
} /* iForLoop */
return(pointer) ;
return(pointer);
}


void InitializeFdtd ()
// template <size_t ie, size_t je>
// template <typename TwoD>
// void InitializeFdtd (TwoD& material_matrix)
void InitializeFdtd (std::vector<std::vector<int>> &material_matrix, std::vector<double> &eps,
std::vector<double> &mur,
std::vector<double> &sig, int src_position_row, int src_position_col)
// void InitializeFdtd (size_t (&material_matrix)[ie][je])

// void InitializeFdtd ()
{


Expand All @@ -484,8 +501,8 @@ double **AllocateMemory (int imax, int jmax, double initialValue)
minimumValue = -0.1;
maximumValue = 0.1;
plottingInterval = 0;
centery = 25 ;
centerx = 15 ;
centery = 25;
centerx = 15;

//***********************************************************************
// Fundamental constants
Expand All @@ -500,6 +517,17 @@ double **AllocateMemory (int imax, int jmax, double initialValue)
omega = 2.0 * pi * freq; //center frequency in radians


double eta1 = sqrt(mur[0]*muz/(eps[0]*epsz));
double eta2 = sqrt(mur[1]*muz/(eps[1]*epsz));
double reflection_coef = abs((eta2-eta1) / (eta2 + eta1));
double transmittance_coef = 1-reflection_coef;

std::cout << "eta1: " << eta1 << std::endl;
std::cout << "eta2: " << eta2 << std::endl;
std::cout << "reflection_coef: " << reflection_coef << std::endl;
std::cout << "transmittance_coef: " << 1-abs(reflection_coef) << std::endl;


//***********************************************************************
// Grid parameters
//***********************************************************************
Expand All @@ -510,8 +538,12 @@ double **AllocateMemory (int imax, int jmax, double initialValue)
ib = ie + 1; // one extra is needed for fields on the boundaries (ie Ex on top boundary, Ey on right boundary)
jb = je + 1; // ditto

is = 15; //location of z-directed hard source
js = je / 2; //location of z-directed hard source
// is = 15; //location of z-directed hard source
// js = je / 2; //location of z-directed hard source

// swap later
js = src_position_row; //location of z-directed hard source
is = src_position_col; //location of z-directed hard source

dx = 3.0e-3; //space increment of square lattice (meters)
dt = dx / (2.0 * cc); //time step, seconds, courant limit, Taflove1995 page 177
Expand Down Expand Up @@ -621,19 +653,30 @@ double **AllocateMemory (int imax, int jmax, double initialValue)
jcenter = je / 2; // j-coordinate of cylinder's center
for (i = 0; i < ie; i++) {
for (j = 0; j < je; j++) {
temporaryi = (double )(i - icenter) ;
temporaryj = (double )(j - jcenter) ;
temporaryi = (double )(i - icenter);
temporaryj = (double )(j - jcenter);
dist2 = (temporaryi + 0.5) * (temporaryi + 0.5) + (temporaryj) * (temporaryj);
if (dist2 <= (rad * rad)) {
caex[i][j] = ca[1];
cbex[i][j] = cb[1];
} /* if */
// This looks tricky! Why can't caey/cbey use the same 'if' statement as caex/cbex above ??
dist2 = (temporaryj + 0.5) * (temporaryj + 0.5) + (temporaryi) * (temporaryi);
if (dist2 <= (rad * rad)) {
caey[i][j] = ca[1];
cbey[i][j] = cb[1];
} /* if */

caex[i][j] = ca[material_matrix[i][j]];
cbex[i][j] = cb[material_matrix[i][j]];

caey[i][j] = ca[material_matrix[i][j]];
cbey[i][j] = cb[material_matrix[i][j]];

// std::cout << i*je + j << ": " << ca[material_matrix[i][j]] << std::endl;


// if (dist2 <= (rad * rad)) {
// caex[i][j] = ca[1];
// cbex[i][j] = cb[1];
// } /* if */

// // This looks tricky! Why can't caey/cbey use the same 'if' statement as caex/cbex above ??
// dist2 = (temporaryj + 0.5) * (temporaryj + 0.5) + (temporaryi) * (temporaryi);
// if (dist2 <= (rad * rad)) {
// caey[i][j] = ca[1];
// cbey[i][j] = cb[1];
// } /* if */
} /* jForLoop */
} /* iForLoop */

Expand Down
22 changes: 11 additions & 11 deletions src/fdtd/2d-upml-tf-sf/fdtd-2d-upml-tf-sf.h
Original file line number Diff line number Diff line change
Expand Up @@ -110,20 +110,20 @@ class TFSF {
// static constexpr size_t ny_a = 6;
// static constexpr size_t ny_b = 24;

// static constexpr size_t nx_a = 60;
// static constexpr size_t nx_b = 160;
// static constexpr size_t ny_a = 60;
// static constexpr size_t ny_b = 160;
static constexpr size_t nx_a = 60;
static constexpr size_t nx_b = 160;
static constexpr size_t ny_a = 60;
static constexpr size_t ny_b = 160;

// static constexpr size_t nx_a = 15;
// static constexpr size_t nx_b = 205;
// static constexpr size_t ny_a = 15;
// static constexpr size_t ny_b = 205;

static constexpr size_t nx_a = 10;
static constexpr size_t nx_b = 209;
static constexpr size_t ny_a = 10;
static constexpr size_t ny_b = 209;
// static constexpr size_t nx_a = 10;
// static constexpr size_t nx_b = 209;
// static constexpr size_t ny_a = 10;
// static constexpr size_t ny_b = 209;

// Pre-allocate 1D fields for TF/SF interface
// TM components
Expand Down Expand Up @@ -275,9 +275,9 @@ class TFSF {
// }};

std::array<std::array<size_t, 4>, 3> mat_conf = {{
{100, 25, 10, 85},
{100, 25, 105, 10},
{100, 25, 125, 80},
{100, 25, 80, 50},
{100, 25, 80, 50},
{100, 25, 80, 50},
// {120, 10, 100, 20},
// {140, 10, 100, 20}

Expand Down
43 changes: 35 additions & 8 deletions src/transform-to-js/fdtd-2d/fdtd-2d.cpp
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
#include <iostream>
#include <cmath>
#include <array>
#include "fdtd-2d.h"

Napi::Object Fdtd2D::Init(Napi::Env env, Napi::Object exports) {
Expand Down Expand Up @@ -157,7 +158,17 @@ Fdtd2D::Fdtd2D(const Napi::CallbackInfo& info)
// std::cout << "before ";// << material_matrix_size;

// Temporary matrix.

// const size_t rows = FdtdPml2D::GetRows();
// const size_t cols = FdtdPml2D::GetCols();
// const size_t nx = 220;
// static const size_t ny = 220;
static const size_t rows = 220;
static const size_t cols = 220;

// std::array<std::array<size_t, cols>, rows> temp_matrix;
std::vector<std::vector<int>> temp_matrix;
std::vector<std::vector<int>> temp_matrix_2;

// std::cout << "--" << material_matrix_size << "--";

Expand All @@ -169,20 +180,27 @@ Fdtd2D::Fdtd2D(const Napi::CallbackInfo& info)
temp_matrix[i].push_back(
(int)material_matrix_js[i * material_matrix_size + j]
.As<Napi::Number>());
// std::cout << ", " << i;
// temp_matrix[i][j] =
// (int)material_matrix_js[i * material_matrix_size + j]
// .As<Napi::Number>();
}
}

// std::cout << "after";

// const size_t rows = FdtdPml2D::GetRows();
// const size_t cols = FdtdPml2D::GetCols();
const size_t rows = 220;
const size_t cols = 220;

// Matrix size coefficient.
size_t coeff = rows / material_matrix_size;


std::vector<double> eps;
std::vector<double> mu;
std::vector<double> sigma;

for( size_t i = 0; i < eps_js.Length(); ++i) {
eps.push_back(static_cast<double>(eps_js[i].As<Napi::Number>()));
mu.push_back(static_cast<double>(mu_js[i].As<Napi::Number>()));
sigma.push_back(static_cast<double>(sigma_js[i].As<Napi::Number>()));
}

// Initialization eps, mu, sigma matrixes.
std::vector<std::vector<double>> eps_matrix;
std::vector<std::vector<double>> mu_matrix;
Expand All @@ -191,10 +209,15 @@ Fdtd2D::Fdtd2D(const Napi::CallbackInfo& info)
eps_matrix.push_back(std::vector<double>());
mu_matrix.push_back(std::vector<double>());
sigma_matrix.push_back(std::vector<double>());

temp_matrix_2.push_back(std::vector<int>());

for (int j = 0; j < cols; j++) {
eps_matrix[i].push_back(0);
mu_matrix[i].push_back(0);
sigma_matrix[i].push_back(0);

temp_matrix_2[i].push_back(0);
}
}

Expand All @@ -205,6 +228,9 @@ Fdtd2D::Fdtd2D(const Napi::CallbackInfo& info)
for (int f = 0; f < coeff; f++) {
int index = temp_matrix[i][j];

// temp_matrix_2[i * coeff + k][j * coeff + f] = index;
temp_matrix_2[j * coeff + f][i * coeff + k] = index;

// Rotate matrix on 90 degree for correctness in numerical method.
// eps_matrix[j * coeff + f][i * coeff + k] = static_cast<double>(eps_js[index].As<Napi::Number>());
// mu_matrix[j * coeff + f][i * coeff + k] = static_cast<double>(mu_js[index].As<Napi::Number>());
Expand Down Expand Up @@ -237,7 +263,8 @@ Fdtd2D::Fdtd2D(const Napi::CallbackInfo& info)
size_t src_position_col = static_cast<size_t>(relative_src_position_x * cols);

// fdtd.SetParams(eps_matrix, mu_matrix, sigma_matrix, src_position_row, src_position_col);
fdtd.InitializeFdtd();
// fdtd.InitializeFdtd();
fdtd.InitializeFdtd(temp_matrix_2, eps, mu, sigma, src_position_row, src_position_col);
}

// Fdtd method in 1D case.
Expand Down
Loading

0 comments on commit 619e421

Please sign in to comment.