Skip to content

Commit

Permalink
add solid body rotation in 3d using coarse grid
Browse files Browse the repository at this point in the history
  • Loading branch information
@JuntaoHuang
JuntaoHuang committed Apr 1, 2024
1 parent bd3b827 commit 28606f7
Show file tree
Hide file tree
Showing 3 changed files with 203 additions and 31 deletions.
111 changes: 80 additions & 31 deletions example/02_hyperbolic_02_var_coeff_coarse_grid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@
int main(int argc, char *argv[])
{
// constant variable
const int DIM = 2;
const int DIM = 3;

// static variables
AlptBasis::PMAX = 1;
Expand Down Expand Up @@ -85,6 +85,12 @@ int main(int argc, char *argv[])
// filter coefficient
double filter_coef = 1.0;

// numerical test number
// 1: example 4.2 (solid body rotation in 2D) in Guo and Cheng, SISC (2016)
// 2: example 4.2 (solid body rotation in 3D) in Guo and Cheng, SISC (2016)
// 3: example 4.3 (deformational flow in 2D) in Guo and Cheng, SISC (2016)
int test_num = 2;

OptionsParser args(argc, argv);
args.AddOption(&NMAX, "-NM", "--max-mesh-level", "Maximum mesh level");
args.AddOption(&N_init, "-N0", "--initial-mesh-level", "Mesh level in initialization");
Expand All @@ -98,6 +104,7 @@ int main(int argc, char *argv[])
args.AddOption(&NMAX_coarse_grid_stage_1, "-N1", "--max-mesh-stage-1", "Maximum mesh level in stage 1");
args.AddOption(&NMAX_coarse_grid_stage_2, "-N2", "--max-mesh-stage-2", "Maximum mesh level in stage 2");
args.AddOption(&filter_coef, "-filter", "--filter-coefficient", "Filter coefficient");
args.AddOption(&test_num, "-test", "--test-number", "Numerical test number");

args.Parse();
if (!args.Good())
Expand All @@ -110,6 +117,11 @@ int main(int argc, char *argv[])
// check mesh level in initialization should be less or equal to maximum mesh level
if (N_init>NMAX) { std::cout << "Mesh level in initialization should not be larger than Maximum mesh level" << std::endl; return 1; }
bool sparse = ((is_init_sparse==1) ? true : false);

// check if the test number is valid
if (test_num < 1 || test_num > 3) { std::cout << "Invalid test number" << std::endl; return 1; }
if (((test_num == 1) || (test_num == 3)) && DIM != 2) { std::cout << "Test number 1 and 3 is only valid for 2D" << std::endl; return 1; }
if (test_num == 2 && DIM != 3) { std::cout << "Test number 2 is only valid for 3D" << std::endl; return 1; }

// initialize hash key
Hash hash;
Expand Down Expand Up @@ -141,14 +153,28 @@ int main(int argc, char *argv[])
// adaptive interpolation for initial function
auto init_func = [=](std::vector<double> x, int i)
{
// example 4.2 (solid body rotation) in Guo and Cheng. "A sparse grid discontinuous Galerkin method for high-dimensional transport equations and its application to kinetic simulations." SISC (2016)
const std::vector<double> xc{0.75, 0.5};
const double b = 0.23;

// // example 4.3 (deformational flow)
// const std::vector<double> xc{0.65, 0.5};
// const double b = 0.35;

std::vector<double> xc;
double b = 0.;

// example 4.2 (deformational flow in 2D)
if (test_num == 1)
{
xc = {0.75, 0.5};
b = 0.23;
}
// example 4.2 (deformational flow in 3D)
else if (test_num == 2)
{
xc = {0.5, 0.55, 0.5};
b = 0.45;
}
// example 4.3 (deformational flow in 2D)
else if (test_num == 3)
{
xc = {0.65, 0.5};
b = 0.35;
}

double r_sqr = 0.;
for (int d = 0; d < DIM; d++) { r_sqr += pow(x[d] - xc[d], 2.); };
double r = pow(r_sqr, 0.5);
Expand All @@ -168,16 +194,23 @@ int main(int argc, char *argv[])
IO inout(dg_f);

// ------------------------------
HyperbolicDiffFluxLagrRHS fast_rhs_lagr(dg_f, oper_matx_lagr);
HyperbolicLagrRHS fast_rhs_lagr(dg_f, oper_matx_lagr);

// fast Lagrange interpolation
LagrInterpolation interp_lagr(dg_f);
FastLagrIntp fast_lagr_intp(dg_f, interp_lagr.Lag_pt_Alpt_1D, interp_lagr.Lag_pt_Alpt_1D_d1);

// constant in global Lax-Friedrich flux
const double lxf_alpha = 0.5;
// wave speed in x and y direction
const std::vector<double> wave_speed{0.5, 0.5};
// wave speed
std::vector<double> wave_speed;
if (test_num == 1)
{
wave_speed = {0.5, 0.5};
}
else if (test_num == 2)
{
wave_speed = {sqrt(2.)/4., sqrt(2.)/2., sqrt(2.)/4.};
}
const std::vector<double> lxf_alpha = wave_speed;

const int max_mesh = dg_f.max_mesh_level();
const double dx = 1./pow(2., max_mesh);
Expand All @@ -196,20 +229,24 @@ int main(int argc, char *argv[])
// --- Part 3: time evolution ---
// // linear operator
HyperbolicAlpt linear_stage_1(dg_f, oper_matx_alpt);
// x direction: u^- * [v] * alpha / 2
linear_stage_1.assemble_matrix_flx_scalar_coarse_grid(0, -1, NMAX_coarse_grid_stage_1, lxf_alpha/2);
// x direction: - u^+ * [v] * alpha / 2
linear_stage_1.assemble_matrix_flx_scalar_coarse_grid(0, 1, NMAX_coarse_grid_stage_1, -lxf_alpha/2);
// y direction: u^- * [v] * alpha / 2
linear_stage_1.assemble_matrix_flx_scalar_coarse_grid(1, -1, NMAX_coarse_grid_stage_1, lxf_alpha/2);
// y direction: - u^+ * [v] * alpha / 2
linear_stage_1.assemble_matrix_flx_scalar_coarse_grid(1, 1, NMAX_coarse_grid_stage_1, -lxf_alpha/2);
for (int d = 0; d < DIM; d++)
{
// u^- * [v] * alpha / 2
linear_stage_1.assemble_matrix_flx_scalar_coarse_grid(d, -1, NMAX_coarse_grid_stage_1, lxf_alpha[d]/2);

// - u^+ * [v] * alpha / 2
linear_stage_1.assemble_matrix_flx_scalar_coarse_grid(d, 1, NMAX_coarse_grid_stage_1, -lxf_alpha[d]/2);
}

HyperbolicAlpt linear_stage_2(dg_f, oper_matx_alpt);
linear_stage_2.assemble_matrix_flx_scalar_coarse_grid(0, -1, NMAX_coarse_grid_stage_2, lxf_alpha/2);
linear_stage_2.assemble_matrix_flx_scalar_coarse_grid(0, 1, NMAX_coarse_grid_stage_2, -lxf_alpha/2);
linear_stage_2.assemble_matrix_flx_scalar_coarse_grid(1, -1, NMAX_coarse_grid_stage_2, lxf_alpha/2);
linear_stage_2.assemble_matrix_flx_scalar_coarse_grid(1, 1, NMAX_coarse_grid_stage_2, -lxf_alpha/2);
for (int d = 0; d < DIM; d++)
{
// u^- * [v] * alpha / 2
linear_stage_2.assemble_matrix_flx_scalar_coarse_grid(d, -1, NMAX_coarse_grid_stage_2, lxf_alpha[d]/2);

// - u^+ * [v] * alpha / 2
linear_stage_2.assemble_matrix_flx_scalar_coarse_grid(d, 1, NMAX_coarse_grid_stage_2, -lxf_alpha[d]/2);
}

std::cout << "--- evolution started ---" << std::endl;
Timer record_time;
Expand All @@ -230,13 +267,25 @@ int main(int argc, char *argv[])
odeSolver.init();

for ( int stage = 0; stage < odeSolver.num_stage; ++stage )
{
// solid body rotation: f_t + (-x2 + 0.5) * f_x1 + (x1 - 0.5) * f_x2 = 0
{
auto coe_func = [&](std::vector<double> x, int d) -> double
{
if (d==0) { return -x[1] + 0.5; }
else { return x[0] - 0.5; }
// return 1.0;
// a = (-x2 + 0.5, x1 - 0.5)
if (test_num == 1)
{
if (d==0) { return -x[1] + 0.5; }
else { return x[0] - 0.5; }
}

// a = (-sqrt(2)/2 * (x2 - 0.5),
// sqrt(2)/2 * (x1 - 0.5) + sqrt(2)/2 * (x3 - 0.5)
// -sqrt(2)/2 * (x2 - 0.5))
else if (test_num == 2)
{
if (d==0) { return -sqrt(2.)/2. * (x[1] - 0.5); }
else if (d==1) { return sqrt(2.)/2. * (x[0] - 0.5) + sqrt(2.)/2. * (x[2] - 0.5); }
else { return -sqrt(2.)/2. * (x[1] - 0.5); }
}
};

if (stage == 0)
Expand Down
13 changes: 13 additions & 0 deletions include/FastMultiplyLU.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -391,6 +391,8 @@ class FastRHS:

void transform_fucoe_to_rhs(const std::vector<const VecMultiD<double>*> & mat_1D, const std::vector<std::string> & operator_type, const int dim_interp, const double coefficient = 1.0, const int vec_index = 0);

void transform_fucoe_to_rhs_coarse_grid(const std::vector<const VecMultiD<double>*> & mat_1D, const std::vector<std::string> & operator_type, const int dim_interp, const int mesh_nmax, const double coefficient = 1.0, const int vec_index = 0);

// transform Element::ucoe_alpt to Element::rhs
// use this function to compute rhs for linear equations
void transform_ucoe_alpt_to_rhs(const std::vector<const VecMultiD<double>*> & mat_1D, const std::vector<std::string> & operator_type, const int dim_operator, const double coefficient = 1.0, const int vec_index = 0);
Expand Down Expand Up @@ -420,9 +422,16 @@ class FastRHS:
const std::vector<int> & size_trans_from, const std::vector<int> & size_trans_to, const int dim_transform_1D,
const bool is_first_step, const bool is_last_step, const int dim_interp, const int vec_index, const double coefficient = 1.0);

void transform_1D_from_fucoe_to_rhs_coarse_grid(const VecMultiD<double> & mat_1D, const std::string & LU, const std::string & operator_type,
const std::vector<int> & size_trans_from, const std::vector<int> & size_trans_to, const int dim_transform_1D,
const bool is_first_step, const bool is_last_step, const int dim_interp, const int vec_index, const int mesh_nmax, const double coefficient = 1.0);

void transform_multiD_partial_sum(const std::vector<const VecMultiD<double>*> & mat_1D_array, const std::vector<int> & dim_order_transform,
const std::vector<std::string> & LU_order_transform, const std::vector<std::string> & operator_type, const int dim_interp, const int vec_index, const double coefficient = 1.0);

void transform_multiD_partial_sum_coarse_grid(const std::vector<const VecMultiD<double>*> & mat_1D_array, const std::vector<int> & dim_order_transform,
const std::vector<std::string> & LU_order_transform, const std::vector<std::string> & operator_type, const int dim_interp, const int vec_index, const int mesh_nmax, const double coefficient = 1.0);

void transform_1D_from_ucoe_alpt_to_rhs(const VecMultiD<double> & mat_1D, const std::string & LU, const std::string & operator_type,
const std::vector<int> & size_trans_from, const std::vector<int> & size_trans_to, const int dim_transform_1D,
const bool is_first_step, const bool is_last_step, const int dim_interp, const int vec_index, const double coefficient = 1.0);
Expand Down Expand Up @@ -570,9 +579,13 @@ class HyperbolicLagrRHS:
// volume integral
void rhs_vol_scalar();

void rhs_vol_scalar_coarse_grid(const int mesh_nmax);

// flux integral on interpolation part
void rhs_flx_intp_scalar();

void rhs_flx_intp_scalar_coarse_grid(const int mesh_nmax);

private:
OperatorMatrix1D<LagrBasis, AlptBasis>* const oper_matx_lagr_ptr;
};
Expand Down
110 changes: 110 additions & 0 deletions source/FastMultiplyLU.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,20 @@ void FastRHS::transform_fucoe_to_rhs(const std::vector<const VecMultiD<double>*>
}
}

void FastRHS::transform_fucoe_to_rhs_coarse_grid(const std::vector<const VecMultiD<double>*> & mat_1D, const std::vector<std::string> & operator_type, const int dim_interp, const int mesh_nmax, const double coefficient, const int vec_index)
{
const int dim = dgsolution_ptr->DIM;

std::vector<std::vector<int>> dim_order_transform;
std::vector<std::vector<std::string>> LU_order_transform;
generate_transform_order(dim, dim_order_transform, LU_order_transform);

for (size_t sum_num = 0; sum_num < dim_order_transform.size(); sum_num++)
{
transform_multiD_partial_sum_coarse_grid(mat_1D, dim_order_transform[sum_num], LU_order_transform[sum_num], operator_type, dim_interp, vec_index, mesh_nmax, coefficient);
}
}

void FastRHS::transform_ucoe_alpt_to_rhs(const std::vector<const VecMultiD<double>*> & mat_1D, const std::vector<std::string> & operator_type, const int dim_interp, const double coefficient, const int vec_index)
{
const int dim = dgsolution_ptr->DIM;
Expand Down Expand Up @@ -78,6 +92,32 @@ void FastRHS::transform_1D_from_fucoe_to_rhs(const VecMultiD<double> & mat_1D, c
if (is_last_step) { add_transto_to_rhs(vec_index); }
}

void FastRHS::transform_1D_from_fucoe_to_rhs_coarse_grid(const VecMultiD<double> & mat_1D, const std::string & LU, const std::string & operator_type,
const std::vector<int> & size_trans_from, const std::vector<int> & size_trans_to, const int dim_transform_1D,
const bool is_first_step, const bool is_last_step, const int dim_interp, const int vec_index, const int mesh_nmax, const double coefficient)
{
// transform from interpolation basis to alpert basis
const int pmax_trans_from = Element::PMAX_intp;
const int pmax_trans_to = Element::PMAX_alpt;

resize_ucoe_transfrom(size_trans_from, vec_index);

// if first step, then copy value in ucoe_alpt to ucoe_trans_from
// else, copy value ucoe_trans_to to ucoe_trans_from
if (is_first_step) { copy_fucoe_to_transfrom(dim_interp, vec_index); }
else { copy_transto_to_transfrom(vec_index); }

resize_ucoe_transto(size_trans_to, vec_index);

// do transformation in 1D
// only multiply coefficient in the first step
if (is_first_step) { transform_1D_coarse_grid(mat_1D, LU, operator_type, dim_transform_1D, pmax_trans_from, pmax_trans_to, mesh_nmax, coefficient, vec_index, vec_index); }
else { transform_1D_coarse_grid(mat_1D, LU, operator_type, dim_transform_1D, pmax_trans_from, pmax_trans_to, mesh_nmax, 1.0, vec_index, vec_index); }

// add transform_to to rhs in the last step
if (is_last_step) { add_transto_to_rhs(vec_index); }
}

void FastRHS::transform_multiD_partial_sum(const std::vector<const VecMultiD<double>*> & mat_1D_array, const std::vector<int> & dim_order_transform,
const std::vector<std::string> & LU_order_transform, const std::vector<std::string> & operator_type, const int dim_interp, const int vec_index, const double coefficient)
{
Expand All @@ -101,6 +141,28 @@ void FastRHS::transform_multiD_partial_sum(const std::vector<const VecMultiD<dou
}
}

void FastRHS::transform_multiD_partial_sum_coarse_grid(const std::vector<const VecMultiD<double>*> & mat_1D_array, const std::vector<int> & dim_order_transform, const std::vector<std::string> & LU_order_transform, const std::vector<std::string> & operator_type, const int dim_interp, const int vec_index, const int mesh_nmax, const double coefficient)
{
// transform from alpert basis to interpolation basis
const int pmax_trans_from = Element::PMAX_intp;
const int pmax_trans_to = Element::PMAX_alpt;
const int dim = Element::DIM;

std::vector<std::vector<int>> series_vec = series_vec_transform_size(pmax_trans_from+1, pmax_trans_to+1, dim, dim_order_transform);

// variable control first or last step
std::vector<bool> is_first_step(dim, false);
is_first_step[0] = true;

std::vector<bool> is_last_step(dim, false);
is_last_step[dim-1] = true;

for (size_t d = 0; d < dim; d++)
{
transform_1D_from_fucoe_to_rhs_coarse_grid(*(mat_1D_array[dim_order_transform[d]]), LU_order_transform[d], operator_type[dim_order_transform[d]], series_vec[d], series_vec[d+1], dim_order_transform[d], is_first_step[d], is_last_step[d], dim_interp, vec_index, mesh_nmax, coefficient);
}
}

void FastRHS::transform_multiD_partial_sum_ucoe_alpt(const std::vector<const VecMultiD<double>*> & mat_1D_array, const std::vector<int> & dim_order_transform,
const std::vector<std::string> & LU_order_transform, const std::vector<std::string> & operator_type, const int dim_interp, const int vec_index, const double coefficient)
{
Expand Down Expand Up @@ -1078,6 +1140,24 @@ void HyperbolicLagrRHS::rhs_vol_scalar()
}
}

void HyperbolicLagrRHS::rhs_vol_scalar_coarse_grid(const int mesh_nmax)
{
const int dim = dgsolution_ptr->DIM;

std::vector<const VecMultiD<double>*> oper_matx_1D;
std::vector<std::string> operator_type(dim, "vol");
for (size_t d = 0; d < dim; d++)
{
oper_matx_1D.clear();
for (size_t dim_derivative = 0; dim_derivative < dim; dim_derivative++)
{
if (dim_derivative == d) { oper_matx_1D.push_back(&(oper_matx_lagr_ptr->u_vx)); }
else { oper_matx_1D.push_back(&(oper_matx_lagr_ptr->u_v)); }
}
transform_fucoe_to_rhs_coarse_grid(oper_matx_1D, operator_type, d, mesh_nmax);
}
}

void HyperbolicLagrRHS::rhs_flx_intp_scalar()
{
const int dim = dgsolution_ptr->DIM;
Expand Down Expand Up @@ -1108,6 +1188,36 @@ void HyperbolicLagrRHS::rhs_flx_intp_scalar()
}
}

void HyperbolicLagrRHS::rhs_flx_intp_scalar_coarse_grid(const int mesh_nmax)
{
const int dim = dgsolution_ptr->DIM;

const VecMultiD<double> oper_matx_uave_vjp = oper_matx_lagr_ptr->ulft_vjp + oper_matx_lagr_ptr->urgt_vjp;

std::vector<const VecMultiD<double>*> oper_matx_1D;
std::vector<std::string> operator_type;
for (size_t d = 0; d < dim; d++)
{
oper_matx_1D.clear();
operator_type.clear();

for (size_t dim_derivative = 0; dim_derivative < dim; dim_derivative++)
{
if (dim_derivative == d)
{
oper_matx_1D.push_back(&oper_matx_uave_vjp);
operator_type.push_back("flx");
}
else
{
oper_matx_1D.push_back(&oper_matx_lagr_ptr->u_v);
operator_type.push_back("vol");
}
}
transform_fucoe_to_rhs_coarse_grid(oper_matx_1D, operator_type, d, mesh_nmax, 0.5);
}
}

void HyperbolicHermRHS::rhs_vol_scalar()
{
const int dim = dgsolution_ptr->DIM;
Expand Down

0 comments on commit 28606f7

Please sign in to comment.