Test block matrix and Schur complement matrix.
static char help[] =
"...\n\n";
};
};
template <>
return MatSetValues<AssemblyTypeSelector<PETSC>>(
};
template <>
return MatSetValues<AssemblyTypeSelector<BLOCK_MAT>>(
};
template <>
return MatSetValues<AssemblyTypeSelector<BLOCK_SCHUR>>(
};
template <>
return MatSetValues<AssemblyTypeSelector<BLOCK_PRECONDITIONER_SCHUR>>(
};
constexpr
bool debug =
false;
int main(
int argc,
char *argv[]) {
try {
DMType dm_name = "DMMOFEM";
auto core_log = logging::core::get();
core_log->add_sink(
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"VOL",
"T", 0, 1);
auto shell_data =
{{simple->getDomainFEName(),
{{"V", "V"},
{"T", "T"},
{"S", "S"},
{"O", "O"},
{"V", "T"},
{"T", "V"},
{"S", "T"},
{"T", "S"},
{"T", "O"},
{"O", "T"},
{"S", "O"},
{"O", "S"}
}}}
)
);
auto [mat, block_data_ptr] = shell_data;
std::vector<std::string> fields{"T", "S", "O"};
{schur_dm, block_dm}, block_data_ptr,
fields, {nullptr, nullptr, nullptr}, true
)
);
using OpMassBlockPreconditionerAssemble =
[](int, int, int o) { return 2 * o; });
auto &pip_lhs = pip_mng->getOpDomainLhsPipeline();
pip_lhs.push_back(new OpMassPETSCAssemble("V", "V"));
pip_lhs.push_back(new OpMassPETSCAssemble("V", "T"));
pip_lhs.push_back(new OpMassPETSCAssemble("T", "V"));
pip_lhs.push_back(new OpMassPETSCAssemble("S", "S", close_zero));
pip_lhs.push_back(new OpMassPETSCAssemble("S", "T", beta));
pip_lhs.push_back(new OpMassPETSCAssemble("T", "S", beta));
pip_lhs.push_back(new OpMassPETSCAssemble("O", "O", close_zero));
pip_lhs.push_back(new OpMassPETSCAssemble("T", "O", beta));
pip_lhs.push_back(new OpMassPETSCAssemble("O", "T", beta));
pip_lhs.push_back(new OpMassPETSCAssemble("S", "O", gamma));
pip_lhs.push_back(new OpMassPETSCAssemble("O", "S", gamma));
pip_lhs.push_back(new OpMassBlockAssemble("V", "V"));
pip_lhs.push_back(new OpMassBlockAssemble("V", "T"));
pip_lhs.push_back(new OpMassBlockAssemble("T", "V"));
pip_lhs.push_back(new OpMassBlockAssemble("S", "S", close_zero));
pip_lhs.push_back(new OpMassBlockAssemble("S", "T", beta));
pip_lhs.push_back(new OpMassBlockAssemble("T", "S", beta));
pip_lhs.push_back(new OpMassBlockAssemble("O", "O", close_zero));
pip_lhs.push_back(new OpMassBlockAssemble("T", "O", beta));
pip_lhs.push_back(new OpMassBlockAssemble("O", "T", beta));
pip_lhs.push_back(new OpMassBlockAssemble("S", "O", gamma));
pip_lhs.push_back(new OpMassBlockAssemble("O", "S", gamma));
pip_lhs.push_back(new OpMassBlockPreconditionerAssemble("T", "T"));
fields,
{nullptr, nullptr, nullptr},
{nullptr, nullptr, ao_up}, {nullptr, nullptr, S},
{true, true, true}, true, block_data_ptr)
);
{
BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
MOFEM_LOG(
"Timeline", Sev::inform) <<
"Assemble start";
pip_mng->getDomainLhsFE());
MOFEM_LOG(
"Timeline", Sev::inform) <<
"Assemble end";
}
{
BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
MOFEM_LOG(
"Timeline", Sev::inform) <<
"Mat assemble start";
MOFEM_LOG(
"Timeline", Sev::inform) <<
"Mat assemble end";
}
auto get_random_vector = [&](auto dm) {
PetscRandom rctx;
PetscRandomCreate(PETSC_COMM_WORLD, &rctx);
PetscRandomDestroy(&rctx);
};
auto v = get_random_vector(
simple->getDM());
auto test = [](auto msg, auto y, double norm0) {
PetscReal norm;
CHKERR VecNorm(y, NORM_2, &norm);
norm = norm / norm0;
<< msg << ": norm of difference: " << norm;
if (norm >
eps || std::isnan(norm) || std::isinf(norm)) {
SETERRQ(PETSC_COMM_WORLD, 1, "norm of difference is too big");
}
};
std::vector<int> zero_rows_and_cols = {
0, 1, 10, 20,
500};
&*zero_rows_and_cols.begin(), 1, PETSC_NULL,
PETSC_NULL);
&*zero_rows_and_cols.begin(), 1, PETSC_NULL,
PETSC_NULL);
{
BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
<< "MatMult(petsc_mat, v, y_petsc) star";
<< "MatMult(petsc_mat, v, y_petsc) end";
}
double nrm0;
CHKERR VecNorm(y_petsc, NORM_2, &nrm0);
{
BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
<< "MatMult(block_mat, v, y_block) star";
<< "MatMult(block_mat, v, y_block) end";
}
CHKERR VecAXPY(y_petsc, -1.0, y_block);
CHKERR test(
"mult", y_petsc, nrm0);
CHKERR VecAXPY(y_petsc, -1.0, y_block);
CHKERR test(
"mult add", y_petsc, nrm0);
CHKERR MatMult(nested_mat,
v, y_nested);
CHKERR VecAXPY(y_petsc, -1.0, y_nested);
CHKERR test(
"mult nested", y_petsc, nrm0);
auto diag_mat = std::get<0>(*nested_data_ptr)[3];
auto diag_block_x = get_random_vector(block_dm);
CHKERR MatMult(diag_mat, diag_block_x, diag_block_f);
CHKERR DMSetMatType(block_dm, MATSHELL);
CHKERR DMKSPSetComputeOperators(
block_dm,
[](KSP, Mat, Mat, void *) {
MOFEM_LOG(
"WORLD", Sev::inform) <<
"empty operator";
return 0;
},
nullptr);
CHKERR KSPSetDM(ksp, block_dm);
CHKERR KSPSetFromOptions(ksp);
auto get_pc = [](auto ksp) {
PC pc_raw;
CHKERR KSPGetPC(ksp, &pc_raw);
};
CHKERR VecZeroEntries(block_solved_x);
CHKERR KSPSolve(ksp, diag_block_f, block_solved_x);
CHKERR MatMult(diag_mat, block_solved_x, diag_block_f_test);
CHKERR VecAXPY(diag_block_f_test, -1.0, diag_block_f);
CHKERR test(
"diag solve", diag_block_f_test, nrm0);
}
}
return 0;
}