228 const string default_options =
"-ksp_type fgmres \n"
230 "-pc_factor_mat_solver_type mumps \n"
231 "-mat_mumps_icntl_20 0 \n"
235 "-snes_type newtonls \n"
236 "-snes_linesearch_type basic \n"
237 "-snes_max_it 100 \n"
243 string param_file =
"param_file.petsc";
244 if (!
static_cast<bool>(ifstream(param_file))) {
245 std::ofstream file(param_file.c_str(), std::ios::ate);
246 if (file.is_open()) {
247 file << default_options;
255 auto core_log = logging::core::get();
257 LogManager::createSink(LogManager::getStrmWorld(),
"DYNAMIC"));
258 LogManager::setLog(
"DYNAMIC");
266 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
267 auto moab_comm_wrap =
268 boost::make_shared<WrapMPIComm>(PETSC_COMM_WORLD,
false);
270 pcomm =
new ParallelComm(&moab, moab_comm_wrap->get_comm());
274 PetscBool is_partitioned = PETSC_FALSE;
275 PetscBool linear = PETSC_TRUE;
276 PetscInt disp_order = 1;
277 PetscInt vel_order = 1;
278 PetscBool is_solve_at_time_zero = PETSC_FALSE;
280 auto read_command_line_parameters = [&]() {
282 PetscBool flg = PETSC_TRUE;
285 if (flg != PETSC_TRUE)
286 SETERRQ(PETSC_COMM_SELF, 1,
"Error -my_file (mesh file needed)");
291 &is_partitioned, &flg);
296 enum bases { LEGENDRE, LOBATTO, BERNSTEIN_BEZIER, LASBASETOP };
297 const char *list_bases[] = {
"legendre",
"lobatto",
"bernstein_bezier"};
298 PetscInt choice_base_value = BERNSTEIN_BEZIER;
300 LASBASETOP, &choice_base_value, PETSC_NULL);
301 if (choice_base_value == LEGENDRE)
303 else if (choice_base_value == LOBATTO)
305 else if (choice_base_value == BERNSTEIN_BEZIER)
310 if (flg != PETSC_TRUE)
315 if (flg != PETSC_TRUE)
316 vel_order = disp_order;
319 "-my_solve_at_time_zero",
320 &is_solve_at_time_zero, &flg);
325 auto read_mesh = [&]() {
327 if (is_partitioned == PETSC_TRUE) {
330 option =
"PARALLEL=BCAST_DELETE;"
331 "PARALLEL_RESOLVE_SHARED_ENTS;"
332 "PARTITION=PARALLEL_PARTITION;";
342 CHKERR read_command_line_parameters();
352 CHKERR moab.create_meshset(MESHSET_SET, meshset_level0);
367 bool check_if_spatial_field_exist = m_field.
check_field(
"DISPLACEMENT");
388 fluid_pressure_fe.addNeumannFluidPressureBCElements(
"DISPLACEMENT");
391 fluid_pressure_fe.setNeumannFluidPressureFiniteElementOperators(
392 "DISPLACEMENT", PETSC_NULL,
false,
true);
432 CHKERR elastic_materials.setBlocks(elastic.setOfBlocks);
438 CHKERR elastic.addElement(
"ELASTIC",
"DISPLACEMENT");
444 false,
false,
false);
445 CHKERR elastic.setOperators(
"DISPLACEMENT",
"MESH_NODE_POSITIONS",
false,
451 CHKERR elastic_materials.setBlocks(inertia.setOfBlocks);
452 CHKERR inertia.addConvectiveMassElement(
"MASS_ELEMENT",
"VELOCITY",
454 CHKERR inertia.addVelocityElement(
"VELOCITY_ELEMENT",
"VELOCITY",
459 string name =
"-my_accelerogram";
460 char time_file_name[255];
463 time_file_name, 255, &flg);
464 if (flg == PETSC_TRUE) {
471 CHKERR elastic_materials.setBlocks(damper.blockMaterialDataMap);
486 "DAMPER",
"MESH_NODE_POSITIONS");
488 std::map<int, KelvinVoigtDamper::BlockMaterialData>::iterator
bit =
489 damper.blockMaterialDataMap.begin();
490 for (;
bit != damper.blockMaterialDataMap.end();
bit++) {
491 bit->second.lInear = linear;
494 damper.constitutiveEquationMap.insert(
500 CHKERR damper.setOperators(3);
504 elastic.getLoopFeEnergy(),
505 inertia.getLoopFeEnergy());
520 if (!check_if_spatial_field_exist) {
522 "MESH_NODE_POSITIONS");
538 "FLUID_PRESSURE_FE");
543 if (is_partitioned) {
562 "FLUID_PRESSURE_FE");
572 if (is_partitioned) {
591 CHKERR TSCreate(PETSC_COMM_WORLD, &ts);
592 CHKERR TSSetType(ts, TSBEULER);
598 ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>(
"Kuu",
600 CHKERR MatDuplicate(shellAij_ctx->
K, MAT_DO_NOT_COPY_VALUES,
604 D,
"DYNAMICS",
COL, shellAij_ctx->
u,
"Kuu",
COL,
607 D,
"DYNAMICS",
"VELOCITY",
COL, shellAij_ctx->
v,
"Kuu",
"DISPLACEMENT",
615 problem_ptr->
getNbDofsRow(), (
void *)shellAij_ctx, &shell_Aij);
616 CHKERR MatShellSetOperation(shell_Aij, MATOP_MULT,
618 CHKERR MatShellSetOperation(
619 shell_Aij, MATOP_ZERO_ENTRIES,
622 ConvectiveMassElement::ShellMatrixElement shell_matrix_element(m_field);
624 m_field,
"DISPLACEMENT", shellAij_ctx->
barK, PETSC_NULL, PETSC_NULL);
627 shell_matrix_element.problemName =
"Kuu";
628 shell_matrix_element.shellMatCtx = shellAij_ctx;
629 shell_matrix_element.DirichletBcPtr = &shell_dirichlet_bc;
630 shell_matrix_element.loopK.push_back(
631 ConvectiveMassElement::ShellMatrixElement::PairNameFEMethodPtr(
632 "ELASTIC", &elastic.getLoopFeLhs()));
634 shell_matrix_element.loopK.push_back(
635 ConvectiveMassElement::ShellMatrixElement::PairNameFEMethodPtr(
636 "ELASTIC", &damper.feLhs));
638 CHKERR inertia.addHOOpsVol();
639 CHKERR inertia.setShellMatrixMassOperators(
"VELOCITY",
"DISPLACEMENT",
640 "MESH_NODE_POSITIONS", linear);
644 shell_matrix_element.loopM.push_back(
645 ConvectiveMassElement::ShellMatrixElement::PairNameFEMethodPtr(
646 "ELASTIC", &inertia.getLoopFeMassLhs()));
649 shell_matrix_element.loopAuxM.push_back(
650 ConvectiveMassElement::ShellMatrixElement::PairNameFEMethodPtr(
651 "ELASTIC", &inertia.getLoopFeMassAuxLhs()));
655 shell_matrix_residual.shellMatCtx = shellAij_ctx;
658 boost::ptr_map<std::string, NeumannForcesSurface> surface_forces;
660 string fe_name_str =
"FORCE_FE";
663 surface_forces.at(fe_name_str).getLoopFe(),
false,
667 CHKERR surface_forces.at(fe_name_str)
668 .addForce(
"DISPLACEMENT", PETSC_NULL, it->getMeshsetId(),
true);
669 surface_forces.at(fe_name_str)
674 boost::ptr_map<std::string, NeumannForcesSurface> surface_pressure;
676 string fe_name_str =
"PRESSURE_FE";
679 surface_pressure.at(fe_name_str).getLoopFe(),
false,
683 CHKERR surface_pressure.at(fe_name_str)
684 .addPressure(
"DISPLACEMENT", PETSC_NULL, it->getMeshsetId(),
true);
685 surface_pressure.at(fe_name_str)
691 boost::ptr_map<std::string, EdgeForce> edge_forces;
693 string fe_name_str =
"FORCE_FE";
694 edge_forces.insert(fe_name_str,
new EdgeForce(m_field));
697 CHKERR edge_forces.at(fe_name_str)
698 .addForce(
"DISPLACEMENT", PETSC_NULL, it->getMeshsetId(),
true);
699 edge_forces.at(fe_name_str).methodsOp.push_back(
new TimeForceScale());
704 boost::ptr_map<std::string, NodalForce> nodal_forces;
706 string fe_name_str =
"FORCE_FE";
707 nodal_forces.insert(fe_name_str,
new NodalForce(m_field));
710 CHKERR nodal_forces.at(fe_name_str)
711 .addForce(
"DISPLACEMENT",
F, it->getMeshsetId(),
true);
712 nodal_forces.at(fe_name_str).methodsOp.push_back(
new TimeForceScale());
718 m_field, ts,
"VELOCITY",
"DISPLACEMENT");
731 auto add_static_rhs = [&](
auto &loops_to_do_Rhs) {
733 loops_to_do_Rhs.push_back(
735 for (
auto fit = surface_forces.begin(); fit != surface_forces.end();
737 loops_to_do_Rhs.push_back(
740 for (
auto fit = surface_pressure.begin(); fit != surface_pressure.end();
742 loops_to_do_Rhs.push_back(
745 for (
auto fit = edge_forces.begin(); fit != edge_forces.end(); fit++) {
746 loops_to_do_Rhs.push_back(
749 for (
auto fit = nodal_forces.begin(); fit != nodal_forces.end(); fit++) {
750 loops_to_do_Rhs.push_back(
754 "FLUID_PRESSURE_FE", &fluid_pressure_fe.getLoopFe()));
758 CHKERR add_static_rhs(loops_to_do_Rhs);
761 loops_to_do_Rhs.push_back(
778 loopsMonitor.push_back(
780 loopsMonitor.push_back(
789 CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
791 CHKERR TSSetFromOptions(ts);
794 CHKERR TSGetSNES(ts, &snes);
797 CHKERR SNESGetKSP(snes, &ksp);
798 CHKERR KSPSetFromOptions(ksp);
800 CHKERR KSPGetPC(ksp, &pc);
801 CHKERR PCSetType(pc, PCSHELL);
803 CHKERR PCShellSetContext(pc, (
void *)&pc_shell_ctx);
809 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
810 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
812 "DYNAMICS",
COL,
D, INSERT_VALUES, SCATTER_REVERSE);
815 if (is_solve_at_time_zero) {
817 Mat Aij = shellAij_ctx->
K;
825 "Kuu",
COL,
D, INSERT_VALUES, SCATTER_FORWARD);
826 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
827 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
829 SnesCtx snes_ctx(m_field,
"Kuu");
832 CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
833 CHKERR SNESSetApplicationContext(snes, &snes_ctx);
836 CHKERR SNESSetFromOptions(snes);
842 snes_ctx.getComputeRhs();
843 snes_ctx.getPreProcComputeRhs().push_back(&my_dirichlet_bc);
844 fluid_pressure_fe.getLoopFe().ts_t = 0;
845 CHKERR add_static_rhs(loops_to_do_Rhs);
846 snes_ctx.getPostProcComputeRhs().push_back(&my_dirichlet_bc);
849 snes_ctx.getSetOperators();
850 snes_ctx.getPreProcSetOperators().push_back(&my_dirichlet_bc);
851 loops_to_do_Mat.push_back(
853 snes_ctx.getPostProcSetOperators().push_back(&my_dirichlet_bc);
861 "Kuu",
COL,
D, INSERT_VALUES, SCATTER_FORWARD);
863 CHKERR SNESSolve(snes, PETSC_NULL,
D);
865 CHKERR SNESGetIterationNumber(snes, &its);
866 MOFEM_LOG_C(
"DYNAMIC", Sev::inform,
"number of Newton iterations = %d\n",
871 "Kuu",
COL,
D, INSERT_VALUES, SCATTER_REVERSE);
875 CHKERR SNESDestroy(&snes);
878 if (is_solve_at_time_zero) {
880 "DYNAMICS",
COL,
D, INSERT_VALUES, SCATTER_FORWARD);
881 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
882 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
886 #if PETSC_VERSION_GE(3, 7, 0)
887 CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
890 CHKERR TSGetTime(ts, &ftime);
892 PetscInt steps, snesfails, rejects, nonlinits, linits;
893 CHKERR TSGetTimeStepNumber(ts, &steps);
894 CHKERR TSGetSNESFailures(ts, &snesfails);
895 CHKERR TSGetStepRejections(ts, &rejects);
896 CHKERR TSGetSNESIterations(ts, &nonlinits);
897 CHKERR TSGetKSPIterations(ts, &linits);
899 "steps %d (%d rejected, %D SNES fails), ftime %g, nonlinits "
901 steps, rejects, snesfails, ftime, nonlinits, linits);
906 CHKERR MatDestroy(&shellAij_ctx->
K);
907 CHKERR MatDestroy(&shellAij_ctx->
M);
910 CHKERR MatDestroy(&shell_Aij);