v0.14.0
Loading...
Searching...
No Matches
thermal_unsteady.cpp
Go to the documentation of this file.
1/** \file thermal_unsteady.cpp
2 \ingroup mofem_thermal_elem
3 \brief Example of thermal unsteady analyze.
4
5 TODO:
6 \todo Make it work in distributed meshes with multigird solver. At the moment
7 it is not working efficient as can.
8*/
9
10
11
13using namespace MoFEM;
14
15#ifdef __GROUND_SURFACE_TEMPERATURE_HPP
16
18 #include <GroundSurfaceTemperature.hpp>
19
20 #include <time.h>
21 extern "C" {
22 #include <spa.h>
23 }
24 #include <CrudeClimateModel.hpp>
25
26#endif // __GROUND_SURFACE_TEMPERATURE_HPP
27
30
31static char help[] =
32 "-my_file mesh file\n"
33 "-order set approx. order to all blocks\n"
34 "-my_block_config set block data\n"
35 "-my_ground_analysis_data data for crude climate model\n"
36 "\n";
37
38struct BlockOptionData {
39 int oRder;
41 double cApacity;
42 double initTemp;
44 oRder(-1),
45 cOnductivity(-1),
46 cApacity(-1),
47 initTemp(0) {}
48};
49
50struct MonitorPostProc : public FEMethod {
51
55
56 bool iNit;
57 int pRT;
58 PetscBool saveSkin;
59
61 : FEMethod(), mField(m_field), postProc(m_field), skinPostProc(m_field),
62 iNit(false) {
63
64 PetscBool flg = PETSC_TRUE;
65 CHKERR PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-my_output_prt", &pRT,
66 &flg);
67 CHKERRABORT(PETSC_COMM_WORLD, ierr);
68 if (flg != PETSC_TRUE) {
69 pRT = 1;
70 }
71 saveSkin = PETSC_TRUE;
72 CHKERR PetscOptionsGetBool(PETSC_NULL, PETSC_NULL, "-my_save_skin",
73 &saveSkin, PETSC_NULL);
74 }
75
79 }
80
84 }
85
88
89 if (!iNit) {
90 CHKERR addHOOpsVol("MESH_NODE_POSITIONS", postProc, true, false, false,
91 false);
96 CHKERR postProc.addFieldValuesPostProc("MESH_NODE_POSITIONS");
97
98 CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", skinPostProc, false, false);
101
102 iNit = true;
103 }
104 int step;
105 CHKERR TSGetTimeStepNumber(ts, &step);
106
107 if (pRT && (step) % pRT == 0) {
108 // CHKERR mField.loop_finite_elements("DMTHERMAL","THERMAL_FE",postProc);
109 // std::ostringstream sss;
110 // sss << "out_thermal_" << step << ".h5m";
111 // CHKERR postProc.writeFile(sss.str().c_str());
112 if (saveSkin) {
113 CHKERR mField.loop_finite_elements("DMTHERMAL", "POST_PROC_SKIN",
115 std::ostringstream sss;
116 sss << "out_skin_" << step << ".h5m";
117 CHKERR skinPostProc.writeFile(sss.str().c_str());
118 }
119 }
121 }
122};
123
124int main(int argc, char *argv[]) {
125
126 const string default_options = "-ksp_type fgmres \n"
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
141 string param_file = "param_file.petsc";
142 if (!static_cast<bool>(ifstream(param_file))) {
143 std::ofstream file(param_file.c_str(), std::ios::ate);
144 if (file.is_open()) {
145 file << default_options;
146 file.close();
147 }
148 }
149
150 MoFEM::Core::Initialize(&argc, &argv, param_file.c_str(), help);
151
152 auto core_log = logging::core::get();
153 core_log->add_sink(
155 LogManager::setLog("THERMALSYNC");
156 MOFEM_LOG_TAG("THERMALSYNC", "thermal");
157
158 try {
159
160 PetscBool flg = PETSC_TRUE;
161 char mesh_file_name[255];
162 CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-my_file",
163 mesh_file_name, 255, &flg);
164 if(flg != PETSC_TRUE) {
165 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_FOUND,
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;
171 CHKERR PetscOptionsGetString(PETSC_NULL,PETSC_NULL,"-my_ground_analysis_data",
172 time_data_file_for_ground_surface,255,&ground_temperature_analysis);
173 if(ground_temperature_analysis) {
174#ifndef WITH_ADOL_C
175 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_INSTALLED,
176 "*** ERROR to do ground thermal analysis MoFEM need to be compiled "
177 "with ADOL-C");
178#endif // WITH_ADOL_C
179 }
180
181 //create MoAB database
182 moab::Core mb_instance;
183 moab::Interface& moab = mb_instance;
184 const char *option;
185 option = "";
186 CHKERR moab.load_file(mesh_file_name, 0, option);
187 //create MoFEM database
188 MoFEM::Core core(moab);
189 MoFEM::Interface& m_field = core;
190
191 DMType dm_name = "DMTHERMAL";
192 CHKERR DMRegister_MoFEM(dm_name);
193 // create dm instance
194 DM dm;
195 CHKERR DMCreate(PETSC_COMM_WORLD, &dm);
196 CHKERR DMSetType(dm, dm_name);
197
198 //set entities bit level (this allow to set refinement levels for h-adaptivity)
199 //only one level is used in this example
200 BitRefLevel bit_level0;
201 bit_level0.set(0);
202 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(0, 3,
203 bit_level0);
204
205 //Fields H1 space rank 1
206 CHKERR m_field.add_field("TEMP", H1, AINSWORTH_LEGENDRE_BASE, 1,
207 MB_TAG_SPARSE, MF_ZERO);
208 CHKERR m_field.add_field("TEMP_RATE", H1, AINSWORTH_LEGENDRE_BASE, 1,
209 MB_TAG_SPARSE, MF_ZERO);
210
211 //Add field H1 space rank 3 to approximate geometry using hierarchical basis
212 //For 10 node tets, before use, geometry is projected on that field (see below)
213 CHKERR m_field.add_field(
214 "MESH_NODE_POSITIONS",H1,AINSWORTH_LEGENDRE_BASE,3,MB_TAG_SPARSE,MF_ZERO
215 );
216
217 //meshset consisting all entities in mesh
218 EntityHandle root_set = moab.get_root_set();
219 //add entities to field (root_mesh, i.e. on all mesh etities fields are approx.)
220 CHKERR m_field.add_ents_to_field_by_type(root_set,MBTET,"TEMP");
221 CHKERR m_field.add_ents_to_field_by_type(root_set,MBTET,"TEMP_RATE");
222
223 int order;
224 CHKERR PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-my_order", &order, &flg);
225
226 if (flg != PETSC_TRUE) {
227 order = 1;
228 }
229 // set app. order
230 // see Hierarchic Finite Element Bases on Unstructured Tetrahedral Meshes
231 // (Mark Ainsworth & Joe Coyle) for simplicity of example to all entities is
232 // applied the same order
233 CHKERR m_field.set_field_order(root_set,MBTET,"TEMP",order);
234 CHKERR m_field.set_field_order(root_set,MBTRI,"TEMP",order);
235 CHKERR m_field.set_field_order(root_set,MBEDGE,"TEMP",order);
236 CHKERR m_field.set_field_order(root_set,MBVERTEX,"TEMP",1);
237
238 CHKERR m_field.set_field_order(root_set,MBTET,"TEMP_RATE",order);
239 CHKERR m_field.set_field_order(root_set,MBTRI,"TEMP_RATE",order);
240 CHKERR m_field.set_field_order(root_set,MBEDGE,"TEMP_RATE",order);
241 CHKERR m_field.set_field_order(root_set,MBVERTEX,"TEMP_RATE",1);
242
243 //geometry approximation is set to 2nd oreder
244 CHKERR m_field.add_ents_to_field_by_type(root_set, MBTET,
245 "MESH_NODE_POSITIONS");
246 CHKERR m_field.set_field_order(0,MBTET,"MESH_NODE_POSITIONS",2);
247 CHKERR m_field.set_field_order(0,MBTRI,"MESH_NODE_POSITIONS",2);
248 CHKERR m_field.set_field_order(0,MBEDGE,"MESH_NODE_POSITIONS",2);
249 CHKERR m_field.set_field_order(0,MBVERTEX,"MESH_NODE_POSITIONS",1);
250
251 // configure blocks by parsing config file
252 // it allow to set approximation order for each block independently
253 PetscBool block_config;
254 char block_config_file[255];
255 CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-my_block_config",
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 // std::cerr << block_config_file << std::endl;
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);
312 Range block_ents;
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);
322 CHKERR m_field.set_field_order(ents_to_set_order, "TEMP",
323 block_data[it->getMeshsetId()].oRder);
324 CHKERR m_field.set_field_order(ents_to_set_order, "TEMP_RATE",
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;
340 SETERRQ(PETSC_COMM_SELF,MOFEM_STD_EXCEPTION_THROW,ss.str().c_str());
341 }
342 }
343
344 // this default class to calculate thermal elements
345 ThermalElement thermal_elements(m_field);
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 // add rate of temperature to data field of finite element
351 CHKERR m_field.modify_finite_element_add_field_data("THERMAL_FE",
352 "TEMP_RATE");
353 // and temperature element default element operators at integration (gauss)
354 // points
355 CHKERR addHOOpsVol("MESH_NODE_POSITIONS", thermal_elements.getLoopFeRhs(), true,
356 false, false, false);
357 CHKERR addHOOpsVol("MESH_NODE_POSITIONS", thermal_elements.getLoopFeLhs(), true,
358 false, false, false);
359 CHKERR thermal_elements.setTimeSteppingProblem("TEMP", "TEMP_RATE");
360
361 // set block material data from option file
362 std::map<int, ThermalElement::BlockData>::iterator mit;
363 mit = thermal_elements.setOfBlocks.begin();
364 for (; mit != thermal_elements.setOfBlocks.end(); mit++) {
365 // std::cerr << mit->first << std::endl;
366 // std::cerr << block_data[mit->first].cOnductivity << " " <<
367 // block_data[mit->first].cApacity << std::endl;
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
384 GroundSurfaceTemperature ground_surface(m_field);
385 CrudeClimateModel time_data(time_data_file_for_ground_surface);
386 GroundSurfaceTemperature::PreProcess exectuteGenericClimateModel(&time_data);
387 if (ground_temperature_analysis) {
388 CHKERR ground_surface.addSurfaces("TEMP");
389 CHKERR ground_surface.setOperators(&time_data, "TEMP");
390 }
391#endif //__GROUND_SURFACE_TEMPERATURE_HPP
392
393 //build database, i.e. declare dofs, elements and adjacencies
394
395 // build field
396 CHKERR m_field.build_fields();
397 // project 10 node tet approximation of geometry on hierarchical basis
398 Projection10NodeCoordsOnField ent_method_material(m_field,
399 "MESH_NODE_POSITIONS");
400 CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
401
402 // set initial temperature from Cubit blocksets
403 mit = thermal_elements.setOfBlocks.begin();
404 for (; mit != thermal_elements.setOfBlocks.end(); mit++) {
405 if (mit->second.initTemp != 0) {
406 Range vertices;
407 CHKERR moab.get_connectivity(mit->second.tEts, vertices, true);
408 CHKERR m_field.getInterface<FieldBlas>()->setField(
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);
422 CHKERR m_field.getInterface<FieldBlas>()->setField(data[0], MBVERTEX,
423 block_verts, "TEMP");
424 }
425 }
426
428 if (block_data[it->getMeshsetId()].initTemp != 0) {
429 Range block_ents;
430 CHKERR moab.get_entities_by_handle(it->meshset, block_ents, true);
431 Range vertices;
432 CHKERR moab.get_connectivity(block_ents, vertices, true);
433 CHKERR m_field.getInterface<FieldBlas>()->setField(
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;
445 CHKERR PetscOptionsGetBool(PETSC_NULL, PETSC_NULL, "-dm_is_partitioned",
446 &is_partitioned, PETSC_NULL);
447
448 Range thermal_element_ents;
449 CHKERR m_field.get_finite_element_entities_by_dimension("THERMAL_FE", 3,
450 thermal_element_ents);
451
452 PetscBool save_skin = PETSC_TRUE;
453 CHKERR PetscOptionsGetBool(PETSC_NULL, PETSC_NULL, "-my_save_skin",
454 &save_skin, PETSC_NULL);
455
456 Skinner skin(&m_field.get_moab());
457 Range skin_faces; // skin faces from 3d ents
458 CHKERR skin.find_skin(0, thermal_element_ents, false, skin_faces);
459 Range proc_skin;
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) {
469 CHKERR m_field.add_finite_element("POST_PROC_SKIN");
470 CHKERR m_field.modify_finite_element_add_field_row("POST_PROC_SKIN", "TEMP");
471 CHKERR m_field.modify_finite_element_add_field_col("POST_PROC_SKIN", "TEMP");
472 CHKERR m_field.modify_finite_element_add_field_data("POST_PROC_SKIN", "TEMP");
473 CHKERR m_field.modify_finite_element_add_field_data("POST_PROC_SKIN",
474 "MESH_NODE_POSITIONS");
475 CHKERR m_field.add_ents_to_finite_element_by_dim(proc_skin, 2,
476 "POST_PROC_SKIN");
477 }
478
479 // build finite elemnts
481 // build adjacencies
482 CHKERR m_field.build_adjacencies(bit_level0);
483
484 // delete old temperature recorded series
485 SeriesRecorder *recorder_ptr;
486 CHKERR m_field.getInterface(recorder_ptr);
487 if (recorder_ptr->check_series("THEMP_SERIES")) {
488 /*for(_IT_SERIES_STEPS_BY_NAME_FOR_LOOP_(recorder_ptr,"THEMP_SERIES",sit)) {
489 CHKERR
490 recorder_ptr->load_series_data("THEMP_SERIES",sit->get_step_number());
491 }*/
492 CHKERR recorder_ptr->delete_recorder_series("THEMP_SERIES");
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];
499 CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-my_eval_points_file",
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()) {
504 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Cannot open file %s",
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 // set dm data structure which created mofem data structures
514 CHKERR DMMoFEMCreateMoFEM(dm, &m_field, dm_name, bit_level0);
515 CHKERR DMSetFromOptions(dm);
516 // add elements to dm
517 CHKERR DMMoFEMAddElement(dm, "THERMAL_FE");
518 CHKERR DMMoFEMAddElement(dm, "THERMAL_FLUX_FE");
519 CHKERR DMMoFEMAddElement(dm, "THERMAL_CONVECTION_FE");
520 CHKERR DMMoFEMAddElement(dm, "THERMAL_RADIATION_FE");
521 if (save_skin)
522 CHKERR DMMoFEMAddElement(dm, "POST_PROC_SKIN");
523
524#ifdef __GROUND_SURFACE_TEMPERATURE_HPP
525 if (ground_temperature_analysis) {
526 CHKERR DMMoFEMAddElement(dm, "GROUND_SURFACE_FE");
527 }
528#endif //__GROUND_SURFACE_TEMPERATURE_HPP
529
530 CHKERR DMSetUp(dm);
531
532 // create matrices
533 Vec T, F;
535 CHKERR VecDuplicate(T, &F);
536 Mat A;
538
539 DirichletTemperatureBc dirichlet_bc(m_field, "TEMP", A, T, F);
540 ThermalElement::UpdateAndControl update_velocities(m_field, "TEMP",
541 "TEMP_RATE");
542 ThermalElement::TimeSeriesMonitor monitor(m_field, "THEMP_SERIES", "TEMP",
543 eval_points);
544 MonitorPostProc post_proc(m_field);
545
546 // Initialize data with values save of on the field
547 CHKERR VecZeroEntries(T);
548 CHKERR DMoFEMMeshToLocalVector(dm, T, INSERT_VALUES, SCATTER_FORWARD);
549 CHKERR DMoFEMPreProcessFiniteElements(dm, &dirichlet_bc);
550 CHKERR DMoFEMMeshToGlobalVector(dm, T, INSERT_VALUES, SCATTER_REVERSE);
551
552 // preprocess
553 CHKERR DMMoFEMTSSetIFunction(dm, DM_NO_ELEMENT, NULL, &update_velocities,
554 NULL);
555 CHKERR DMMoFEMTSSetIFunction(dm, DM_NO_ELEMENT, NULL, &dirichlet_bc, NULL);
556 CHKERR DMMoFEMTSSetIJacobian(dm, DM_NO_ELEMENT, NULL, &dirichlet_bc, NULL);
557#ifdef __GROUND_SURFACE_TEMPERATURE_HPP
559 &exectuteGenericClimateModel, NULL);
560 { // add preprocessor, calculating angle on which sun ray on the surface
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++) {
568 CHKERR DMMoFEMTSSetIFunction(dm, DM_NO_ELEMENT, NULL, &*it, NULL);
569 }
570 }
571 }
572#endif //__GROUND_SURFACE_TEMPERATURE_HPP
573
574 // loops rhs
575 CHKERR DMMoFEMTSSetIFunction(dm, "THERMAL_FE", &thermal_elements.feRhs, NULL,
576 NULL);
577 CHKERR DMMoFEMTSSetIFunction(dm, "THERMAL_FLUX_FE", &thermal_elements.feFlux,
578 NULL, NULL);
579 CHKERR DMMoFEMTSSetIFunction(dm, "THERMAL_CONVECTION_FE",
580 &thermal_elements.feConvectionRhs, NULL, NULL);
581 CHKERR DMMoFEMTSSetIFunction(dm, "THERMAL_RADIATION_FE",
582 &thermal_elements.feRadiationRhs, NULL, NULL);
583#ifdef __GROUND_SURFACE_TEMPERATURE_HPP
584 if (ground_temperature_analysis) {
585 CHKERR DMMoFEMTSSetIFunction(dm, "GROUND_SURFACE_FE",
586 &ground_surface.getFeGroundSurfaceRhs(), NULL,
587 NULL);
588 }
589#endif //__GROUND_SURFACE_TEMPERATURE_HPP
590
591 // loops lhs
592 CHKERR DMMoFEMTSSetIJacobian(dm, "THERMAL_FE", &thermal_elements.feLhs, NULL,
593 NULL);
594 CHKERR DMMoFEMTSSetIJacobian(dm, "THERMAL_CONVECTION_FE",
595 &thermal_elements.feConvectionLhs, NULL, NULL);
596 CHKERR DMMoFEMTSSetIJacobian(dm, "THERMAL_RADIATION_FE",
597 &thermal_elements.feRadiationLhs, NULL, NULL);
598#ifdef __GROUND_SURFACE_TEMPERATURE_HPP
599 if (ground_temperature_analysis) {
600 CHKERR DMMoFEMTSSetIJacobian(dm, "GROUND_SURFACE_FE",
601 &ground_surface.getFeGroundSurfaceLhs(), NULL,
602 NULL);
603 }
604#endif //__GROUND_SURFACE_TEMPERATURE_HPP
605
606 //postprocess
607 CHKERR DMMoFEMTSSetIFunction(dm,DM_NO_ELEMENT,NULL,NULL,&dirichlet_bc);
608 CHKERR DMMoFEMTSSetIJacobian(dm,DM_NO_ELEMENT,NULL,NULL,&dirichlet_bc);
609
610 TsCtx *ts_ctx;
612 //add monitor operator
613 ts_ctx->getPostProcessMonitor().push_back(&monitor);
614 ts_ctx->getPostProcessMonitor().push_back(&post_proc);
615
616 //create time solver
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 //add monitor to TS solver
624 CHKERR TSMonitorSet(ts,TsMonitorSet,ts_ctx,PETSC_NULL); // !!!
625
626 CHKERR recorder_ptr->add_series_recorder("THEMP_SERIES");
627 //start to record
628 CHKERR recorder_ptr->initialize_series_recorder("THEMP_SERIES");
629
630 double ftime = 1;
631 CHKERR TSSetDuration(ts,PETSC_DEFAULT,ftime);
632 CHKERR TSSetFromOptions(ts);
633 CHKERR TSSetDM(ts,dm);
634
635 CHKERR TSSolve(ts,T);
636 CHKERR TSGetTime(ts,&ftime);
637
638 //end recoder
639 CHKERR recorder_ptr->finalize_series_recorder("THEMP_SERIES");
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 // save solution, if boundary conditions are defined you can use that file in
654 // mechanical problem to calculate thermal stresses
655 PetscBool save_solution = PETSC_TRUE;
656 CHKERR PetscOptionsGetBool(PETSC_NULL, PETSC_NULL, "-my_save_solution",
657 &save_solution, PETSC_NULL);
658
659 if (save_solution) {
660 if (is_partitioned) {
661 CHKERR moab.write_file("solution.h5m");
662 } else {
663 if (m_field.get_comm_rank() == 0) {
664 CHKERR moab.write_file("solution.h5m");
665 }
666 }
667 }
668
669 CHKERR MatDestroy(&A);
670 CHKERR VecDestroy(&T);
671 CHKERR VecDestroy(&F);
672
673 CHKERR TSDestroy(&ts);
674
675 }
677
678 return MoFEM::Core::Finalize();
679}
#define DM_NO_ELEMENT
Definition: DMMoFEM.hpp:10
const std::string default_options
std::string param_file
int main()
Definition: adol-c_atom.cpp:46
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ MF_ZERO
Definition: definitions.h:98
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ H1
continuous field
Definition: definitions.h:85
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ BLOCKSET
Definition: definitions.h:148
@ MOFEM_NOT_FOUND
Definition: definitions.h:33
@ MOFEM_NOT_INSTALLED
Definition: definitions.h:37
@ MOFEM_STD_EXCEPTION_THROW
Definition: definitions.h:39
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
constexpr int order
@ F
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:483
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
Definition: DMMoFEM.cpp:786
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.
Definition: DMMoFEM.cpp:118
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:509
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMoFEM.cpp:1183
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
Definition: DMMoFEM.cpp:1128
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
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
Definition: DMMoFEM.cpp:839
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMoFEM.cpp:1153
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMoFEM.cpp:521
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:532
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.
MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2 or H1 field gradient.
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
MoFEMErrorCode writeFile(const std::string file_name, const char *file_type="MOAB", const char *file_options="PARALLEL=WRITE_PART")
wrote results in (MOAB) format, use "file_name.h5m"
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
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.
virtual MoFEMErrorCode loop_finite_elements(const std::string problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
#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)
MoFEMErrorCode addThermalElements(const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
add thermal element on tets
MoFEMErrorCode addThermalFluxElement(const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
add heat flux element
MoFEMErrorCode addThermalConvectionElement(const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
add convection element
MoFEMErrorCode addThermalRadiationElement(const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
add Non-linear Radiation element
MoFEMErrorCode setTimeSteppingProblem(string field_name, string rate_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
set up operators for unsteady heat flux; convection; radiation problem
MoFEM::TsCtx * ts_ctx
Definition: level_set.cpp:1932
char mesh_file_name[255]
const double T
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
Definition: TsCtx.cpp:259
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)
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
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)
constexpr AssemblyType A
PetscErrorCode monitor(TS ts, PetscInt step, PetscReal crtime, Vec u, void *ctx)
Implementation of ground surface temperature.
MoFEMErrorCode setOperators(GenericClimateModel *time_data_ptr, string field_name, const string mesh_nodals_positions="MESH_NODE_POSITIONS")
boost::ptr_vector< SolarRadiationPreProcessor > preProcessShade
MoFEMErrorCode addSurfaces(const string field_name, const string mesh_nodals_positions="MESH_NODE_POSITIONS")
Managing BitRefLevels.
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
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Deprecated interface functions.
structure for User Loop Methods on finite elements
Basic algebra on fields.
Definition: FieldBlas.hpp:21
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:298
static boost::shared_ptr< std::ostream > getStrmSync()
Get the strm sync object.
Definition: LogManager.cpp:348
Projection of edge entities with one mid-node on hierarchical basis.
TS ts
time solver
Interface for Time Stepping (TS) solver.
Definition: TsCtx.hpp:17
BasicMethodsSequence & getPostProcessMonitor()
Get the postProcess to do Monitor object.
Definition: TsCtx.hpp:144
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
MoFEMErrorCode postProcess()
function is run at the end of loop
MoFEMErrorCode operator()()
function is run for every finite element
MoFEMErrorCode preProcess()
function is run at the beginning of loop
MoFEM::Interface & mField
PetscBool saveSkin
bool iNit
int pRT
PostProcFaceOnRefinedMesh skinPostProc
MonitorPostProc(MoFEM::Interface &m_field)
PostProcVolumeOnRefinedMesh postProc
int * step
Postprocess on face.
MoFEMErrorCode generateReferenceElementMesh()
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
Post processing.
TS monitore it records temperature at time steps.
this calass is to control time stepping
structure grouping operators and data used for thermal problems
MyTriFE feRadiationLhs
MyTriFE feConvectionLhs
MyVolumeFE feRhs
cauclate right hand side for tetrahedral elements
MyVolumeFE & getLoopFeLhs()
get lhs volume element
MyTriFE feConvectionRhs
MyVolumeFE & getLoopFeRhs()
get rhs volume element
MyVolumeFE feLhs
MyTriFE feRadiationRhs
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
static char help[]