#include <boost/math/quadrature/gauss_kronrod.hpp>
using namespace boost::math::quadrature;
};
PETSC>::LinearForm<GAUSS>::OpBaseTimesScalar<1>;
constexpr
double omega = 7.292 * 1e-5;
constexpr
double g = 9.80616;
constexpr
double mu = 1e4;
constexpr
double h0 = 1e4;
constexpr
double h_hat = 120;
constexpr
double u_max = 80;
constexpr
double phi_0 = M_PI / 7;
constexpr
double phi_2 = M_PI / 4;
boost::shared_ptr<MatrixDouble> u_dot_ptr,
boost::shared_ptr<MatrixDouble> grad_u_ptr,
boost::shared_ptr<MatrixDouble> grad_h_ptr)
uPtr(u_ptr), uDotPtr(u_dot_ptr), uGradPtr(grad_u_ptr),
hGradPtr(grad_h_ptr) {}
const double vol = getMeasure();
auto t_w = getFTensor0IntegrationWeight();
auto t_dot_u = getFTensor1FromMat<3>(*uDotPtr);
auto t_u = getFTensor1FromMat<3>(*uPtr);
auto t_grad_u = getFTensor2FromMat<3, 3>(*uGradPtr);
auto t_grad_h = getFTensor1FromMat<3>(*hGradPtr);
auto t_coords = getFTensor1CoordsAtGaussPts();
auto t_normal = getFTensor1NormalsAtGaussPts();
for (int gg = 0; gg != nbIntegrationPts; gg++) {
const double alpha = t_w * vol;
auto t_nf = getFTensor1FromArray<3, 3>(locF);
const auto a = std::sqrt(t_coords(
i) * t_coords(
i));
const auto sin_fi = t_coords(2) /
a;
const auto f = 2 *
omega * sin_fi;
t_P(
i,
j) = t_r(
i) * t_r(
j);
t_rhs(
m) = t_Q(
m,
i) * (t_dot_u(
i) + t_grad_u(
i,
j) * t_u(
j) +
f * t_A(
i,
j) * t_u(
j) +
g * t_grad_h(
i));
t_rhs_grad(
m,
j) = t_Q(
m,
i) * (
mu * t_grad_u(
i,
j));
t_rhs(
m) += t_P(
m,
j) * t_u(
j);
int rr = 0;
for (; rr != nbRows / 3; ++rr) {
t_nf(
i) += alpha * t_row_base * t_rhs(
i);
t_nf(
i) += alpha * t_row_diff_base(
j) * t_rhs_grad(
i,
j);
++t_row_base;
++t_row_diff_base;
++t_nf;
}
for (; rr < nbRowBaseFunctions; ++rr) {
++t_row_base;
++t_row_diff_base;
}
++t_w;
++t_u;
++t_dot_u;
++t_grad_u;
++t_grad_h;
++t_coords;
++t_normal;
}
}
private:
boost::shared_ptr<MatrixDouble> uPtr;
boost::shared_ptr<MatrixDouble> uDotPtr;
boost::shared_ptr<MatrixDouble> uGradPtr;
boost::shared_ptr<MatrixDouble> hGradPtr;
};
OpULhs_dU(
const std::string field_name_row,
const std::string field_name_col,
boost::shared_ptr<MatrixDouble> u_ptr,
boost::shared_ptr<MatrixDouble> grad_u_ptr)
uPtr(u_ptr), uGradPtr(grad_u_ptr) {
this->sYmm = false;
}
const double vol = getMeasure();
auto t_w = getFTensor0IntegrationWeight();
auto t_coords = getFTensor1CoordsAtGaussPts();
auto t_normal = getFTensor1NormalsAtGaussPts();
auto t_u = getFTensor1FromMat<3>(*uPtr);
auto t_grad_u = getFTensor2FromMat<3, 3>(*uGradPtr);
auto get_t_mat = [&](const int rr) {
&locMat(rr + 0, 0), &locMat(rr + 0, 1), &locMat(rr + 0, 2),
&locMat(rr + 1, 0), &locMat(rr + 1, 1), &locMat(rr + 1, 2),
&locMat(rr + 2, 0), &locMat(rr + 2, 1), &locMat(rr + 2, 2)};
};
const auto ts_a = getFEMethod()->ts_a;
for (int gg = 0; gg != nbIntegrationPts; gg++) {
const auto a = std::sqrt(t_coords(
i) * t_coords(
i));
const auto sin_fi = t_coords(2) /
a;
const auto f = 2 *
omega * sin_fi;
t_P(
i,
j) = t_r(
i) * t_r(
j);
t_Q(
m,
i) * (ts_a *
t_kd(
i,
j) + t_grad_u(
i,
j) +
f * t_A(
i,
j)) +
const double alpha = t_w * vol;
int rr = 0;
for (; rr != nbRows / 3; rr++) {
auto t_mat = get_t_mat(3 * rr);
for (int cc = 0; cc != nbCols / 3; cc++) {
t_mat(
i,
j) += (alpha * t_row_base * t_col_base) * t_rhs_du(
i,
j);
t_mat(
i,
j) += (alpha *
mu) * t_Q(
i,
j) *
(t_row_diff_base(
m) * t_col_diff_base(
m));
++t_col_diff_base;
++t_col_base;
++t_mat;
}
++t_row_base;
++t_row_diff_base;
}
for (; rr < nbRowBaseFunctions; ++rr) {
++t_row_base;
++t_row_diff_base;
}
++t_w;
++t_coords;
++t_normal;
++t_u;
++t_grad_u;
}
}
private:
boost::shared_ptr<MatrixDouble> uPtr;
boost::shared_ptr<MatrixDouble> uGradPtr;
};
OpULhs_dH(
const std::string field_name_row,
const std::string field_name_col)
this->sYmm = false;
}
const double vol = getMeasure();
auto t_w = getFTensor0IntegrationWeight();
auto t_normal = getFTensor1NormalsAtGaussPts();
auto get_t_vec = [&](const int rr) {
&locMat(rr + 0, 0),
&locMat(rr + 1, 0),
&locMat(rr + 2, 0)};
};
for (int gg = 0; gg != nbIntegrationPts; gg++) {
t_P(
i,
j) = t_r(
i) * t_r(
j);
const double alpha = t_w * vol;
int rr = 0;
for (; rr != nbRows / 3; rr++) {
auto t_vec = get_t_vec(3 * rr);
const double a = alpha *
g * t_row_base;
for (int cc = 0; cc != nbCols; cc++) {
t_vec(
i) +=
a * (t_Q(
i,
m) * t_col_diff_base(
m));
++t_vec;
++t_col_diff_base;
}
++t_row_base;
}
for (; rr < nbRowBaseFunctions; ++rr)
++t_row_base;
++t_w;
++t_normal;
}
}
};
private:
boost::shared_ptr<FEMethod> domianLhsFEPtr;
boost::shared_ptr<FEMethod> domianRhsFEPtr;
};
}
}
}
}
PetscBool is_restart = PETSC_FALSE;
PETSC_NULL);
auto restart_vector = [&]() {
<< "reading vector in binary from vector.dat ...";
PetscViewer viewer;
PetscViewerBinaryOpen(PETSC_COMM_WORLD, "vector.dat", FILE_MODE_READ,
&viewer);
VecLoad(T, viewer);
};
if (is_restart) {
} else {
const double e_n = exp(-4 / pow(
phi_1 -
phi_0, 2));
auto get_phi = [&](const double x, const double y, const double z) {
const double r = std::sqrt(t_r(
i) * t_r(
i));
};
auto init_u_phi = [&](
const double phi) {
} else {
return 0.;
}
};
auto init_u = [&](
const double x,
const double y,
const double z) {
const double u_phi = init_u_phi(get_phi(x, y, z));
if (u_phi > 0) {
t_u(
i) = ((t_A(
i,
j) * t_r(
j)) * u_phi);
}
return t_u;
};
auto init_h = [&](
const double x,
const double y,
const double z) {
const double a = std::sqrt(x * x + y * y + z * z);
auto integral = [&](const double fi) {
const double u_phi = init_u_phi(fi);
const auto f = 2 *
omega * sin(fi);
return a * u_phi * (
f + (tan(fi) /
a) * u_phi);
};
auto montain = [&](
const double lambda,
const double fi) {
else
return 0.;
};
const double fi = get_phi(x, y, z);
const double lambda = atan2(x, y);
double h1 = 0;
h1 = gauss_kronrod<double, 32>::integrate(
integral,
phi_0, fi, 0, std::numeric_limits<float>::epsilon());
return h0 + (montain(
lambda, fi) - (h1 /
g));
};
auto set_domain_general = [&](auto &pipeline) {
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
};
auto set_domain_rhs = [&](auto &pipeline) {
};
auto set_domain_lhs = [&](auto &pipeline) {
pipeline.push_back(
new OpMassUU(
"U",
"U", [](
double,
double,
double) {
return 1; }));
pipeline.push_back(
new OpMassHH(
"H",
"H", [](
double,
double,
double) {
return 1; }));
};
auto post_proc = [&]() {
auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
auto u_ptr = boost::make_shared<MatrixDouble>();
auto h_ptr = boost::make_shared<VectorDouble>();
auto pos_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(post_proc_fe->getPostProcMesh(),
post_proc_fe->getMapGaussPts(),
{{"H", h_ptr}},
{{"U", u_ptr}, {"HO_POSITIONS", pos_ptr}},
{}, {}
)
);
CHKERR post_proc_fe->writeFile(
"out_init.h5m");
};
auto solve_init = [&]() {
auto solver = pipeline_mng->createKSP();
CHKERR KSPSetFromOptions(solver);
PC pc;
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
if (is_pcfs == PETSC_TRUE) {
auto bc_mng = mField.getInterface<
BcManager>();
auto name_prb =
simple->getProblemName();
name_prb,
ROW,
"U", 0, 3, is_u);
CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_u);
CHKERR PCFieldSplitSetType(pc, PC_COMPOSITE_ADDITIVE);
}
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
};
};
set_domain_general(pipeline_mng->getOpDomainRhsPipeline());
set_domain_rhs(pipeline_mng->getOpDomainRhsPipeline());
set_domain_general(pipeline_mng->getOpDomainLhsPipeline());
set_domain_lhs(pipeline_mng->getOpDomainLhsPipeline());
pipeline_mng->getOpDomainRhsPipeline().clear();
pipeline_mng->getOpDomainLhsPipeline().clear();
}
}
auto set_domain_general = [&](auto &pipeline) {
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
};
auto set_domain_rhs = [&](auto &pipeline) {
auto dot_u_ptr = boost::make_shared<MatrixDouble>();
auto u_ptr = boost::make_shared<MatrixDouble>();
auto grad_u_ptr = boost::make_shared<MatrixDouble>();
auto div_u_ptr = boost::make_shared<VectorDouble>();
auto dot_h_ptr = boost::make_shared<VectorDouble>();
auto grad_h_ptr = boost::make_shared<MatrixDouble>();
pipeline.push_back(
pipeline.push_back(
"H", dot_h_ptr, [](double, double, double) { return 1.; }));
"H", div_u_ptr, [](
double,
double,
double) {
return h0; }));
pipeline.push_back(
new OpURhs(
"U", u_ptr, dot_u_ptr, grad_u_ptr, grad_h_ptr));
};
auto set_domain_lhs = [&](auto &pipeline) {
auto u_ptr = boost::make_shared<MatrixDouble>();
auto grad_u_ptr = boost::make_shared<MatrixDouble>();
auto grad_h_ptr = boost::make_shared<MatrixDouble>();
pipeline.push_back(
pipeline.push_back(
new OpMassHH(
"H",
"H", [&](
double,
double,
double) {
return domianLhsFEPtr->ts_a;
}));
"H",
"U", [](
const double,
const double,
const double) {
return h0; },
false, false));
pipeline.push_back(
pipeline.push_back(
pipeline.push_back(
new OpULhs_dU(
"U",
"U", u_ptr, grad_u_ptr));
};
};
set_domain_general(pipeline_mng->getOpDomainRhsPipeline());
set_domain_general(pipeline_mng->getOpDomainLhsPipeline());
set_domain_rhs(pipeline_mng->getOpDomainRhsPipeline());
set_domain_lhs(pipeline_mng->getOpDomainLhsPipeline());
domianLhsFEPtr = pipeline_mng->getDomainLhsFE();
domianRhsFEPtr = pipeline_mng->getDomainRhsFE();
}
: dM(dm), postProc(post_proc){};
"out_step_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
<< "writing vector in binary to vector.dat ...";
PetscViewer viewer;
PetscViewerBinaryOpen(PETSC_COMM_WORLD, "vector.dat", FILE_MODE_WRITE,
&viewer);
VecView(ts_u, viewer);
PetscViewerDestroy(&viewer);
}
}
private:
boost::shared_ptr<PostProcEle> postProc;
};
auto set_initial_step = [&](auto ts) {
int step = 0;
PETSC_NULL);
CHKERR TSSetStepNumber(ts, step);
};
auto set_fieldsplit_preconditioner_ksp = [&](auto ksp) {
PC pc;
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
if (is_pcfs == PETSC_TRUE) {
auto bc_mng = mField.getInterface<
BcManager>();
auto name_prb =
simple->getProblemName();
name_prb,
ROW,
"U", 0, 3, is_u);
CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_u);
}
};
auto get_fe_post_proc = [&]() {
auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
auto u_ptr = boost::make_shared<MatrixDouble>();
auto h_ptr = boost::make_shared<VectorDouble>();
auto pos_ptr = boost::make_shared<MatrixDouble>();
auto grad_u_ptr = boost::make_shared<MatrixDouble>();
auto grad_h_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{{"H", h_ptr}},
{{"U", u_ptr}, {"HO_POSITIONS", pos_ptr}, {"GRAD_H", grad_h_ptr}},
{{"GRAD_U", grad_u_ptr}}, {}
)
);
return post_proc_fe;
};
auto set_fieldsplit_preconditioner_ts = [&](auto solver) {
SNES snes;
CHKERR TSGetSNES(solver, &snes);
KSP ksp;
CHKERR SNESGetKSP(snes, &ksp);
CHKERR set_fieldsplit_preconditioner_ksp(ksp);
};
ts = pipeline_mng->createTSIM();
boost::shared_ptr<FEMethod> null_fe;
auto monitor_ptr = boost::make_shared<Monitor>(dm, get_fe_post_proc());
null_fe, monitor_ptr);
double ftime = 1;
CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
SCATTER_FORWARD);
CHKERR set_fieldsplit_preconditioner_ts(ts);
}
}
}
static char help[] =
"...\n\n";
int main(
int argc,
char *argv[]) {
const char param_file[] = "param_file.petsc";
auto core_log = logging::core::get();
core_log->add_sink(LogManager::createSink(LogManager::getStrmWorld(), "SW"));
LogManager::setLog("SW");
try {
DMType dm_name = "DMMOFEM";
}
}