v0.13.0
Unsaturated2DFlowUDOP.hpp
Go to the documentation of this file.
1 #ifndef __UFOPERATORS2D_HPP__
2 #define __UFOPERATORS2D_HPP__
3 
4 #include <stdlib.h>
5 #include <BasicFiniteElements.hpp>
6 
7 namespace UFOperators2D {
8 
11 
14 
16 
18 
19 int nth_step = 4;
20 constexpr double natural_bc_values = -1;
21 constexpr double essential_bc_values = 0.0;
22 constexpr double eps = 1e-16;
23 
24 // const int 2 = 3;
26 struct BlockData {
27  int block_id;
28  double K_s; // K_s
29  double ll; // l
30  double theta_r; // theta_r
31  double theta_s; // theta_s
32  double theta_m; // theta_m
33  double nn; // n
34  double alpha; // alpha
35  double h_s; // h_s
36 
37  Range block_ents;
38 
40  : K_s(1.00), ll(0.5000), theta_r(0.0450), theta_s(0.4300),
41  theta_m(0.4300000), nn(5.3800000), alpha(1.812500), h_s(0.0000) {}
42 
45  ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "", "Material options", "none");
46  CHKERR PetscOptionsScalar("-K_s", "K_s", "", K_s, &K_s, PETSC_NULL);
47  CHKERR PetscOptionsScalar("-saturated_conductivity", "K_s", "", K_s, &K_s,
48  PETSC_NULL);
49  CHKERR PetscOptionsScalar("-l", "l", "", ll, &ll, PETSC_NULL);
50  CHKERR PetscOptionsScalar("-theta_r", "theta_r", "", theta_r, &theta_r,
51  PETSC_NULL);
52  CHKERR PetscOptionsScalar("-theta_s", "theta_s", "", theta_s, &theta_s,
53  PETSC_NULL);
54  CHKERR PetscOptionsScalar("-theta_m", "theta_m", "", theta_m, &theta_m,
55  PETSC_NULL);
56  CHKERR PetscOptionsScalar("-n", "n", "", nn, &nn, PETSC_NULL);
57  CHKERR PetscOptionsScalar("-alpha", "alpha", "", alpha, &alpha, PETSC_NULL);
58  CHKERR PetscOptionsScalar("-h_s", "h_s", "", h_s, &h_s, PETSC_NULL);
59  ierr = PetscOptionsEnd();
60  CHKERRG(ierr);
62  };
63 
64  template <typename T> T get_waterContent(T head) {
65  T ret_val;
66  T m = 1 - 1.0 / nn;
67  return ret_val = theta_r +
68  (theta_m - theta_r) / pow(1 + pow(-alpha * head, nn), m);
69  }
70 
71  template <typename T> T get_conductivity(T head, double head_norm) {
72  T ret_val;
73  if (head_norm < h_s) {
74  ret_val = K_s * get_relativeConductivity(head);
75  } else {
76  ret_val = K_s;
77  }
78  return ret_val;
79  }
80 
81  template <typename T> T get_relativeConductivity(T head) {
82  const T S_e = get_effSaturation(head);
83  const T F_e = get_Fe(S_e);
84  const T F_1 = get_Fe(1.0);
85  return pow(S_e, ll) * pow(((1 - F_e) / (1 - F_1)), 2);
86  }
87 
88  template <typename T> T get_Fe(T eff_satu) {
89  const T m = 1 - 1.0 / nn;
90  const T S_eStar = get_eff_satuStar(eff_satu);
91  return pow(1 - pow(S_eStar, 1.0 / m), m);
92  }
93  template <typename T> T get_eff_satuStar(T eff_satu) {
94  return (theta_s - theta_r) / (theta_m - theta_r) * eff_satu;
95  }
96 
97  template <typename T> T get_effSaturation(T head) {
98  const T theta = get_waterContent(head);
99  return (theta - theta_r) / (theta_s - theta_r);
100  }
101 
102  template <typename T> T get_capacity(T head, double head_norm) {
103  const T m = 1 - 1.0 / nn;
104  T ret_val;
105  if (head_norm < h_s) {
106  ret_val = alpha * m * nn * (theta_m - theta_r) *
107  pow(-alpha * head, nn - 1) /
108  pow(1 + pow(-alpha * head, nn), m + 1);
109  } else {
110  ret_val = 0;
111  }
112  return ret_val;
113  }
114 
115  double get_diffCapacity(double head) {
116  std::complex<double> cx_head = head + eps * 1i;
117  auto capacity = get_capacity(cx_head, head);
118  return capacity.imag() / eps;
119  }
120 
121  double get_diffConductivity(double head) {
122  std::complex<double> cx_head = head + eps * 1i;
123  auto conductivity = get_conductivity(cx_head, head);
124  return conductivity.imag() / eps;
125  }
126 };
127 
129 
130 struct PreviousData {
131  MatrixDouble grads; ///< Gradients of field "u" at integration points
132  VectorDouble values; ///< Values of field "u" at integration points
134  dot_values; ///< Rate of values of field "u" at integration points
136 };
137 
139 
140  OpAssembleStiffRhs(std::string field, boost::shared_ptr<PreviousData> &data)
141  : OpVolEle(field, OpVolEle::OPROW), commonData(data) {}
142  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
144  const int nb_dofs = data.getIndices().size();
145  if (nb_dofs) {
146  vecF.resize(nb_dofs, false);
147  vecF.clear();
148  const int nb_integration_pts = getGaussPts().size2();
149 
150  auto t_dot_val = getFTensor0FromVec(commonData->dot_values);
151  auto t_val = getFTensor0FromVec(commonData->values);
152  auto t_grad = getFTensor1FromMat<2>(commonData->grads);
153 
154  auto t_base = data.getFTensor0N();
155  auto t_diff_base = data.getFTensor1DiffN<2>();
156  auto t_w = getFTensor0IntegrationWeight();
157 
159 
160  const double vol = getMeasure();
161 
162  for (int gg = 0; gg != nb_integration_pts; ++gg) {
163  const double a = vol * t_w;
164 
165  const double K_h = block_map.get_conductivity<double>(t_val, t_val);
166  const double C_h = block_map.get_capacity<double>(t_val, t_val);
167 
168  for (int rr = 0; rr != nb_dofs; ++rr) {
169  vecF[rr] +=
170  a * (t_base * C_h * t_dot_val + K_h * t_diff_base(i) * t_grad(i));
171  ++t_diff_base;
172  ++t_base;
173  }
174 
175  ++t_dot_val;
176  ++t_grad;
177  ++t_val;
178  ++t_w;
179  }
180  CHKERR VecSetOption(getFEMethod()->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
181  PETSC_TRUE);
182  CHKERR VecSetValues(getFEMethod()->ts_F, data, &*vecF.begin(),
183  ADD_VALUES);
184  }
186  }
187 
188 private:
189  boost::shared_ptr<PreviousData> commonData;
191  std::string field;
192 };
193 
195 
196  OpAssembleStiffLhs(std::string fieldu, boost::shared_ptr<PreviousData> &data)
197  : OpVolEle(fieldu, fieldu, OpVolEle::OPROWCOL), commonData(data) {
198  sYmm = false;
199  }
200  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
201  EntityType col_type, EntData &row_data,
202  EntData &col_data) {
204 
205  const int nb_row_dofs = row_data.getIndices().size();
206  const int nb_col_dofs = col_data.getIndices().size();
207  // cerr << "In doWork() : (row, col) = (" << nb_row_dofs << ", " <<
208  // nb_col_dofs << ")" << endl;
209  if (nb_row_dofs && nb_col_dofs) {
210  mat.resize(nb_row_dofs, nb_col_dofs, false);
211  mat.clear();
212  const int nb_integration_pts = getGaussPts().size2();
213  auto t_row_base = row_data.getFTensor0N();
214 
215  auto t_val = getFTensor0FromVec(commonData->values);
216  auto t_dot_val = getFTensor0FromVec(commonData->dot_values);
217  auto t_grad = getFTensor1FromMat<2>(commonData->grads);
218 
219  auto t_row_diff_base = row_data.getFTensor1DiffN<2>();
220  auto t_w = getFTensor0IntegrationWeight();
221  const double ts_a = getFEMethod()->ts_a;
222  const double vol = getMeasure();
223 
225 
226  for (int gg = 0; gg != nb_integration_pts; ++gg) {
227  const double a = vol * t_w;
228  const double K_h = block_map.get_conductivity<double>(t_val, t_val);
229  const double C_h = block_map.get_capacity<double>(t_val, t_val);
230  const double DK_h = block_map.get_diffConductivity(t_val);
231  const double DC_h = block_map.get_diffCapacity(t_val);
232 
233  for (int rr = 0; rr != nb_row_dofs; ++rr) {
234  auto t_col_base = col_data.getFTensor0N(gg, 0);
235  auto t_col_diff_base = col_data.getFTensor1DiffN<2>(gg, 0);
236  for (int cc = 0; cc != nb_col_dofs; ++cc) {
237 
238  mat(rr, cc) +=
239  a * (t_row_base * (DC_h * t_dot_val + C_h * ts_a) * t_col_base +
240  DK_h * t_grad(i) * t_row_diff_base(i) * t_col_base +
241  K_h * t_row_diff_base(i) * t_col_diff_base(i));
242 
243  ++t_col_base;
244  ++t_col_diff_base;
245  }
246  ++t_row_base;
247  ++t_row_diff_base;
248  }
249  ++t_w;
250  ++t_val;
251  ++t_grad;
252  ++t_dot_val;
253  }
254  CHKERR MatSetValues(getFEMethod()->ts_B, row_data, col_data, &mat(0, 0),
255  ADD_VALUES);
256  // if (row_side != col_side || row_type != col_type) {
257  // transMat.resize(nb_col_dofs, nb_row_dofs, false);
258  // noalias(transMat) = trans(mat);
259  // CHKERR MatSetValues(getFEMethod()->ts_B, col_data, row_data,
260  // &transMat(0, 0), ADD_VALUES);
261  // }
262  }
264  }
265 
266 private:
267  boost::shared_ptr<PreviousData> commonData;
269 };
270 
272 {
273  OpAssembleNaturalBCRhs(std::string mass_field, Range &natural_bd_ents)
274  : OpFaceEle(mass_field, OpFaceEle::OPROW),
276 
277  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
279  const int nb_dofs = data.getIndices().size();
280 
281  if (nb_dofs) {
282  EntityHandle row_side_ent = getFEEntityHandle();
283  bool is_natural =
284  (natural_bd_ents.find(row_side_ent) != natural_bd_ents.end());
285  if (is_natural) {
286  // cerr << "In NaturalBCRhsTau..." << endl;
287  vecF.resize(nb_dofs, false);
288  vecF.clear();
289  const int nb_integration_pts = getGaussPts().size2();
290  auto t_row_base = data.getFTensor0N();
291 
292  auto t_w = getFTensor0IntegrationWeight();
293  const double vol = getMeasure();
294 
295  for (int gg = 0; gg != nb_integration_pts; ++gg) {
296  const double a = vol * t_w;
297 
298  double h = natural_bc_values;
299  for (int rr = 0; rr != nb_dofs; ++rr) {
300  vecF[rr] += t_row_base * h * a;
301  ++t_row_base;
302  }
303  // cout << "vecF : " << vecF << endl;
304  ++t_w;
305  }
306  CHKERR VecSetOption(getFEMethod()->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
307  PETSC_TRUE);
308  CHKERR VecSetValues(getFEMethod()->ts_F, data, &*vecF.begin(),
309  ADD_VALUES);
310  }
311  }
313  }
314 
315 private:
318 };
319 
320 struct Monitor : public FEMethod {
321  double &eRror;
322  Monitor(MPI_Comm &comm, const int &rank, SmartPetscObj<DM> &dm,
323  boost::shared_ptr<PostProc> &post_proc, double &err)
324  : cOmm(comm), rAnk(rank), dM(dm), postProc(post_proc), eRror(err){};
325  MoFEMErrorCode preProcess() { return 0; }
326  MoFEMErrorCode operator()() { return 0; }
327 
330  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-save_every_nth_step", &nth_step,
331  PETSC_NULL);
332  if (ts_step % nth_step == 0) {
334  CHKERR postProc->writeFile(
335  "out_level_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
336  }
337 
339  }
340 
341 private:
342  SmartPetscObj<DM> dM;
343  boost::shared_ptr<PostProc> postProc;
344  MPI_Comm cOmm;
345  const int rAnk;
346 };
347 
348 }; // namespace UFOperators2D
349 
350 #endif
constexpr double a
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:496
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:544
EdgeElementForcesAndSourcesCoreSwitch< 0 > EdgeElementForcesAndSourcesCore
Edge finite element default.
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
const double T
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
UBlasVector< double > VectorDouble
Definition: Types.hpp:79
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:149
constexpr double natural_bc_values
MoFEM::FaceElementForcesAndSourcesCore VolEle
DataForcesAndSourcesCore::EntData EntData
FTensor::Index< 'i', 3 > i
constexpr double essential_bc_values
constexpr double eps
double h
convective heat coefficient
FTensor::Index< 'm', 3 > m
bool sYmm
If true assume that matrix is symmetric structure.
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
auto getFTensor0IntegrationWeight()
Get integration weights.
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
PetscReal ts_a
shift for U_t (see PETSc Time Solver)
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
Postprocess on face.
double get_diffCapacity(double head)
T get_conductivity(T head, double head_norm)
T get_capacity(T head, double head_norm)
double get_diffConductivity(double head)
Monitor(MPI_Comm &comm, const int &rank, SmartPetscObj< DM > &dm, boost::shared_ptr< PostProc > &post_proc, double &err)
boost::shared_ptr< PostProc > postProc
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpAssembleNaturalBCRhs(std::string mass_field, Range &natural_bd_ents)
boost::shared_ptr< PreviousData > commonData
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
OpAssembleStiffLhs(std::string fieldu, boost::shared_ptr< PreviousData > &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::shared_ptr< PreviousData > commonData
OpAssembleStiffRhs(std::string field, boost::shared_ptr< PreviousData > &data)
MatrixDouble grads
Gradients of field "u" at integration points.
VectorDouble values
Values of field "u" at integration points.
VectorDouble dot_values
Rate of values of field "u" at integration points.