118 {
119
121 "-pc_type lu \n"
122 "-pc_factor_mat_solver_type mumps \n"
123 "-mat_mumps_icntl_20 0 \n"
124 "-ksp_monitor \n"
125 "-snes_type newtonls \n"
126 "-snes_linesearch_type basic \n"
127 "-snes_max_it 100 \n"
128 "-snes_atol 1e-8 \n"
129 "-snes_rtol 1e-8 \n"
130 "-snes_monitor \n"
131 "-ts_monitor \n"
132 "-ts_type beuler \n"
133 "-ts_exact_final_time stepover \n";
134
136 if (!
static_cast<bool>(ifstream(
param_file))) {
137 std::ofstream file(
param_file.c_str(), std::ios::ate);
138 if (file.is_open()) {
140 file.close();
141 }
142 }
143
145
146 auto core_log = logging::core::get();
147 core_log->add_sink(
151
152 try {
153
154 PetscBool flg = PETSC_TRUE;
158 if(flg != PETSC_TRUE) {
160 "*** ERROR -my_file (MESH FILE NEEDED)");
161 }
162
163 char time_data_file_for_ground_surface[255];
164 PetscBool ground_temperature_analysis = PETSC_FALSE;
166 time_data_file_for_ground_surface,255,&ground_temperature_analysis);
167 if(ground_temperature_analysis) {
168#ifndef WITH_ADOL_C
170 "*** ERROR to do ground thermal analysis MoFEM need to be compiled "
171 "with ADOL-C");
172#endif
173 }
174
175
176 moab::Core mb_instance;
177 moab::Interface& moab = mb_instance;
178 const char *option;
179 option = "";
181
184
185 DMType dm_name = "DMTHERMAL";
187
188 DM dm;
189 CHKERR DMCreate(PETSC_COMM_WORLD, &dm);
190 CHKERR DMSetType(dm, dm_name);
191
192
193
195 bit_level0.set(0);
197 bit_level0);
198
199
204
205
206
209 );
210
211
213
216
219
220 if (flg != PETSC_TRUE) {
222 }
223
224
225
226
231
236
237
239 "MESH_NODE_POSITIONS");
244
245
246
247 PetscBool block_config;
248 char block_config_file[255];
250 block_config_file, 255, &block_config);
251 std::map<int,BlockOptionData> block_data;
252 bool solar_radiation = false;
253 if (block_config) {
254 try {
255 ifstream ini_file(block_config_file);
256
257 po::variables_map vm;
258 po::options_description config_file_options;
260
261 std::ostringstream str_order;
262 str_order << "block_" << it->getMeshsetId() << ".temperature_order";
263 config_file_options.add_options()(
264 str_order.str().c_str(),
265 po::value<int>(&block_data[it->getMeshsetId()].oRder)
266 ->default_value(
order));
267
268 std::ostringstream str_cond;
269 str_cond << "block_" << it->getMeshsetId() << ".heat_conductivity";
270 config_file_options.add_options()(
271 str_cond.str().c_str(),
272 po::value<double>(&block_data[it->getMeshsetId()].cOnductivity)
273 ->default_value(-1));
274
275 std::ostringstream str_capa;
276 str_capa << "block_" << it->getMeshsetId() << ".heat_capacity";
277 config_file_options.add_options()(
278 str_capa.str().c_str(),
279 po::value<double>(&block_data[it->getMeshsetId()].cApacity)
280 ->default_value(-1));
281
282 std::ostringstream str_init_temp;
283 str_init_temp << "block_" << it->getMeshsetId()
284 << ".initial_temperature";
285 config_file_options.add_options()(
286 str_init_temp.str().c_str(),
287 po::value<double>(&block_data[it->getMeshsetId()].initTemp)
288 ->default_value(0));
289 }
290 config_file_options.add_options()(
291 "climate_model.solar_radiation",
292 po::value<bool>(&solar_radiation)->default_value(false));
293
294 po::parsed_options parsed =
295 parse_config_file(ini_file, config_file_options, true);
296 store(parsed,vm);
297 po::notify(vm);
298
300 if (block_data[it->getMeshsetId()].oRder == -1)
301 continue;
302 if (block_data[it->getMeshsetId()].oRder ==
order)
303 continue;
304 PetscPrintf(PETSC_COMM_WORLD, "Set block %d oRder to %d\n",
305 it->getMeshsetId(), block_data[it->getMeshsetId()].oRder);
307 CHKERR moab.get_entities_by_handle(it->meshset, block_ents,
true);
308 Range ents_to_set_order;
309 CHKERR moab.get_adjacencies(block_ents, 3,
false, ents_to_set_order,
310 moab::Interface::UNION);
311 ents_to_set_order = ents_to_set_order.subset_by_type(MBTET);
312 CHKERR moab.get_adjacencies(block_ents, 2,
false, ents_to_set_order,
313 moab::Interface::UNION);
314 CHKERR moab.get_adjacencies(block_ents, 1,
false, ents_to_set_order,
315 moab::Interface::UNION);
317 block_data[it->getMeshsetId()].oRder);
319 block_data[it->getMeshsetId()].oRder);
320 }
321 std::vector<std::string> additional_parameters;
322 additional_parameters =
323 collect_unrecognized(parsed.options, po::include_positional);
324 for (std::vector<std::string>::iterator vit =
325 additional_parameters.begin();
326 vit != additional_parameters.end(); vit++) {
327 CHKERR PetscPrintf(PETSC_COMM_WORLD,
328 "** WARRING Unrecognised option %s\n", vit->c_str());
329 }
330
331 } catch (const std::exception& ex) {
332 std::ostringstream ss;
333 ss << ex.what() << std::endl;
335 }
336 }
337
338
340 CHKERR thermal_elements.addThermalElements(
"TEMP");
341 CHKERR thermal_elements.addThermalFluxElement(
"TEMP");
342 CHKERR thermal_elements.addThermalConvectionElement(
"TEMP");
343 CHKERR thermal_elements.addThermalRadiationElement(
"TEMP");
344
346 "TEMP_RATE");
347
348
350 false, false, false);
352 false, false, false);
353 CHKERR thermal_elements.setTimeSteppingProblem(
"TEMP",
"TEMP_RATE");
354
355
356 std::map<int, ThermalElement::BlockData>::iterator mit;
357 mit = thermal_elements.setOfBlocks.begin();
358 for (; mit != thermal_elements.setOfBlocks.end(); mit++) {
359
360
361
362 if (block_data[mit->first].cOnductivity != -1) {
363 PetscPrintf(PETSC_COMM_WORLD, "Set block %d heat conductivity to %3.2e\n",
364 mit->first, block_data[mit->first].cOnductivity);
365 for (
int dd = 0;
dd < 3;
dd++) {
366 mit->second.cOnductivity_mat(dd, dd) =
367 block_data[mit->first].cOnductivity;
368 }
369 }
370 if (block_data[mit->first].cApacity != -1) {
371 PetscPrintf(PETSC_COMM_WORLD, "Set block %d heat capacity to %3.2e\n",
372 mit->first, block_data[mit->first].cApacity);
373 mit->second.cApacity = block_data[mit->first].cApacity;
374 }
375 }
376
377#ifdef __GROUND_SURFACE_TEMPERATURE_HPP
381 if (ground_temperature_analysis) {
382 CHKERR ground_surface.addSurfaces(
"TEMP");
383 CHKERR ground_surface.setOperators(&time_data,
"TEMP");
384 }
385#endif
386
387
388
389
391
393 "MESH_NODE_POSITIONS");
395
396
397 mit = thermal_elements.setOfBlocks.begin();
398 for (; mit != thermal_elements.setOfBlocks.end(); mit++) {
399 if (mit->second.initTemp != 0) {
401 CHKERR moab.get_connectivity(mit->second.tEts, vertices,
true);
403 mit->second.initTemp, MBVERTEX, vertices, "TEMP");
404 }
405 }
406
408 if (std::regex_match(it->getName(), std::regex("INT_THERMAL(.*)"))) {
409 std::vector<double> data;
410 CHKERR it->getAttributes(data);
411 if (data.size() != 1)
412 SETERRQ(PETSC_COMM_SELF, 1, "Data inconsistency");
413 Range block_ents, block_verts;
414 CHKERR moab.get_entities_by_handle(it->getMeshset(), block_ents,
true);
415 CHKERR moab.get_connectivity(block_ents, block_verts,
true);
417 block_verts, "TEMP");
418 }
419 }
420
422 if (block_data[it->getMeshsetId()].initTemp != 0) {
424 CHKERR moab.get_entities_by_handle(it->meshset, block_ents,
true);
426 CHKERR moab.get_connectivity(block_ents, vertices,
true);
428 block_data[it->getMeshsetId()].initTemp, MBVERTEX, vertices, "TEMP");
429 }
430 }
431
432 MPI_Comm moab_comm_world;
433 MPI_Comm_dup(PETSC_COMM_WORLD, &moab_comm_world);
434 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
435 if (pcomm == NULL)
436 pcomm = new ParallelComm(&moab, moab_comm_world);
437
438 PetscBool is_partitioned = PETSC_FALSE;
440 &is_partitioned, PETSC_NULL);
441
442 Range thermal_element_ents;
444 thermal_element_ents);
447 CHKERR skin.find_skin(0, thermal_element_ents,
false, skin_faces);
449 if (is_partitioned) {
450 CHKERR pcomm->filter_pstatus(skin_faces,
451 PSTATUS_SHARED | PSTATUS_MULTISHARED,
452 PSTATUS_NOT, -1, &proc_skin);
453 } else {
454 proc_skin = skin_faces;
455 }
461 "MESH_NODE_POSITIONS");
463 "POST_PROC_SKIN");
464
465
467
469
470
474
475
476
477
479 }
480
481 std::vector<std::array<double, 3>> eval_points;
482 eval_points.resize(0);
483 PetscBool eval_points_flg = PETSC_FALSE;
484 char eval_points_file[255];
486 eval_points_file, 255, &eval_points_flg);
487 if (eval_points_flg) {
488 std::ifstream in_file(eval_points_file, std::ios::in);
489 if (!in_file.is_open()) {
491 eval_points_file);
492 }
493 double x, y, z;
494 while (in_file >> x >> y >> z) {
495 eval_points.push_back({x, y, z});
496 }
497 }
498
499
501 CHKERR DMSetFromOptions(dm);
502
508#ifdef __GROUND_SURFACE_TEMPERATURE_HPP
509 if (ground_temperature_analysis) {
511 }
512#endif
513
515
516
522
525 "TEMP_RATE");
527 eval_points);
529
530
535
536
538 NULL);
541#ifdef __GROUND_SURFACE_TEMPERATURE_HPP
543 &exectuteGenericClimateModel, NULL);
544 {
545 if (solar_radiation) {
546 boost::ptr_vector<
548 hi_it;
549 it = ground_surface.preProcessShade.begin();
550 hi_it = ground_surface.preProcessShade.end();
551 for (; it != hi_it; it++) {
553 }
554 }
555 }
556#endif
557
558
560 NULL);
562 NULL, NULL);
564 &thermal_elements.feConvectionRhs, NULL, NULL);
566 &thermal_elements.feRadiationRhs, NULL, NULL);
567#ifdef __GROUND_SURFACE_TEMPERATURE_HPP
568 if (ground_temperature_analysis) {
570 &ground_surface.getFeGroundSurfaceRhs(), NULL,
571 NULL);
572 }
573#endif
574
575
577 NULL);
579 &thermal_elements.feConvectionLhs, NULL, NULL);
581 &thermal_elements.feRadiationLhs, NULL, NULL);
582#ifdef __GROUND_SURFACE_TEMPERATURE_HPP
583 if (ground_temperature_analysis) {
585 &ground_surface.getFeGroundSurfaceLhs(), NULL,
586 NULL);
587 }
588#endif
589
590
593
596
599
600
601 TS ts;
602 CHKERR TSCreate(PETSC_COMM_WORLD,&ts);
603 CHKERR TSSetType(ts,TSBEULER);
604
605 CHKERR TSSetIFunction(ts,
F,PETSC_NULL,PETSC_NULL);
606 CHKERR TSSetIJacobian(ts,A,A,PETSC_NULL,PETSC_NULL);
607
609
611
613
614 double ftime = 1;
615 CHKERR TSSetDuration(ts,PETSC_DEFAULT,ftime);
616 CHKERR TSSetFromOptions(ts);
618
620 CHKERR TSGetTime(ts,&ftime);
621
622
624
625 PetscInt steps,snesfails,rejects,nonlinits,linits;
626 CHKERR TSGetTimeStepNumber(ts,&steps);
627 CHKERR TSGetSNESFailures(ts,&snesfails);
628 CHKERR TSGetStepRejections(ts,&rejects);
629 CHKERR TSGetSNESIterations(ts,&nonlinits);
630 CHKERR TSGetKSPIterations(ts,&linits);
631
632 PetscPrintf(PETSC_COMM_WORLD,
633 "steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits %D, "
634 "linits %D\n",
635 steps, rejects, snesfails, ftime, nonlinits, linits);
636
637
638
639 PetscBool save_solution = PETSC_TRUE;
641 &save_solution, PETSC_NULL);
642
643 if (save_solution) {
644 if (is_partitioned) {
645 CHKERR moab.write_file(
"solution.h5m");
646 } else {
648 CHKERR moab.write_file(
"solution.h5m");
649 }
650 }
651 }
652
656
658
659 }
661
663}
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