v0.9.0
GroundSurfaceTemerature.hpp
Go to the documentation of this file.
1 /** \file GroundSurfaceTemperature.hpp
2  * \brief Operators and data structures for thermal analys
3  *
4  * Implementation of boundary conditions for ground temerature
5  *
6  */
7 
8 /*
9  * This file is part of MoFEM.
10  * MoFEM is free software: you can redistribute it and/or modify it under
11  * the terms of the GNU Lesser General Public License as published by the
12  * Free Software Foundation, either version 3 of the License, or (at your
13  * option) any later version.
14  *
15  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
16  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
18  * License for more details.
19  *
20  * You should have received a copy of the GNU Lesser General Public
21  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
22 
23 #ifndef __GROUND_SURFACE_TEMPERATURE_HPP
24 #define __GROUND_SURFACE_TEMPERATURE_HPP
25 
26 #ifndef WITH_ADOL_C
27 #error "MoFEM need to be compiled with ADOL-C"
28 #endif
29 
30 /** \brief Implementation of ground surface temperature
31  *
32  * Ground surface temperature simulation for different land covers
33  * William R. Herb *, Ben Janke, Omid Mohseni, Heinz G. Stefan
34  * Journal of Hydrology (2008) 356, 327– 343
35  */
37 
39 
40  /** \brief common data used by volume elements
41  * \infroup mofem_thermal_elem
42  */
43  struct CommonData {
45  };
47 
48  /** \brief define surface element
49  *
50  * This element is used to integrate heat fluxes; convection and radiation
51  */
54  int getRule(int order) { return order; };
58  inline const MyTriFE *getFriendFaceFE() {
59  return static_cast<MyTriFE *>(ptrFE);
60  }
61  };
62  using FaceElementForcesAndSourcesCore::dataH1;
63  };
64 
65  MyTriFE feGroundSurfaceRhs; //< radiation element
69 
71  : mField(m_field), feGroundSurfaceRhs(m_field),
72  feGroundSurfaceLhs(m_field){};
73 
74  struct Parameters {
75 
76  double alpha; //< Solar albedo
77  double d; //< Constatnt used to clulate albedo for urtain angle
78  double Cfc; //< Surface heat/moisture transfer coefficient for forced
79  // convection
80  double Cnc; //< Coefficient for natural convection
81  double CSh; //< Wind sheltering coefficient
82  double eps; //< Pavement emissivity
83  double rhoCp; //< Density specific heat pavement (J/m3/°C)
84  Range tRis; //< Triangles on which parameters are defined
85 
86  Parameters(Range &tris) : tRis(tris) {}
87  };
88 
89  struct Asphalt : public Parameters {
90  Asphalt(Range &tris) : Parameters(tris) {
91  alpha = 0.12;
92  d = 0.25; // not estimated, some goods given number
93  Cfc = 0.0015;
94  Cnc = 0.0015;
95  CSh = 1.;
96  eps = 0.94;
97  rhoCp = 2.0e06;
98  }
99  };
100 
101  struct Concrete : public Parameters {
102  Concrete(Range &tris) : Parameters(tris) {
103  alpha = 0.20;
104  d = 0.25; // not estimated, some goods given number
105  Cfc = 0.0015;
106  Cnc = 0.0015;
107  CSh = 1.;
108  eps = 0.94;
109  rhoCp = 2.0e06;
110  }
111  };
112 
113  struct BareSoil : public Parameters {
114  BareSoil(Range &tris) : Parameters(tris) {
115  alpha = 0.15;
116  d = 0.25; // not estimated, some goods given number
117  Cfc = 0.003;
118  Cnc = 0.0015;
119  CSh = 1.;
120  eps = 0.94;
121  rhoCp = 2.0e06;
122  }
123  };
124 
126  addSurfaces(const string field_name,
127  const string mesh_nodals_positions = "MESH_NODE_POSITIONS") {
129 
130  CHKERR mField.add_finite_element("GROUND_SURFACE_FE", MF_ZERO);
132  field_name);
134  field_name);
136  field_name);
137  if (mField.check_field(mesh_nodals_positions))
139  mesh_nodals_positions);
140 
142  if (it->getName().compare(0, 7, "ASPHALT") == 0) {
143  Range tris;
144  CHKERR mField.get_moab().get_entities_by_type(it->meshset, MBTRI, tris,
145  true);
146  blockData.push_back(new Asphalt(tris));
148  "GROUND_SURFACE_FE");
149  }
150  }
151 
153  if (it->getName().compare(0, 8, "CONCRETE") == 0) {
154  Range tris;
155  CHKERR mField.get_moab().get_entities_by_type(it->meshset, MBTRI, tris,
156  true);
157  blockData.push_back(new Concrete(tris));
159  "GROUND_SURFACE_FE");
160  }
161  }
162 
164  if (it->getName().compare(0, 8, "BARESOIL") == 0) {
165  Range tris;
166  CHKERR mField.get_moab().get_entities_by_type(it->meshset, MBTRI, tris,
167  true);
168  blockData.push_back(new BareSoil(tris));
170  "GROUND_SURFACE_FE");
171  }
172  }
173 
175  }
176 
177  static double netSolarRadiation(double alpha, double d, double cos_omega,
178  GenericClimateModel *time_data_ptr) {
179  // Parametrizing the Dependence of Surface Albedo on Solar Zenith Angle
180  // Using Atmospheric Radiation Measurement Program Observations F. Yang
181  // Environmental Modeling Center National Centers for Environmental
182  // Prediction Camp Springs, Maryland
183  alpha = alpha * (1 + d / (1 + 2 * d * cos_omega));
184  return (1 - alpha) * time_data_ptr->Rs; // net solar radiation (W/m2)
185  }
186 
187  static double incomingLongWaveRadiation(double eps,
188  GenericClimateModel *time_data_ptr) {
189  const double sigma = 5.67037321e-8;
190  double sigma_eps = eps * sigma;
191  double ea = time_data_ptr->calculateVapourPressure(
192  time_data_ptr->calculateVapourPressure(time_data_ptr->Td));
193  ea *= 1e-3; // this need to be expressed in kPa
194  // equation taken from William R. Herb but need too look as well at, Klok
195  // and Oerlemans, 2002 incoming longwave radiation (W/m2)
196  return sigma_eps *
197  (time_data_ptr->CR +
198  0.67 * (1 - time_data_ptr->CR) * pow(ea, 0.08)) *
199  pow((time_data_ptr->Ta + 273.15), 4);
200  }
201 
202  static double outgoingLongWaveRadiation(double eps, double T) {
203  const double sigma = 5.67037321e-8;
204  double sigma_eps = eps * sigma;
205  return -sigma_eps * pow(T + 273.15, 4);
206  }
207 
208  static double outgoingLongWaveRadiation_dT(double eps, double T) {
209  const double sigma = 5.67037321e-8;
210  double sigma_eps = eps * sigma;
211  return -4 * sigma_eps * pow(T + 273.15, 3);
212  }
213 
214  /**
215  * \brief convection between ladn or water and atmosphere
216  *
217  * \param T land tamerature
218  * \param u10 wind at 10 m
219  * \param CSh wind sheltering coefficient
220  * \param rhoCp density specific heat pavement (J/m3/°C)
221  * \param Cfc surface heat/moisture transfer coefficient for forced convection
222  * \param Cnc coefficient for natural convection
223  * \param time_data_ptr
224  */
225  static double
226  convectinBetweenLandOrWaterAndAtmosphere(double T, double CSh, double rhoCp,
227  double Cfc, double Cnc,
228  GenericClimateModel *time_data_ptr) {
229  double us =
230  CSh * time_data_ptr->u10; // win speed with sheltering coeeficient
231  double h_conv1 = -rhoCp * Cfc * us * (T - time_data_ptr->Ta);
232  double Tv = time_data_ptr->calculateAbsoluteVirtualTempertaure(
233  T, time_data_ptr->Td, time_data_ptr->P);
234  double Tv_a = time_data_ptr->calculateAbsoluteVirtualTempertaure(
235  time_data_ptr->Ta, time_data_ptr->Td, time_data_ptr->P);
236  double delta_phi = Tv - Tv_a;
237  if (fabs(delta_phi) > 1) {
238  double A = pow(fabs(delta_phi), 0.33);
239  double h_conv2 = -rhoCp * Cnc * A * (T - time_data_ptr->Ta);
240  return h_conv1 + h_conv2;
241  } else {
242  double A = 0.835 + 0.165 * pow(delta_phi, 2);
243  double h_conv2 = -rhoCp * Cnc * A * (T - time_data_ptr->Ta);
244  return h_conv1 + h_conv2;
245  }
246  return h_conv1;
247  }
248 
249  /**
250  * \brief convection between ladn or water and atmosphere
251  *
252  * \param T land tamerature
253  * \param CSh wind sheltering coefficient
254  * \param rhoCp density specific heat pavement (J/m3/°C)
255  * \param Cfc surface heat/moisture transfer coefficient for forced convection
256  * \param Cnc coefficient for natural convection
257  * \param time_data_ptr
258  */
260  double T, double CSh, double rhoCp, double Cfc, double Cnc,
261  GenericClimateModel *time_data_ptr) {
262  double us =
263  CSh * time_data_ptr->u10; // win speed with sheltering coeeficient
264  double h_conv1_dT = -rhoCp * Cfc * us;
265  double Tv = time_data_ptr->calculateAbsoluteVirtualTempertaure(
266  T, time_data_ptr->Td, time_data_ptr->P);
267  double Tv_a = time_data_ptr->calculateAbsoluteVirtualTempertaure(
268  time_data_ptr->Ta, time_data_ptr->Td, time_data_ptr->P);
269  double delta_phi = Tv - Tv_a;
270  if (fabs(delta_phi) > 1) {
271  double A = pow(fabs(delta_phi), 0.33);
272  double Tv_dT = time_data_ptr->calculateAbsoluteVirtualTempertaure_dT(
273  T, time_data_ptr->Td, time_data_ptr->P);
274  double A_dT =
275  ::copysign(1, delta_phi) * Tv_dT * 0.33 / pow(fabs(delta_phi), 0.67);
276  double h_conv2_dT =
277  -(rhoCp * Cnc * A_dT * (T - time_data_ptr->Ta) + rhoCp * Cnc * A);
278  return h_conv1_dT + h_conv2_dT;
279  } else {
280  double A = 0.835 + 0.165 * pow(delta_phi, 2);
281  double Tv_dT = time_data_ptr->calculateAbsoluteVirtualTempertaure_dT(
282  T, time_data_ptr->Td, time_data_ptr->P);
283  double A_dT = 2 * 0.165 * delta_phi * Tv_dT;
284  double h_conv2_dT =
285  -(rhoCp * Cnc * A_dT * (T - time_data_ptr->Ta) + rhoCp * Cnc * A);
286  return h_conv1_dT + h_conv2_dT;
287  }
288  return h_conv1_dT;
289  }
290 
291  struct PreProcess : public MoFEM::FEMethod {
292 
295  : timeDataPtr(time_data_ptr){};
296 
301  }
302  };
303 
305 
309  AdaptiveKDTree kdTree;
310  double ePs;
311  bool iNit;
312 
314  GenericClimateModel *time_data_ptr,
315  Parameters *parameters_ptr, double eps = 1e-6)
316  : mField(m_field), timeDataPtr(time_data_ptr),
317  pArametersPtr(parameters_ptr), kdTree(&m_field.get_moab()), ePs(eps),
318  iNit(false) {
319  azimuth = zenith = 0;
320  };
322  if (kdTree_rootMeshset) {
323  mField.get_moab().delete_entities(&kdTree_rootMeshset, 1);
324  }
325  }
326 
327  Range sKin, skinNodes;
329 
330  MoFEMErrorCode getSkin(Range &tets) {
332  Skinner skin(&mField.get_moab());
333  CHKERR skin.find_skin(0, tets, false, sKin);
334  CHKERR mField.get_moab().create_meshset(MESHSET_SET, kdTree_rootMeshset);
335  CHKERR kdTree.build_tree(sKin, &kdTree_rootMeshset);
337  }
338 
339  double azimuth, zenith;
340 
343 
344  int def_VAL = 0;
345  Tag th_solar_exposure;
346  CHKERR mField.get_moab().tag_get_handle(
347  "SOLAR_EXPOSURE", 1, MB_TYPE_INTEGER, th_solar_exposure,
348  MB_TAG_CREAT | MB_TAG_SPARSE, &def_VAL);
349 
350  double zero[3] = {0, 0, 0};
351  Tag th_solar_radiation;
352  CHKERR mField.get_moab().tag_get_handle(
353  "SOLAR_RADIATION", 1, MB_TYPE_DOUBLE, th_solar_radiation,
354  MB_TAG_CREAT | MB_TAG_SPARSE, zero);
355 
356  Tag th_ray_direction;
357  CHKERR mField.get_moab().tag_get_handle(
358  "SUN_RAY", 3, MB_TYPE_DOUBLE, th_ray_direction,
359  MB_TAG_CREAT | MB_TAG_SPARSE, zero);
360 
361  if (!iNit)
362  CHKERR mField.get_moab().get_connectivity(pArametersPtr->tRis,
363  skinNodes, true);
364 
365  if (iNit) {
368  }
369  }
370  iNit = true;
371 
374  // assume that X pointing to North
375  double ray_unit_dir[] = {
376  cos(azimuth * M_PI / 180) * sin(zenith * M_PI / 180),
377  sin(azimuth * M_PI / 180) * sin(zenith * M_PI / 180),
378  cos(zenith * M_PI / 180)};
379 
380  vector<EntityHandle> triangles_out;
381  vector<double> distance_out;
382 
383  double diffN[6];
384  CHKERR ShapeDiffMBTRI(diffN);
385  Range::iterator tit = pArametersPtr->tRis.begin();
386  for (; tit != pArametersPtr->tRis.end(); tit++) {
387 
388  if (ray_unit_dir[2] <= 0) {
389  CHKERR mField.get_moab().tag_set_data(th_solar_radiation, &*tit, 1,
390  zero);
391  continue;
392  }
393 
394  int num_nodes;
395  const EntityHandle *conn;
396  CHKERR mField.get_moab().get_connectivity(*tit, conn, num_nodes, true);
397  double coords[9];
398  CHKERR mField.get_moab().get_coords(conn, 3, coords);
399 
400  double normal[3];
401  CHKERR ShapeFaceNormalMBTRI(diffN, coords, normal);
402 
403  for (int nn = 1; nn < 3; nn++) {
404  for (int dd = 0; dd < 3; dd++) {
405  coords[dd] += coords[3 * nn + dd];
406  }
407  }
408  for (int dd = 0; dd < 3; dd++) {
409  coords[dd] /= 3;
410  }
411 
412  triangles_out.resize(0);
413  distance_out.resize(0);
414  CHKERR kdTree.ray_intersect_triangles(kdTree_rootMeshset, 1e-12,
415  ray_unit_dir, coords,
416  triangles_out, distance_out);
417 
418  double exposed = 0;
419  if (triangles_out.size() > 0) {
420  for (unsigned int nn = 0; nn < triangles_out.size(); nn++) {
421  if (exposed < distance_out[nn])
422  exposed = distance_out[nn];
423  }
424  }
425 
426  double hsol;
427  if (exposed > ePs) {
428  hsol = 0;
429  } else {
430  double cos_phi = 0;
431  for (int nn = 0; nn < 3; nn++) {
432  cos_phi += normal[nn] * ray_unit_dir[nn];
433  }
434  cos_phi /=
435  sqrt(pow(normal[0], 2) + pow(normal[1], 2) + pow(normal[2], 2));
436  cos_phi = fabs(cos_phi);
437 
439  cos_phi, timeDataPtr);
440  }
441 
442  CHKERR mField.get_moab().tag_set_data(th_solar_radiation, &*tit, 1,
443  &hsol);
444  }
445 
446  Range::iterator nit = skinNodes.begin();
447  for (; nit != skinNodes.end(); nit++) {
448 
449  if (ray_unit_dir[2] <= 0) {
450  int set = 0;
451  CHKERR mField.get_moab().tag_set_data(th_solar_exposure, &*nit, 1,
452  &set);
453  CHKERR mField.get_moab().tag_set_data(th_ray_direction, &*nit, 1,
454  zero);
455  continue;
456  }
457 
458  double coords[3];
459  CHKERR mField.get_moab().get_coords(&*nit, 1, coords);
460 
461  triangles_out.resize(0);
462  distance_out.resize(0);
463  CHKERR kdTree.ray_intersect_triangles(kdTree_rootMeshset, 1e-12,
464  ray_unit_dir, coords,
465  triangles_out, distance_out);
466 
467  double exposed = 0;
468  if (triangles_out.size() > 0) {
469  for (unsigned int nn = 0; nn < triangles_out.size(); nn++) {
470  if (exposed < distance_out[nn])
471  exposed = distance_out[nn];
472  }
473  }
474 
475  int set = 1;
476  if (exposed > ePs) {
477  set = 0;
478  }
479  CHKERR mField.get_moab().tag_set_data(th_solar_exposure, &*nit, 1,
480  &set);
481  CHKERR mField.get_moab().tag_set_data(th_ray_direction, &*nit, 1,
482  ray_unit_dir);
483  }
484 
486  }
487  };
488 
489  boost::ptr_vector<Parameters> blockData;
490  boost::ptr_vector<SolarRadiationPreProcessor> preProcessShade;
491 
492  /** \brief opearator to caulate temperature at Gauss points
493  * \infroup mofem_thermal_elem
494  */
497 
499  OpGetTriTemperatureAtGaussPts(const string field_name,
500  VectorDouble &field_at_gauss_pts)
502  field_name, ForcesAndSourcesCore::UserDataOperator::OPROW),
503  fieldAtGaussPts(field_at_gauss_pts) {}
504 
505  /** \brief operator calculating temperature and rate of temperature
506  *
507  * temperature temperature or rate of temperature is calculated
508  * multiplyingshape functions by degrees of freedom
509  */
510  MoFEMErrorCode doWork(int side, EntityType type,
513 
514  if (data.getFieldData().size() == 0)
516  int nb_dofs = data.getFieldData().size();
517  int nb_gauss_pts = data.getN().size1();
518 
519  // initialize
520  fieldAtGaussPts.resize(nb_gauss_pts);
521  if (type == MBVERTEX) {
522  // loop over shape functions on entities allways start from
523  // vertices, so if nodal shape functions are processed, vector of
524  // field values is zeroad at initialization
525  fill(fieldAtGaussPts.begin(), fieldAtGaussPts.end(), 0);
526  }
527 
528  for (int gg = 0; gg < nb_gauss_pts; gg++) {
529  fieldAtGaussPts[gg] +=
530  inner_prod(data.getN(gg, nb_dofs), data.getFieldData());
531  }
532 
534  }
535  };
536 
538 
543 
544  OpRhs(const string field_name, GenericClimateModel *time_data_ptr,
545  Parameters *parameters_ptr, CommonData &common_data,
546  bool _ho_geometry = false)
548  field_name, ForcesAndSourcesCore::UserDataOperator::OPROW),
549  commonData(common_data), timeDataPtr(time_data_ptr),
550  pArametersPtr(parameters_ptr), ho_geometry(_ho_geometry) {}
551 
555 
556  EntityHandle element_ent = getNumeredEntFiniteElementPtr()->getEnt();
557  const EntityHandle *conn;
558  int num_nodes;
559  CHKERR getFaceFE()->mField.get_moab().get_connectivity(element_ent, conn,
560  num_nodes, true);
561  int def_VAL[] = {0};
562  Tag th_solar_exposure;
563  CHKERR getFaceFE()->mField.get_moab().tag_get_handle(
564  "SOLAR_EXPOSURE", 1, MB_TYPE_INTEGER, th_solar_exposure,
565  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
566  CHKERR getFaceFE()->mField.get_moab().tag_get_data(
567  th_solar_exposure, conn, num_nodes, nodalExposure);
569  }
570 
572  MoFEMErrorCode doWork(int side, EntityType type,
575 
576  if (data.getIndices().size() == 0)
578  if (pArametersPtr->tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
579  pArametersPtr->tRis.end())
581 
582  if (type == MBVERTEX)
584 
585  const FENumeredDofEntity *dof_ptr;
586  CHKERR getNumeredEntFiniteElementPtr()->getRowDofsByPetscGlobalDofIdx(
587  data.getIndices()[0], &dof_ptr);
588  int rank = dof_ptr->getNbOfCoeffs();
589  int nb_row_dofs = data.getIndices().size() / rank;
590 
591  Nf.resize(data.getIndices().size());
592  Nf.clear();
593 
594  for (unsigned int gg = 0; gg < data.getN().size1(); gg++) {
595 
596  double val = getGaussPts()(2, gg);
597 
598  VectorDouble normal;
599  if (ho_geometry) {
600  normal = getNormalsAtGaussPts(gg);
601  } else {
602  normal = getNormal();
603  }
604  val *= norm_2(normal) * 0.5;
605 
606  double T = commonData.temperatureAtGaussPts[gg];
607 
608  double azimuth = timeDataPtr->azimuth;
609  double zenith = timeDataPtr->zenith;
610  // assume that X pointing to North
611  double ray_unit_dir[] = {
612  cos(azimuth * M_PI / 180) * sin(zenith * M_PI / 180),
613  sin(azimuth * M_PI / 180) * sin(zenith * M_PI / 180),
614  cos(zenith * M_PI / 180)};
615  double cos_phi = 0;
616  for (int nn = 0; nn < 3; nn++) {
617  cos_phi += normal[nn] * ray_unit_dir[nn];
618  }
619  cos_phi /= norm_2(normal);
620 
621  eXposure = 0;
622  for (int nn = 0; nn < 3; nn++) {
623  eXposure +=
624  getFriendFaceFE()->dataH1.dataOnEntities[MBVERTEX][0].getN()(gg,
625  nn) *
626  nodalExposure[nn];
627  }
628 
629  double hnet = 0;
630 
631  if (eXposure > 0) {
633  cos_phi, timeDataPtr);
634  }
637 
641  hnet /= (double)86400; // number of second in the day
642 
643  ublas::noalias(Nf) -= val * hnet * data.getN(gg, nb_row_dofs);
644  }
645 
646  CHKERR VecSetValues(getFEMethod()->ts_F, data.getIndices().size(),
647  &data.getIndices()[0], &Nf[0], ADD_VALUES);
648 
650  }
651  };
652 
653  struct OpLhs : public OpRhs {
654 
655  OpLhs(const string field_name, GenericClimateModel *time_data_ptr,
656  Parameters *parameters_ptr, CommonData &common_data,
657  bool _ho_geometry = false)
658  : OpRhs(field_name, time_data_ptr, parameters_ptr, common_data,
659  _ho_geometry) {
660  setOpType(ForcesAndSourcesCore::UserDataOperator::OPROWCOL);
661  }
662 
663  MoFEMErrorCode doWork(int side, EntityType type,
667  }
668 
670 
671  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
672  EntityType col_type,
676 
677  if (row_data.getIndices().size() == 0)
679  if (col_data.getIndices().size() == 0)
681 
682  int nb_row = row_data.getN().size2();
683  int nb_col = col_data.getN().size2();
684 
685  if (row_type == MBVERTEX) {
686  // CHKERR getExposure(); (ierr);
687  }
688 
689  NN.resize(nb_row, nb_col);
690  NN.clear();
691 
692  for (unsigned int gg = 0; gg < row_data.getN().size1(); gg++) {
693 
694  double val = getGaussPts()(2, gg);
695  VectorDouble normal;
696  if (ho_geometry) {
697  normal = getNormalsAtGaussPts(gg);
698  } else {
699  normal = getNormal();
700  }
701  val *= norm_2(normal) * 0.5;
702  double T = commonData.temperatureAtGaussPts[gg];
703 
704  double grad[1] = {0};
709  grad[0] /= (double)86400; // number of second in the day
710 
711  noalias(NN) -=
712  val * grad[0] *
713  outer_prod(row_data.getN(gg, nb_row), col_data.getN(gg, nb_col));
714  }
715 
716  CHKERR MatSetValues((getFEMethod()->ts_B), nb_row,
717  &row_data.getIndices()[0], nb_col,
718  &col_data.getIndices()[0], &NN(0, 0), ADD_VALUES);
719  if (row_side != col_side || row_type != col_type) {
720  transNN.resize(nb_col, nb_row);
721  noalias(transNN) = trans(NN);
723  (getFEMethod()->ts_B), nb_col, &col_data.getIndices()[0], nb_row,
724  &row_data.getIndices()[0], &transNN(0, 0), ADD_VALUES);
725  }
726 
728  }
729  };
730 
732  setOperators(GenericClimateModel *time_data_ptr, string field_name,
733  const string mesh_nodals_positions = "MESH_NODE_POSITIONS") {
735 
736  bool ho_geometry = false;
737  if (mField.check_field(mesh_nodals_positions)) {
738  ho_geometry = true;
739  }
740 
741  {
742  boost::ptr_vector<Parameters>::iterator sit = blockData.begin();
743  for (; sit != blockData.end(); sit++) {
744  // add finite element operator
745  feGroundSurfaceRhs.getOpPtrVector().push_back(
747  field_name, commonData.temperatureAtGaussPts));
748  feGroundSurfaceRhs.getOpPtrVector().push_back(new OpRhs(
749  field_name, time_data_ptr, &*sit, commonData, ho_geometry));
750  preProcessShade.push_back(
751  new SolarRadiationPreProcessor(mField, time_data_ptr, &*sit));
752  }
753  }
754  {
755  boost::ptr_vector<Parameters>::iterator sit = blockData.begin();
756  for (; sit != blockData.end(); sit++) {
757  // add finite element operator
758  feGroundSurfaceLhs.getOpPtrVector().push_back(
760  field_name, commonData.temperatureAtGaussPts));
761  feGroundSurfaceLhs.getOpPtrVector().push_back(new OpLhs(
762  field_name, time_data_ptr, &*sit, commonData, ho_geometry));
763  }
764  }
765 
767  }
768 };
769 
770 #endif //__GROUND_SURFACE_TEMPERATURE_HPP
structure for User Loop Methods on finite elementsIt can be used to calculate stiffness matrices,...
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 add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual moab::Interface & get_moab()=0
static double incomingLongWaveRadiation(double eps, GenericClimateModel *time_data_ptr)
OpLhs(const string field_name, GenericClimateModel *time_data_ptr, Parameters *parameters_ptr, CommonData &common_data, bool _ho_geometry=false)
OpGetTriTemperatureAtGaussPts(const string field_name, VectorDouble &field_at_gauss_pts)
PetscErrorCode ShapeFaceNormalMBTRI(double *diffN, const double *coords, double *normal)
Definition: fem_tools.c:241
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
static double outgoingLongWaveRadiation_dT(double eps, double T)
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
OpRhs(const string field_name, GenericClimateModel *time_data_ptr, Parameters *parameters_ptr, CommonData &common_data, bool _ho_geometry=false)
common data used by volume elements \infroup mofem_thermal_elem
static double netSolarRadiation(double alpha, double d, double cos_omega, GenericClimateModel *time_data_ptr)
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
double calculateVapourPressure(double T)
PetscReal ts_t
time
MoFEMErrorCode preProcess()
function is run at the beginning of loop
static double outgoingLongWaveRadiation(double eps, double T)
PreProcess(GenericClimateModel *time_data_ptr)
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
PetscErrorCode ShapeDiffMBTRI(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:206
opearator to caulate temperature at Gauss points \infroup mofem_thermal_elem
double calculateAbsoluteVirtualTempertaure(double T, double Td, double P)
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
DataForcesAndSourcesCore::EntData EntData
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
static double convectinBetweenLandOrWaterAndAtmosphere(double T, double CSh, double rhoCp, double Cfc, double Cnc, GenericClimateModel *time_data_ptr)
convection between ladn or water and atmosphere
MoFEMErrorCode preProcess()
function is run at the beginning of loop
virtual MoFEMErrorCode set(double t=0)=0
MoFEMErrorCode MatSetValues(Mat M, const DataForcesAndSourcesCore::EntData &row_data, const DataForcesAndSourcesCore::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
Implementation of ground surface temperature.
boost::ptr_vector< SolarRadiationPreProcessor > preProcessShade
double calculateAbsoluteVirtualTempertaure_dT(double T, double Td, double P)
MoFEMErrorCode addSurfaces(const string field_name, const string mesh_nodals_positions="MESH_NODE_POSITIONS")
boost::ptr_vector< Parameters > blockData
MoFEMErrorCode VecSetValues(Vec V, const DataForcesAndSourcesCore::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
GroundSurfaceTemperature(MoFEM::Interface &m_field)
#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.
#define CHKERR
Inline error check.
Definition: definitions.h:596
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
ForcesAndSourcesCore::UserDataOperator UserDataOperator
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
operator calculating temperature and rate of temperature
SolarRadiationPreProcessor(MoFEM::Interface &m_field, GenericClimateModel *time_data_ptr, Parameters *parameters_ptr, double eps=1e-6)
virtual bool check_field(const std::string &name) const =0
check if field is in database
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
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
static const double eps
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:72
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
MoFEMErrorCode setOperators(GenericClimateModel *time_data_ptr, string field_name, const string mesh_nodals_positions="MESH_NODE_POSITIONS")
static double convectinBetweenLandOrWaterAndAtmosphere_dT(double T, double CSh, double rhoCp, double Cfc, double Cnc, GenericClimateModel *time_data_ptr)
convection between ladn or water and atmosphere