124 {
125
127 "-pc_type lu \n"
128 "-pc_factor_mat_solver_type mumps \n"
129 "-mat_mumps_icntl_20 0 \n"
130 "-ksp_monitor \n"
131 "-snes_type newtonls \n"
132 "-snes_linesearch_type basic \n"
133 "-snes_max_it 100 \n"
134 "-snes_atol 1e-8 \n"
135 "-snes_rtol 1e-8 \n"
136 "-snes_monitor \n"
137 "-ts_monitor \n"
138 "-ts_type beuler \n"
139 "-ts_exact_final_time stepover \n";
140
142 if (!
static_cast<bool>(ifstream(
param_file))) {
143 std::ofstream file(
param_file.c_str(), std::ios::ate);
144 if (file.is_open()) {
146 file.close();
147 }
148 }
149
151
152 auto core_log = logging::core::get();
153 core_log->add_sink(
157
158 try {
159
160 PetscBool flg = PETSC_TRUE;
164 if(flg != PETSC_TRUE) {
166 "*** ERROR -my_file (MESH FILE NEEDED)");
167 }
168
169 char time_data_file_for_ground_surface[255];
170 PetscBool ground_temperature_analysis = PETSC_FALSE;
172 time_data_file_for_ground_surface,255,&ground_temperature_analysis);
173 if(ground_temperature_analysis) {
174#ifndef WITH_ADOL_C
176 "*** ERROR to do ground thermal analysis MoFEM need to be compiled "
177 "with ADOL-C");
178#endif
179 }
180
181
182 moab::Core mb_instance;
183 moab::Interface& moab = mb_instance;
184 const char *option;
185 option = "";
187
190
191 DMType dm_name = "DMTHERMAL";
193
194 DM dm;
195 CHKERR DMCreate(PETSC_COMM_WORLD, &dm);
196 CHKERR DMSetType(dm, dm_name);
197
198
199
201 bit_level0.set(0);
203 bit_level0);
204
205
210
211
212
215 );
216
217
219
222
225
226 if (flg != PETSC_TRUE) {
228 }
229
230
231
232
237
242
243
245 "MESH_NODE_POSITIONS");
250
251
252
253 PetscBool block_config;
254 char block_config_file[255];
256 block_config_file, 255, &block_config);
257 std::map<int,BlockOptionData> block_data;
258 bool solar_radiation = false;
259 if (block_config) {
260 try {
261 ifstream ini_file(block_config_file);
262
263 po::variables_map vm;
264 po::options_description config_file_options;
266
267 std::ostringstream str_order;
268 str_order << "block_" << it->getMeshsetId() << ".temperature_order";
269 config_file_options.add_options()(
270 str_order.str().c_str(),
271 po::value<int>(&block_data[it->getMeshsetId()].oRder)
272 ->default_value(
order));
273
274 std::ostringstream str_cond;
275 str_cond << "block_" << it->getMeshsetId() << ".heat_conductivity";
276 config_file_options.add_options()(
277 str_cond.str().c_str(),
278 po::value<double>(&block_data[it->getMeshsetId()].cOnductivity)
279 ->default_value(-1));
280
281 std::ostringstream str_capa;
282 str_capa << "block_" << it->getMeshsetId() << ".heat_capacity";
283 config_file_options.add_options()(
284 str_capa.str().c_str(),
285 po::value<double>(&block_data[it->getMeshsetId()].cApacity)
286 ->default_value(-1));
287
288 std::ostringstream str_init_temp;
289 str_init_temp << "block_" << it->getMeshsetId()
290 << ".initial_temperature";
291 config_file_options.add_options()(
292 str_init_temp.str().c_str(),
293 po::value<double>(&block_data[it->getMeshsetId()].initTemp)
294 ->default_value(0));
295 }
296 config_file_options.add_options()(
297 "climate_model.solar_radiation",
298 po::value<bool>(&solar_radiation)->default_value(false));
299
300 po::parsed_options parsed =
301 parse_config_file(ini_file, config_file_options, true);
302 store(parsed,vm);
303 po::notify(vm);
304
306 if (block_data[it->getMeshsetId()].oRder == -1)
307 continue;
308 if (block_data[it->getMeshsetId()].oRder ==
order)
309 continue;
310 PetscPrintf(PETSC_COMM_WORLD, "Set block %d oRder to %d\n",
311 it->getMeshsetId(), block_data[it->getMeshsetId()].oRder);
313 CHKERR moab.get_entities_by_handle(it->meshset, block_ents,
true);
314 Range ents_to_set_order;
315 CHKERR moab.get_adjacencies(block_ents, 3,
false, ents_to_set_order,
316 moab::Interface::UNION);
317 ents_to_set_order = ents_to_set_order.subset_by_type(MBTET);
318 CHKERR moab.get_adjacencies(block_ents, 2,
false, ents_to_set_order,
319 moab::Interface::UNION);
320 CHKERR moab.get_adjacencies(block_ents, 1,
false, ents_to_set_order,
321 moab::Interface::UNION);
323 block_data[it->getMeshsetId()].oRder);
325 block_data[it->getMeshsetId()].oRder);
326 }
327 std::vector<std::string> additional_parameters;
328 additional_parameters =
329 collect_unrecognized(parsed.options, po::include_positional);
330 for (std::vector<std::string>::iterator vit =
331 additional_parameters.begin();
332 vit != additional_parameters.end(); vit++) {
333 CHKERR PetscPrintf(PETSC_COMM_WORLD,
334 "** WARRING Unrecognised option %s\n", vit->c_str());
335 }
336
337 } catch (const std::exception& ex) {
338 std::ostringstream ss;
339 ss << ex.what() << std::endl;
341 }
342 }
343
344
346 CHKERR thermal_elements.addThermalElements(
"TEMP");
347 CHKERR thermal_elements.addThermalFluxElement(
"TEMP");
348 CHKERR thermal_elements.addThermalConvectionElement(
"TEMP");
349 CHKERR thermal_elements.addThermalRadiationElement(
"TEMP");
350
352 "TEMP_RATE");
353
354
356 false, false, false);
358 false, false, false);
359 CHKERR thermal_elements.setTimeSteppingProblem(
"TEMP",
"TEMP_RATE");
360
361
362 std::map<int, ThermalElement::BlockData>::iterator mit;
363 mit = thermal_elements.setOfBlocks.begin();
364 for (; mit != thermal_elements.setOfBlocks.end(); mit++) {
365
366
367
368 if (block_data[mit->first].cOnductivity != -1) {
369 PetscPrintf(PETSC_COMM_WORLD, "Set block %d heat conductivity to %3.2e\n",
370 mit->first, block_data[mit->first].cOnductivity);
371 for (
int dd = 0;
dd < 3;
dd++) {
372 mit->second.cOnductivity_mat(dd, dd) =
373 block_data[mit->first].cOnductivity;
374 }
375 }
376 if (block_data[mit->first].cApacity != -1) {
377 PetscPrintf(PETSC_COMM_WORLD, "Set block %d heat capacity to %3.2e\n",
378 mit->first, block_data[mit->first].cApacity);
379 mit->second.cApacity = block_data[mit->first].cApacity;
380 }
381 }
382
383#ifdef __GROUND_SURFACE_TEMPERATURE_HPP
387 if (ground_temperature_analysis) {
388 CHKERR ground_surface.addSurfaces(
"TEMP");
389 CHKERR ground_surface.setOperators(&time_data,
"TEMP");
390 }
391#endif
392
393
394
395
397
399 "MESH_NODE_POSITIONS");
401
402
403 mit = thermal_elements.setOfBlocks.begin();
404 for (; mit != thermal_elements.setOfBlocks.end(); mit++) {
405 if (mit->second.initTemp != 0) {
407 CHKERR moab.get_connectivity(mit->second.tEts, vertices,
true);
409 mit->second.initTemp, MBVERTEX, vertices, "TEMP");
410 }
411 }
412
414 if (std::regex_match(it->getName(), std::regex("INT_THERMAL(.*)"))) {
415 std::vector<double> data;
416 CHKERR it->getAttributes(data);
417 if (data.size() != 1)
418 SETERRQ(PETSC_COMM_SELF, 1, "Data inconsistency");
419 Range block_ents, block_verts;
420 CHKERR moab.get_entities_by_handle(it->getMeshset(), block_ents,
true);
421 CHKERR moab.get_connectivity(block_ents, block_verts,
true);
423 block_verts, "TEMP");
424 }
425 }
426
428 if (block_data[it->getMeshsetId()].initTemp != 0) {
430 CHKERR moab.get_entities_by_handle(it->meshset, block_ents,
true);
432 CHKERR moab.get_connectivity(block_ents, vertices,
true);
434 block_data[it->getMeshsetId()].initTemp, MBVERTEX, vertices, "TEMP");
435 }
436 }
437
438 MPI_Comm moab_comm_world;
439 MPI_Comm_dup(PETSC_COMM_WORLD, &moab_comm_world);
440 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
441 if (pcomm == NULL)
442 pcomm = new ParallelComm(&moab, moab_comm_world);
443
444 PetscBool is_partitioned = PETSC_FALSE;
446 &is_partitioned, PETSC_NULL);
447
448 Range thermal_element_ents;
450 thermal_element_ents);
451
452 PetscBool save_skin = PETSC_TRUE;
454 &save_skin, PETSC_NULL);
455
458 CHKERR skin.find_skin(0, thermal_element_ents,
false, skin_faces);
460 if (is_partitioned) {
461 CHKERR pcomm->filter_pstatus(skin_faces,
462 PSTATUS_SHARED | PSTATUS_MULTISHARED,
463 PSTATUS_NOT, -1, &proc_skin);
464 } else {
465 proc_skin = skin_faces;
466 }
467
468 if (save_skin) {
474 "MESH_NODE_POSITIONS");
476 "POST_PROC_SKIN");
477 }
478
479
481
483
484
488
489
490
491
493 }
494
495 std::vector<std::array<double, 3>> eval_points;
496 eval_points.resize(0);
497 PetscBool eval_points_flg = PETSC_FALSE;
498 char eval_points_file[255];
500 eval_points_file, 255, &eval_points_flg);
501 if (eval_points_flg) {
502 std::ifstream in_file(eval_points_file, std::ios::in);
503 if (!in_file.is_open()) {
505 eval_points_file);
506 }
507 double x, y, z;
508 while (in_file >> x >> y >> z) {
509 eval_points.push_back({x, y, z});
510 }
511 }
512
513
515 CHKERR DMSetFromOptions(dm);
516
521 if (save_skin)
523
524#ifdef __GROUND_SURFACE_TEMPERATURE_HPP
525 if (ground_temperature_analysis) {
527 }
528#endif
529
531
532
538
541 "TEMP_RATE");
543 eval_points);
545
546
551
552
554 NULL);
557#ifdef __GROUND_SURFACE_TEMPERATURE_HPP
559 &exectuteGenericClimateModel, NULL);
560 {
561 if (solar_radiation) {
562 boost::ptr_vector<
564 hi_it;
565 it = ground_surface.preProcessShade.begin();
566 hi_it = ground_surface.preProcessShade.end();
567 for (; it != hi_it; it++) {
569 }
570 }
571 }
572#endif
573
574
576 NULL);
578 NULL, NULL);
580 &thermal_elements.feConvectionRhs, NULL, NULL);
582 &thermal_elements.feRadiationRhs, NULL, NULL);
583#ifdef __GROUND_SURFACE_TEMPERATURE_HPP
584 if (ground_temperature_analysis) {
586 &ground_surface.getFeGroundSurfaceRhs(), NULL,
587 NULL);
588 }
589#endif
590
591
593 NULL);
595 &thermal_elements.feConvectionLhs, NULL, NULL);
597 &thermal_elements.feRadiationLhs, NULL, NULL);
598#ifdef __GROUND_SURFACE_TEMPERATURE_HPP
599 if (ground_temperature_analysis) {
601 &ground_surface.getFeGroundSurfaceLhs(), NULL,
602 NULL);
603 }
604#endif
605
606
609
612
615
616
617 TS ts;
618 CHKERR TSCreate(PETSC_COMM_WORLD,&ts);
619 CHKERR TSSetType(ts,TSBEULER);
620
621 CHKERR TSSetIFunction(ts,
F,PETSC_NULL,PETSC_NULL);
622 CHKERR TSSetIJacobian(ts,A,A,PETSC_NULL,PETSC_NULL);
623
625
627
629
630 double ftime = 1;
631 CHKERR TSSetDuration(ts,PETSC_DEFAULT,ftime);
632 CHKERR TSSetFromOptions(ts);
634
636 CHKERR TSGetTime(ts,&ftime);
637
638
640
641 PetscInt steps,snesfails,rejects,nonlinits,linits;
642 CHKERR TSGetTimeStepNumber(ts,&steps);
643 CHKERR TSGetSNESFailures(ts,&snesfails);
644 CHKERR TSGetStepRejections(ts,&rejects);
645 CHKERR TSGetSNESIterations(ts,&nonlinits);
646 CHKERR TSGetKSPIterations(ts,&linits);
647
648 PetscPrintf(PETSC_COMM_WORLD,
649 "steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits %D, "
650 "linits %D\n",
651 steps, rejects, snesfails, ftime, nonlinits, linits);
652
653
654
655 PetscBool save_solution = PETSC_TRUE;
657 &save_solution, PETSC_NULL);
658
659 if (save_solution) {
660 if (is_partitioned) {
661 CHKERR moab.write_file(
"solution.h5m");
662 } else {
664 CHKERR moab.write_file(
"solution.h5m");
665 }
666 }
667 }
668
672
674
675 }
677
679}
const std::string default_options
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MYPCOMM_INDEX
default communicator number PCOMM
@ MOFEM_STD_EXCEPTION_THROW
@ MOFEM_DATA_INCONSISTENCY
#define CHKERR
Inline error check.
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
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 DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
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
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
virtual MoFEMErrorCode add_ents_to_finite_element_by_dim(const EntityHandle entities, const int dim, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode get_finite_element_entities_by_dimension(const std::string name, int dim, Range &ents) const =0
get entities in the finite element by dimension
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
virtual MoFEMErrorCode finalize_series_recorder(const std::string &serie_name)
virtual bool check_series(const std::string &name) const
check if series is in database
virtual MoFEMErrorCode initialize_series_recorder(const std::string &serie_name)
virtual MoFEMErrorCode add_series_recorder(const std::string &series_name)
virtual MoFEMErrorCode delete_recorder_series(const std::string &series_name)
const FTensor::Tensor2< T, Dim, Dim > Vec
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
PetscErrorCode monitor(TS ts, PetscInt step, PetscReal crtime, Vec u, void *ctx)
Implementation of ground surface temperature.
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
virtual int get_comm_rank() const =0
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmSync()
Get the strm sync object.
Projection of edge entities with one mid-node on hierarchical basis.
Interface for Time Stepping (TS) solver.
BasicMethodsSequence & getPostProcessMonitor()
Get the postProcess to do Monitor object.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
TS monitore it records temperature at time steps.
this calass is to control time stepping
structure grouping operators and data used for thermal problems