v0.13.1
Loading...
Searching...
No Matches
Unsaturated2DFlowUDOP.hpp
Go to the documentation of this file.
1#ifndef __UFOPERATORS2D_HPP__
2#define __UFOPERATORS2D_HPP__
3
4#include <stdlib.h>
6
7namespace UFOperators2D {
8
11
14
15using EntData = DataForcesAndSourcesCore::EntData;
16
18
19int nth_step = 4;
20constexpr double natural_bc_values = -1;
21constexpr double essential_bc_values = 0.0;
22constexpr double eps = 1e-16;
23
24// const int 2 = 3;
26struct BlockData {
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
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();
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
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
188private:
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
266private:
267 boost::shared_ptr<PreviousData> commonData;
269};
270
272{
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
315private:
318};
319
320struct 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){};
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
341private:
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:346
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
#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
FTensor::Index< 'm', SPACE_DIM > m
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:533
const double T
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
constexpr double natural_bc_values
DataForcesAndSourcesCore::EntData EntData
FTensor::Index< 'i', 3 > i
constexpr double essential_bc_values
constexpr double eps
double h
bool sYmm
If true assume that matrix is symmetric structure.
structure for User Loop Methods on finite elements
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
auto getFTensor0IntegrationWeight()
Get integration weights.
@ OPROW
operator doWork function is executed on FE rows
@ OPROWCOL
operator doWork is executed on FE rows &columns
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
intrusive_ptr for managing petsc objects
PetscReal ts_a
shift for U_t (see PETSc Time Solver)
PetscInt ts_step
time step number
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)
MoFEMErrorCode postProcess()
function is run at the end of loop
MoFEMErrorCode operator()()
function is run for every finite element
boost::shared_ptr< PostProc > postProc
MoFEMErrorCode preProcess()
function is run at the beginning of loop
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.