1935 {
1937
1941
1942 auto preProc =
1944
1945 auto create_custom_ts = [&]() {
1946 auto set_dm_section = [&](auto dm) {
1948 PetscSection section;
1950 simple->getProblemName(), §ion);
1951 CHKERR DMSetDefaultSection(dm, section);
1952 CHKERR DMSetDefaultGlobalSection(dm, section);
1953 CHKERR PetscSectionDestroy(§ion);
1955 };
1956 auto dm =
simple->getDM();
1957
1958 CHKERR set_dm_section(dm);
1959 boost::shared_ptr<FEMethod> null;
1960 preProc->methodsOpForMove = boost::shared_ptr<MethodForForceScaling>(
1962
1963 preProc->methodsOpForRollersDisp = boost::shared_ptr<MethodForForceScaling>(
1965
1966
1967
1968
1969
1970
1971
1972
1973
1974
1975
1976
1977
1978
1979
1980
1981
1982
1983
1984
1985 if (true) {
1986
1987
1991 null);
1992
1993 auto postproc_bound_el = [&]() {
1996 auto &bmc = *bmc_ptr;
1997 bmc.clear();
1998
1999
2000
2001
2002 }
2004 };
2005
2007 pipeline_mng->getBoundaryLhsFE()->postProcessHook = postproc_bound_el;
2008
2009 if (pipeline_mng->getBoundaryLhsFE())
2011 pipeline_mng->getBoundaryLhsFE(), null,
2012 pipeline_mng->getBoundaryLhsFE());
2013
2015 pipeline_mng->getDomainLhsFE(), null, null);
2018 null, null);
2019
2020 if (pipeline_mng->getSkeletonLhsFE())
2022 pipeline_mng->getSkeletonLhsFE(), null,
2023 null);
2024
2025
2026
2028 preProc);
2029
2030 if (pipeline_mng->getDomainRhsFE()) {
2031
2032
2036 preProc, null);
2040
2042 pipeline_mng->getDomainRhsFE(), null,
2043 null);
2044
2045 if (pipeline_mng->getBoundaryRhsFE())
2047 pipeline_mng->getBoundaryRhsFE(), null,
2048 null);
2049 if (pipeline_mng->getSkeletonRhsFE())
2051 pipeline_mng->getSkeletonRhsFE(), null,
2052 null);
2053
2055 fit++) {
2057 &fit->second->getLoopFe(), NULL, NULL);
2058 }
2059
2060
2063 &fit->second->getLoopFe(), NULL, NULL);
2064 }
2065
2066
2068 fit++) {
2070 &fit->second->getLoopFe(), NULL, NULL);
2071 }
2072
2073
2074
2075
2077 preProc);
2078 }
2079 }
2080
2081 auto print_active_set = [&](std::array<int, 2> &lgauss_pts, string name,
2084 std::vector<int> gauss_pts(2, 0);
2086 CHKERR MPIU_Allreduce(lgauss_pts.data(), gauss_pts.data(), 2, MPIU_INT,
2087 MPIU_SUM, PetscObjectComm((PetscObject)dm));
2088
2089 MOFEM_LOG_C(
"WORLD", Sev::inform,
"\t \t Active %s gauss pts: %d / %d",
2090 name.c_str(), (int)gauss_pts[0], (int)gauss_pts[1]);
2091 lgauss_pts[0] = lgauss_pts[1] = 0;
2093 };
2094
2095 auto postproc_dom = [&]() {
2100 };
2101 auto postproc_bound = [&]() {
2105 };
2106
2109 pipeline_mng->getDomainRhsFE()->postProcessHook = postproc_dom;
2110
2112 pipeline_mng->getBoundaryRhsFE()->postProcessHook = postproc_bound;
2113 }
2114
2117 return ts;
2118 };
2119
2120 auto solver = create_custom_ts();
2121
2122 CHKERR TSSetExactFinalTime(solver, TS_EXACTFINALTIME_MATCHSTEP);
2123
2124 auto dm =
simple->getDM();
2126
2128 CHKERR TSSetFromOptions(solver);
2129 CHKERR TSSetSolution(solver,
D);
2130
2132
2133
2134 if (true) {
2135 SNES snes;
2136 CHKERR TSGetSNES(solver, &snes);
2137 KSP ksp;
2138 CHKERR SNESGetKSP(snes, &ksp);
2139 PC pC;
2140 CHKERR KSPGetPC(ksp, &pC);
2141 PetscBool is_pcfs = PETSC_FALSE;
2142 PetscObjectTypeCompare((PetscObject)pC, PCFIELDSPLIT, &is_pcfs);
2144
2145 auto make_is_field_map = [&]() {
2146 PetscSection section;
2147 CHKERR DMGetDefaultSection(dm, §ion);
2148
2149 int num_fields;
2150 CHKERR PetscSectionGetNumFields(section, &num_fields);
2151
2154
2155 map<std::string, SmartPetscObj<IS>> is_map;
2156 for (int ff = 0; ff != num_fields; ff++) {
2157
2160
2164 };
2165 return is_map;
2166 };
2167 auto is_map = make_is_field_map();
2168
2169 auto set_fieldsplit_on_bc = [&](PC &bc_pc, string name_prb) {
2171
2172 PetscBool is_bc_field_split;
2173 PetscObjectTypeCompare((PetscObject)bc_pc, PCFIELDSPLIT,
2174 &is_bc_field_split);
2175
2176 auto block_prefix =
simple->getProblemName();
2178
2179 auto is_all_bc =
2180 bc_mng->
getBlockIS(block_prefix,
"MOVE_Z",
"U", name_prb, 2, 2);
2181 is_all_bc = bc_mng->
getBlockIS(block_prefix,
"MOVE_Y",
"U", name_prb, 1,
2182 1, is_all_bc);
2183 is_all_bc = bc_mng->
getBlockIS(block_prefix,
"MOVE_X",
"U", name_prb, 0,
2184 0, is_all_bc);
2185 is_all_bc = bc_mng->
getBlockIS(block_prefix,
"ROTATE_X",
"U", name_prb, 1,
2186 2, is_all_bc);
2187 is_all_bc = bc_mng->
getBlockIS(block_prefix,
"ROTATE_Z",
"U", name_prb, 0,
2188 1, is_all_bc);
2189 is_all_bc = bc_mng->
getBlockIS(block_prefix,
"ROTATE_ALL",
"U", name_prb,
2190 0, 2, is_all_bc);
2191 is_all_bc = bc_mng->
getBlockIS(block_prefix,
"MOVE_ALL",
"U", name_prb, 0,
2192 2, is_all_bc);
2193
2194 int is_all_bc_size;
2195 CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
2196
2197 if (is_bc_field_split && is_all_bc_size) {
2199 << "Fieldsplit preconditioner for boundary block size: "
2200 << is_all_bc_size;
2201 CHKERR PCFieldSplitSetIS(bc_pc, NULL,
2202 is_all_bc);
2203
2204
2205
2206
2207
2208
2209
2210
2211
2212
2213
2214
2215
2217 }
2218
2219
2220
2221
2222
2223
2225 };
2226
2227 auto set_global_mat_and_vec = [&]() {
2231 CHKERR TSSetIFunction(solver, fe->ts_F, PETSC_NULL, PETSC_NULL);
2232 CHKERR TSSetIJacobian(solver, fe->ts_B, fe->ts_B, PETSC_NULL, PETSC_NULL);
2233 CHKERR KSPSetOperators(ksp, fe->ts_B, fe->ts_B);
2235 };
2236
2237
2239
2240 CHKERR set_global_mat_and_vec();
2241
2242 PetscBool is_l0_field_split = PETSC_TRUE;
2243 PetscBool is_l1_field_split = PETSC_FALSE;
2244
2245
2246
2247 auto set_fieldsplit_contact = [&]() {
2249
2253
2254 is_map["E_IS_SIG"] = sig_data->rowIs;
2255 sigma_ao = sig_data->getSmartRowMap();
2256
2257 CHKERR PCFieldSplitSetIS(pC, NULL, is_map[
"SIGMA"]);
2258 CHKERR PCFieldSplitSetIS(pC, NULL, is_map[
"E_IS_SIG"]);
2259
2260
2261
2263
2264
2265
2266
2268 KSP *ct_ksp;
2269 CHKERR PCFieldSplitGetSubKSP(pC, &
n, &ct_ksp);
2270 PC l1_pc;
2271 CHKERR KSPGetPC(ct_ksp[1], &l1_pc);
2272 PetscObjectTypeCompare((PetscObject)l1_pc, PCFIELDSPLIT,
2273 &is_l1_field_split);
2274
2275
2276 pC = l1_pc;
2277 CHKERR PetscFree(ct_ksp);
2278 }
2280 };
2281
2282
2283
2285 CHKERR set_fieldsplit_contact();
2286
2287 PetscObjectTypeCompare((PetscObject)pC, PCFIELDSPLIT, &is_l1_field_split);
2289
2293
2294 is_map["E_IS_EP"] = ep_data->rowIs;
2295
2296 if (sigma_ao)
2297 CHKERR AOApplicationToPetscIS(sigma_ao, is_map[
"EP"]);
2298
2299
2300 CHKERR PCFieldSplitSetIS(pC,
"ep", is_map[
"EP"]);
2301 CHKERR PCFieldSplitSetIS(pC,
"etau", is_map[
"E_IS_EP"]);
2302
2303 PCCompositeType pc_type;
2304 CHKERR PCFieldSplitGetType(pC, &pc_type);
2305
2306 if (pc_type == PC_COMPOSITE_SCHUR) {
2309 if (sigma_ao)
2311
2313 boost::make_shared<BlockMatContainer>();
2316
2317 CHKERR PCFieldSplitSetSchurPre(pC, PC_FIELDSPLIT_SCHUR_PRE_USER,
2319
2321 }
2322
2324 KSP *tt_ksp;
2325 CHKERR PCFieldSplitGetSubKSP(pC, &
n, &tt_ksp);
2326 PC tau_pc;
2327 CHKERR KSPGetPC(tt_ksp[1], &tau_pc);
2328 PetscBool is_tau_field_split;
2329 PetscObjectTypeCompare((PetscObject)tau_pc, PCFIELDSPLIT,
2330 &is_tau_field_split);
2331
2332 if (is_tau_field_split) {
2336
2337
2338 is_map["E_IS_TAU"] = tau_data->rowIs;
2339
2340 AO ep_ao;
2341 CHKERR ep_data->getRowMap(&ep_ao);
2342
2343 if (sigma_ao)
2344 CHKERR AOApplicationToPetscIS(sigma_ao, is_map[
"TAU"]);
2345 CHKERR AOApplicationToPetscIS(ep_ao, is_map[
"TAU"]);
2346
2347 CHKERR PCFieldSplitSetIS(tau_pc,
"tau", is_map[
"TAU"]);
2348 CHKERR PCFieldSplitSetIS(tau_pc,
"elastic", is_map[
"E_IS_TAU"]);
2349
2350 CHKERR PCFieldSplitGetType(tau_pc, &pc_type);
2351 if (pc_type == PC_COMPOSITE_SCHUR) {
2354
2355 CHKERR PCFieldSplitSetSchurPre(tau_pc,
2356 PC_FIELDSPLIT_SCHUR_PRE_USER,
2358 }
2360
2361 PetscInt bc_n;
2362 KSP *bc_ksp;
2363 CHKERR PCFieldSplitGetSubKSP(tau_pc, &bc_n, &bc_ksp);
2364 PC bc_pc;
2365 CHKERR KSPGetPC(bc_ksp[1], &bc_pc);
2366
2367 CHKERR set_fieldsplit_on_bc(bc_pc, schur_tau_ptr->
getName());
2368 CHKERR PetscFree(bc_ksp);
2369 }
2370 }
2371 CHKERR PetscFree(tt_ksp);
2372 }
2373 }
2375
2376 char mumps_options[] = "-fieldsplit_0_mat_mumps_icntl_14 1200 "
2377 "-fieldsplit_0_mat_mumps_icntl_24 1 "
2378 "-fieldsplit_0_mat_mumps_icntl_20 0 "
2379 "-fieldsplit_0_mat_mumps_icntl_13 1 "
2380 "-fieldsplit_1_mat_mumps_icntl_14 1200 "
2381 "-fieldsplit_1_mat_mumps_icntl_24 1 "
2382 "-fieldsplit_1_mat_mumps_icntl_20 0 "
2383 "-fieldsplit_1_mat_mumps_icntl_13 1";
2384 CHKERR PetscOptionsInsertString(NULL, mumps_options);
2385 CHKERR set_global_mat_and_vec();
2386 CHKERR set_fieldsplit_on_bc(pC,
simple->getProblemName());
2387
2388 }
2389
2390
2391
2392
2393
2394
2395
2396
2397
2398
2399
2400 }
2401
2402 auto set_section_monitor = [&]() {
2404 SNES snes;
2405 CHKERR TSGetSNES(solver, &snes);
2406 PetscViewerAndFormat *vf;
2407 CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
2408 PETSC_VIEWER_DEFAULT, &vf);
2410 snes,
2411 (
MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
void *))SNESMonitorFields,
2414 };
2415
2416 auto create_post_process_element = [&]() {
2425
2431 } else
2443 } else
2446
2451 else
2455 }
2463
2467 }
2472 else
2476
2481
2483
2487 postProcFe->addFieldValuesPostProc(
"U",
"DISPLACEMENT");
2489 postProcFe->addFieldValuesPostProc(
"MESH_NODE_POSITIONS");
2491 };
2492
2493 auto scatter_create = [&](auto coeff) {
2495 CHKERR is_manager->isCreateProblemFieldAndRank(
simple->getProblemName(),
2496 ROW,
"U", coeff, coeff, is);
2497 int loc_size;
2498 CHKERR ISGetLocalSize(is, &loc_size);
2501 VecScatter scatter;
2502 CHKERR VecScatterCreate(
D, is,
v, PETSC_NULL, &scatter);
2505 };
2506
2507 auto set_time_monitor = [&]() {
2509 boost::shared_ptr<MMonitor> monitor_ptr(
2512 boost::shared_ptr<ForcesAndSourcesCore> null;
2514 monitor_ptr, null, null);
2516 };
2517
2518 CHKERR set_section_monitor();
2519 CHKERR create_post_process_element();
2520
2522
2523 CHKERR PetscPrintf(PETSC_COMM_WORLD,
2524 "Reading vector in binary from %s file...\n",
2526 PetscViewer viewer;
2527 CHKERR PetscViewerBinaryOpen(PETSC_COMM_WORLD,
2529 &viewer);
2531
2532 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
2533 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
2535
2536 typedef boost::tokenizer<boost::char_separator<char>> Tokenizer;
2538 Tokenizer tok{s, boost::char_separator<char>("_")};
2539 auto it = tok.begin();
2540 const int restart_step = std::stoi((++it)->c_str());
2542 string restart_tt = *(++it);
2543 restart_tt.erase(restart_tt.length() - 4);
2544 const double restart_time = std::atof(restart_tt.c_str());
2545 CHKERR TSSetStepNumber(solver, restart_step);
2546 CHKERR TSSetTime(solver, restart_time);
2547 }
2548
2552 CHKERR set_time_monitor();
2553
2554
2555
2556
2557
2558
2559
2560
2561
2562
2563
2564
2565
2566
2567
2568
2569
2571
2572 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
2573 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
2575
2577}
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
FTensor::Index< 'n', SPACE_DIM > n
PetscErrorCode DMMoFEMTSSetIFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set TS implicit function evaluation function
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
PetscErrorCode DMMoFEMTSSetIJacobian(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS Jacobian evaluation function
const double v
phase velocity of light in medium (cm/ns)
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
DEPRECATED auto smartCreateDMVector(DM dm)
boost::weak_ptr< BlockMatContainer > block_mat_container_ptr
auto createTS(MPI_Comm comm)
constexpr auto field_name
SmartPetscObj< IS > getBlockIS(const std::string block_prefix, const std::string block_name, const std::string field_name, const std::string problem_name, int lo, int hi, SmartPetscObj< IS > is_expand=SmartPetscObj< IS >())
Get block IS.
Deprecated interface functions.
Section manager is used to create indexes and sections.
keeps basic data about problem
boost::shared_ptr< SubProblemData > subProblemData
intrusive_ptr for managing petsc objects