v0.14.0
Classes | Typedefs | Functions | Variables
thermal_unsteady.cpp File Reference

Example of thermal unsteady analyze. More...

#include <BasicFiniteElements.hpp>

Go to the source code of this file.

Classes

struct  BlockOptionData
 
struct  MonitorPostProc
 

Typedefs

using PostProcFaceEle = PostProcBrokenMeshInMoab< FaceElementForcesAndSourcesCore >
 

Functions

int main (int argc, char *argv[])
 

Variables

static char help []
 

Detailed Description

Example of thermal unsteady analyze.

TODO:

Todo:
Make it work in distributed meshes with multigird solver. At the moment it is not working efficient as can.

Definition in file thermal_unsteady.cpp.

Typedef Documentation

◆ PostProcFaceEle

Definition at line 29 of file thermal_unsteady.cpp.

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)

Definition at line 124 of file thermal_unsteady.cpp.

124  {
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(
154  LogManager::createSink(LogManager::getStrmSync(), "THERMALSYNC"));
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
480  CHKERR m_field.build_finite_elements();
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<
563  GroundSurfaceTemperature::SolarRadiationPreProcessor>::iterator it,
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;
611  DMMoFEMGetTsCtx(dm,&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  }
676  CATCH_ERRORS;
677 
678  return MoFEM::Core::Finalize();
679 }

Variable Documentation

◆ help

char help[]
static
Initial value:
=
"-my_file mesh file\n"
"-order set approx. order to all blocks\n"
"-my_block_config set block data\n"
"-my_ground_analysis_data data for crude climate model\n"
"\n"

Definition at line 31 of file thermal_unsteady.cpp.

MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
MoFEM::CoreInterface::loop_dofs
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.
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
MOFEM_STD_EXCEPTION_THROW
@ MOFEM_STD_EXCEPTION_THROW
Definition: definitions.h:39
EntityHandle
MoFEM::addHOOpsVol
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
Definition: HODataOperators.hpp:764
MoFEM::SeriesRecorder::check_series
virtual bool check_series(const std::string &name) const
check if series is in database
Definition: SeriesRecorder.cpp:301
MoFEM::CoreInterface::modify_finite_element_add_field_row
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
MoFEM::SeriesRecorder
Definition: SeriesRecorder.hpp:25
MoFEM::TsCtx::getPostProcessMonitor
BasicMethodsSequence & getPostProcessMonitor()
Get the postProcess to do Monitor object.
Definition: TsCtx.hpp:144
DirichletTemperatureBc
Definition: DirichletBC.hpp:239
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
MoFEM::DMoFEMMeshToLocalVector
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:523
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
ts_ctx
MoFEM::TsCtx * ts_ctx
Definition: level_set.cpp:1932
MoFEM::CoreInterface::add_ents_to_field_by_type
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.
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:497
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
MoFEM::TsMonitorSet
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
Definition: TsCtx.cpp:259
MoFEM::DMCreateGlobalVector_MoFEM
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMoFEM.cpp:1167
MoFEM::SeriesRecorder::add_series_recorder
virtual MoFEMErrorCode add_series_recorder(const std::string &series_name)
Definition: SeriesRecorder.cpp:90
MoFEM::DMCreateMatrix_MoFEM
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMoFEM.cpp:1197
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::CoreInterface::add_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
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::DMMoFEMTSSetIJacobian
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:853
MoFEM::CoreInterface::modify_finite_element_add_field_col
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
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
help
static char help[]
Definition: thermal_unsteady.cpp:31
MoFEM::CoreInterface::add_field
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.
DM_NO_ELEMENT
#define DM_NO_ELEMENT
Definition: DMMoFEM.hpp:10
MoFEM::SeriesRecorder::finalize_series_recorder
virtual MoFEMErrorCode finalize_series_recorder(const std::string &serie_name)
Definition: SeriesRecorder.cpp:264
MoFEM::DMoFEMMeshToGlobalVector
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMoFEM.cpp:535
MoFEM::TsCtx
Interface for Time Stepping (TS) solver.
Definition: TsCtx.hpp:17
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
MOFEM_NOT_INSTALLED
@ MOFEM_NOT_INSTALLED
Definition: definitions.h:37
MoFEM::DMoFEMPreProcessFiniteElements
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:546
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
MoFEM::CoreInterface::add_ents_to_finite_element_by_dim
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
MoFEM::DMMoFEMCreateMoFEM
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:114
MOFEM_NOT_FOUND
@ MOFEM_NOT_FOUND
Definition: definitions.h:33
ThermalElement::TimeSeriesMonitor
TS monitore it records temperature at time steps.
Definition: ThermalElement.hpp:579
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:23
MoFEM::CoreInterface::get_finite_element_entities_by_dimension
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
MoFEM::SeriesRecorder::delete_recorder_series
virtual MoFEMErrorCode delete_recorder_series(const std::string &series_name)
Definition: SeriesRecorder.cpp:111
ThermalElement::UpdateAndControl
this calass is to control time stepping
Definition: ThermalElement.hpp:562
Range
ThermalElement
structure grouping operators and data used for thermal problems
Definition: ThermalElement.hpp:27
MoFEM::CoreTmp< 0 >::Initialize
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
FTensor::dd
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)
Definition: ddTensor0.hpp:33
MF_ZERO
@ MF_ZERO
Definition: definitions.h:98
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#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.
Definition: MeshsetsManager.hpp:71
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1102
BLOCKSET
@ BLOCKSET
Definition: definitions.h:148
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MoFEM::DMMoFEMTSSetIFunction
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:800
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition: DeprecatedPetsc.hpp:172
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
MoFEM::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
MoFEM::DMMoFEMGetTsCtx
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
Definition: DMMoFEM.cpp:1142
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MoFEM::CoreInterface::build_adjacencies
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
MoFEM::SeriesRecorder::initialize_series_recorder
virtual MoFEMErrorCode initialize_series_recorder(const std::string &serie_name)
Definition: SeriesRecorder.cpp:246
MoFEM::CoreInterface::set_field_order
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.
MonitorPostProc
Definition: nonlinear_dynamics.cpp:30
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEM::FieldBlas
Basic algebra on fields.
Definition: FieldBlas.hpp:21
F
@ F
Definition: free_surface.cpp:394
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182