max/tug
1
0
Fork 0
mirror of https://git.gfz-potsdam.de/naaice/tug.git synced 2025-12-13 17:38:23 +01:00
Code Issues Packages Projects Releases Wiki Activity

Resolved merge conflict (HEAD taken)

Browse Source
This commit is contained in:
Marco De Lucia 2022-04-26 13:38:14 +02:00
commit 722ed240eb
1 changed files with 24 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

40
src/BTCSDiffusion.cpp
View File

@ -141,11 +141,10 @@ void Diffusion::BTCSDiffusion::simulate2D(
int n_rows = this->grid_cells[1];
int n_cols = this->grid_cells[0];
double dx = this->deltas[0];
DMatrixRowMajor t0_c;
double local_dt = this->time_step / BTCS_2D_DT_SIZE;
t0_c = calc_t0_c(c, alpha, bc, local_dt, dx);
DMatrixRowMajor t0_c = calc_t0_c(c, alpha, bc, local_dt, dx);
#pragma omp parallel for schedule(dynamic)
for (int i = 0; i < n_rows; i++) {
@ -184,12 +183,17 @@ auto Diffusion::BTCSDiffusion::calc_t0_c(const DMatrixRowMajor &c,
// first, iterate over first row
for (int j = 0; j < n_cols; j++) {
y_values[0] = getBCFromFlux(bc(0, j), c(0, j), alpha(0, j));
boundary_condition tmp_bc = bc(0,j+1);
if (tmp_bc.type == Diffusion::BC_CLOSED)
continue;
y_values[0] = getBCFromFlux(tmp_bc, c(0, j), alpha(0, j));
y_values[1] = c(0, j);
y_values[2] = c(1, j);
t0_c(0, j) = time_step * alpha(0, j) *
(y_values[0] - 2 * y_values[1] + y_values[2]) / (dx * dx);
(2*y_values[0] - 3 * y_values[1] + y_values[2]) / (dx * dx);
}
// then iterate over inlet
@ -210,12 +214,17 @@ auto Diffusion::BTCSDiffusion::calc_t0_c(const DMatrixRowMajor &c,
// and finally over last row
for (int j = 0; j < n_cols; j++) {
boundary_condition tmp_bc = bc(end+1,j+1);
if (tmp_bc.type == Diffusion::BC_CLOSED)
continue;
y_values[0] = c(end - 1, j);
y_values[1] = c(end, j);
y_values[2] = getBCFromFlux(bc(end, j), c(end, j), alpha(end, j));
y_values[2] = getBCFromFlux(tmp_bc, c(end, j), alpha(end, j));
t0_c(end, j) = time_step * alpha(end, j) *
(y_values[0] - 2 * y_values[1] + y_values[2]) / (dx * dx);
(y_values[0] - 3 * y_values[1] + 2*y_values[2]) / (dx * dx);
}
return t0_c;
@ -240,7 +249,7 @@ void Diffusion::BTCSDiffusion::fillMatrixFromRow(
} else {
double sx = (alpha[0] * time_step) / (dx * dx);
A_matrix.insert(1, 1) = -1. - 3. * sx;
A_matrix.insert(1, 0) = 2.*sx;
A_matrix.insert(1, 0) = 2. * sx;
A_matrix.insert(1, 2) = sx;
}
@ -263,7 +272,7 @@ void Diffusion::BTCSDiffusion::fillMatrixFromRow(
double sx = (alpha[size - 1] * time_step) / (dx * dx);
A_matrix.insert(A_size - 2, A_size - 2) = -1. - 3. * sx;
A_matrix.insert(A_size - 2, A_size - 3) = sx;
A_matrix.insert(A_size - 2, A_size - 1) = 2.*sx;
A_matrix.insert(A_size - 2, A_size - 1) = 2. * sx;
}
A_matrix.insert(A_size - 1, A_size - 1) = 1;
@ -294,15 +303,14 @@ void Diffusion::BTCSDiffusion::fillVectorFromRow(
b_vector[j + 1] = -c[j] - t0_c_j;
}
if (!left_constant) {
// this is not correct currently.We will fix this when we are able to define
// FLUX boundary conditions
b_vector[0] = getBCFromFlux(left, c[0], alpha[0]);
}
// this is not correct currently.We will fix this when we are able to define
// FLUX boundary conditions
b_vector[0] =
(left_constant ? left.value : getBCFromFlux(left, c[0], alpha[0]));
if (!right_constant) {
b_vector[b_size - 1] = getBCFromFlux(right, c[size - 1], alpha[size - 1]);
}
b_vector[b_size - 1] =
(right_constant ? right.value
: getBCFromFlux(right, c[size - 1], alpha[size - 1]));
}
void Diffusion::BTCSDiffusion::setTimestep(double time_step) {