Make some changes

This commit is contained in:
2024-01-01 15:57:03 +01:00
parent 943f3b305b
commit 3cc82dd959
2 changed files with 63 additions and 5 deletions

View File

@@ -27,10 +27,10 @@ void WaveSimulation::initialize_U(double x_c, double y_c, double sigma_x,
double x, y, diff_x, diff_y;
std::complex<double> sum = 0.;
for (size_t j = 0; j < this->U.n_cols; j++) {
x = j * h;
x = (j+1) * h;
diff_x = x - x_c;
for (size_t i = 0; i < this->U.n_rows; i++) {
y = i * h;
y = (i+1) * h;
diff_y = y - y_c;
this->U(i, j) =
std::exp(-(diff_x * diff_x) / (2. * sigma_x * sigma_x)
@@ -117,7 +117,8 @@ void WaveSimulation::initialize_V(double thickness, double pos_x, double ap_sep,
res = arma::join_cols(res, arma::cx_vec(fill + ((this->N - res.size()) % 2),
arma::fill::value(1e10)));
uint32_t start = pos_x / this->h - thickness / this->h / 2;
// Subtract one to compensate for the border being gone
uint32_t start = (pos_x / this->h - 1) - thickness / this->h / 2;
for (size_t i = 0; i < thickness / this->h; i++) {
this->V.col(start + i) = res;
}
@@ -125,32 +126,40 @@ void WaveSimulation::initialize_V(double thickness, double pos_x, double ap_sep,
void WaveSimulation::initialize_A()
{
DEBUG("Inside initialize_A");
// Create the diagonal
arma::cx_vec diagonal(this->N * this->N);
DEBUG("Before setting diagonal values");
// Set diagonal values
std::complex<double> r = (1._i * this->dt) / (2 * h * h);
for (size_t i = 0; i < diagonal.size(); i++) {
diagonal(i) = 1. + 4. * r + (1._i * this->dt / 2.) * this->V(i);
}
DEBUG("Create submatrix");
// Create the submatrix
arma::cx_mat sub_matrix(this->N, this->N, arma::fill::zeros);
sub_matrix.diag(-1).fill(-r);
sub_matrix.diag(1).fill(-r);
DEBUG("Set the size of A");
// Set the size of A
this->A.set_size(this->N * this->N, this->N * this->N);
DEBUG("Fill in values");
// Fill in the values in the submatrix diagonal
for (size_t i = 0; i < this->A.n_cols; i += this->N) {
this->A.submat(i, i, i + this->N - 1, i + this->N - 1) = sub_matrix;
}
DEBUG("Fill in diagonals");
// Fill the last sub/sup-diagonals
this->A.diag() = diagonal;
this->A.diag(-this->N).fill(-r);
DEBUG("After diagonal");
this->A.diag(this->N).fill(-r);
DEBUG("After");
}
void WaveSimulation::initialize_B()
@@ -205,15 +214,21 @@ WaveSimulation::WaveSimulation(double h, double dt, double T, double x_c,
double y_c, double sigma_x, double sigma_y,
double p_x, double p_y)
{
DEBUG("Inside constructor");
this->dt = dt;
this->h = h;
this->T = T;
this->M = 1. / h;
this->N = M - 2;
DEBUG("Before V");
this->initialize_V();
DEBUG("Before U");
this->initialize_U(x_c, y_c, sigma_x, sigma_y, p_x, p_y);
DEBUG("Before A");
this->initialize_A();
DEBUG("Before B");
this->initialize_B();
DEBUG("After B");
}
// Public methods
@@ -226,7 +241,7 @@ void WaveSimulation::step()
arma::spsolve(this->U, this->A, tmp);
}
void WaveSimulation::solve(std::string outfile, bool write_each_step)
void WaveSimulation::simulate(std::string outfile, bool write_each_step)
{
// Create path and proceed if successful.
if (!utils::mkpath(utils::dirname(outfile))) {
@@ -256,7 +271,7 @@ void WaveSimulation::solve(std::string outfile, bool write_each_step)
this->write_U(ofile);
ofile.close();
}
void WaveSimulation::solve(std::string outfile, std::vector<double> &steps)
void WaveSimulation::simulate(std::string outfile, std::vector<double> &steps)
{
// Create path and proceed if successful.
if (!utils::mkpath(utils::dirname(outfile))) {