Flow123d  master-f44eb46
darcy_flow_mh.cc
Go to the documentation of this file.
1 /*!
2  *
3  * Copyright (C) 2015 Technical University of Liberec. All rights reserved.
4  *
5  * This program is free software; you can redistribute it and/or modify it under
6  * the terms of the GNU General Public License version 3 as published by the
7  * Free Software Foundation. (http://www.gnu.org/licenses/gpl-3.0.en.html)
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
12  *
13  *
14  * @file darcy_flow_mh.cc
15  * @ingroup flow
16  * @brief Setup and solve linear system of mixed-hybrid discretization of the linear
17  * porous media flow with possible preferential flow in fractures and chanels.
18  */
19 
20 //#include <limits>
21 #include <vector>
22 //#include <iostream>
23 //#include <iterator>
24 //#include <algorithm>
25 #include <armadillo>
26 
27 #include "petscmat.h"
28 #include "petscviewer.h"
29 #include "petscerror.h"
30 #include "mpi.h"
31 
32 #include "system/system.hh"
33 #include "system/sys_profiler.hh"
34 #include "system/index_types.hh"
35 #include "input/factory.hh"
36 
37 #include "mesh/mesh.h"
38 #include "mesh/bc_mesh.hh"
39 #include "mesh/partitioning.hh"
40 #include "mesh/accessors.hh"
41 #include "mesh/range_wrapper.hh"
42 #include "la/distribution.hh"
43 #include "la/linsys.hh"
44 #include "la/linsys_PETSC.hh"
45 #include "la/linsys_BDDC.hh"
46 #include "la/schur.hh"
47 //#include "la/sparse_graph.hh"
49 #include "la/vector_mpi.hh"
50 
51 #include "flow/assembly_mh_old.hh"
52 #include "flow/darcy_flow_mh.hh"
54 #include "flow/assembly_models.hh"
55 
56 #include "tools/time_governor.hh"
58 #include "fields/field.hh"
59 #include "fields/field_values.hh"
61 #include "fields/field_fe.hh"
62 #include "fields/field_model.hh"
63 #include "fields/field_constant.hh"
64 
65 #include "coupling/balance.hh"
66 
69 
70 #include "fem/fe_p.hh"
71 
72 
73 FLOW123D_FORCE_LINK_IN_CHILD(darcy_flow_mh)
74 
75 
76 
77 
78 
79 namespace it = Input::Type;
80 
82  return it::Selection("MH_MortarMethod")
83  .add_value(NoMortar, "None", "No Mortar method is applied.")
84  .add_value(MortarP0, "P0", "Mortar space: P0 on elements of lower dimension.")
85  .add_value(MortarP1, "P1", "Mortar space: P1 on intersections, using non-conforming pressures.")
86  .close();
87 }
88 
89 
91  return it::Selection("Flow_Darcy_BC_Type")
92  .add_value(none, "none",
93  "Homogeneous Neumann boundary condition\n(zero normal flux over the boundary).")
94  .add_value(dirichlet, "dirichlet",
95  "Dirichlet boundary condition. "
96  "Specify the pressure head through the ``bc_pressure`` field "
97  "or the piezometric head through the ``bc_piezo_head`` field.")
98  .add_value(total_flux, "total_flux", "Flux boundary condition (combines Neumann and Robin type). "
99  "Water inflow equal to (($ \\delta_d(q_d^N + \\sigma_d (h_d^R - h_d) )$)). "
100  "Specify the water inflow by the ``bc_flux`` field, the transition coefficient by ``bc_robin_sigma`` "
101  "and the reference pressure head or piezometric head through ``bc_pressure`` or ``bc_piezo_head`` respectively.")
102  .add_value(seepage, "seepage",
103  "Seepage face boundary condition. Pressure and inflow bounded from above. Boundary with potential seepage flow "
104  "is described by the pair of inequalities: "
105  "(($h_d \\le h_d^D$)) and (($ -\\boldsymbol q_d\\cdot\\boldsymbol n \\le \\delta q_d^N$)), where the equality holds in at least one of them. "
106  "Caution: setting (($q_d^N$)) strictly negative "
107  "may lead to an ill posed problem since a positive outflow is enforced. "
108  "Parameters (($h_d^D$)) and (($q_d^N$)) are given by the fields ``bc_switch_pressure`` (or ``bc_switch_piezo_head``) and ``bc_flux`` respectively."
109  )
110  .add_value(river, "river",
111  "River boundary condition. For the water level above the bedrock, (($H_d > H_d^S$)), the Robin boundary condition is used with the inflow given by: "
112  "(( $ \\delta_d(q_d^N + \\sigma_d(H_d^D - H_d) )$)). For the water level under the bedrock, constant infiltration is used: "
113  "(( $ \\delta_d(q_d^N + \\sigma_d(H_d^D - H_d^S) )$)). Parameters: ``bc_pressure``, ``bc_switch_pressure``, "
114  " ``bc_sigma``, ``bc_flux``."
115  )
116  .close();
117 }
118 
120 
121  const it::Record &field_descriptor =
123  .copy_keys( DarcyMH::EqFields().make_field_descriptor_type(equation_name() + "_Data_aux") )
124  .declare_key("bc_piezo_head", FieldAlgorithmBase< 3, FieldValue<3>::Scalar >::get_input_type_instance(),
125  "Boundary piezometric head for BC types: dirichlet, robin, and river." )
126  .declare_key("bc_switch_piezo_head", FieldAlgorithmBase< 3, FieldValue<3>::Scalar >::get_input_type_instance(),
127  "Boundary switch piezometric head for BC types: seepage, river." )
128  .declare_key("init_piezo_head", FieldAlgorithmBase< 3, FieldValue<3>::Scalar >::get_input_type_instance(),
129  "Initial condition for the pressure given as the piezometric head." )
130  .close();
131  return field_descriptor;
132 }
133 
135  it::Record ns_rec = Input::Type::Record("NonlinearSolver", "Non-linear solver settings.")
136  .declare_key("linear_solver", LinSys::get_input_type(), it::Default("{}"),
137  "Linear solver for MH problem.")
138  .declare_key("tolerance", it::Double(0.0), it::Default("1E-6"),
139  "Residual tolerance.")
140  .declare_key("min_it", it::Integer(0), it::Default("1"),
141  "Minimum number of iterations (linear solutions) to use.\nThis is usefull if the convergence criteria "
142  "does not characterize your goal well enough so it converges prematurely, possibly even without a single linear solution."
143  "If greater then 'max_it' the value is set to 'max_it'.")
144  .declare_key("max_it", it::Integer(0), it::Default("100"),
145  "Maximum number of iterations (linear solutions) of the non-linear solver.")
146  .declare_key("converge_on_stagnation", it::Bool(), it::Default("false"),
147  "If a stagnation of the nonlinear solver is detected the solver stops. "
148  "A divergence is reported by default, forcing the end of the simulation. By setting this flag to 'true', the solver "
149  "ends with convergence success on stagnation, but it reports warning about it.")
150  .close();
151 
153 
154  return it::Record(equation_name(), "Mixed-Hybrid solver for saturated Darcy flow.")
157  .declare_key("gravity", it::Array(it::Double(), 3,3), it::Default("[ 0, 0, -1]"),
158  "Vector of the gravity force. Dimensionless.")
159 // .declare_key("user_fields", it::Array(DarcyMH::EqFields().make_user_field_type(equation_name())),
160 // IT::Default::optional(),
161 // "Input fields of the equation defined by user.")
163  "Input data for Darcy flow model.")
164  .declare_key("nonlinear_solver", ns_rec, it::Default("{}"),
165  "Non-linear solver for MH problem.")
166  .declare_key("output_stream", OutputTime::get_input_type(), it::Default("{}"),
167  "Output stream settings.\n Specify file format, precision etc.")
168 
170  IT::Default("{ \"fields\": [ \"pressure_p0\", \"velocity_p0\" ] }"),
171  "Specification of output fields and output times.")
173  "Output settings specific to Darcy flow model.\n"
174  "Includes raw output and some experimental functionality.")
175  .declare_key("balance", Balance::get_input_type(), it::Default("{}"),
176  "Settings for computing mass balance.")
177  .declare_key("n_schurs", it::Integer(0,2), it::Default("2"),
178  "Number of Schur complements to perform when solving MH system.")
179  .declare_key("mortar_method", get_mh_mortar_selection(), it::Default("\"None\""),
180  "Method for coupling Darcy flow between dimensions on incompatible meshes. [Experimental]" )
181  .close();
182 }
183 
184 
185 const int DarcyMH::registrar =
186  Input::register_class< DarcyMH, Mesh &, const Input::Record >(equation_name()) +
188 
189 
190 
192 {
193  *this += field_ele_pressure.name("pressure_p0")
194  .units(UnitSI().m())
196  .description("Pressure solution - P0 interpolation.");
197 
198  *this += field_edge_pressure.name("pressure_edge")
199  .units(UnitSI().m())
200  .flags(FieldFlag::input_copy)
201  .description("Pressure solution - Crouzeix-Raviart interpolation.");
202 
203  *this += field_ele_piezo_head.name("piezo_head_p0")
204  .units(UnitSI().m())
206  .description("Piezo head solution - P0 interpolation.");
207 
208  *this += field_ele_velocity.name("velocity_p0")
209  .units(UnitSI().m().s(-1))
211  .description("Velocity solution - P0 interpolation.");
212 
213  *this += flux.name("flux")
214  .units(UnitSI().m().s(-1))
216  .description("Darcy flow flux.");
217 
218  *this += anisotropy.name("anisotropy")
219  .description("Anisotropy of the conductivity tensor.")
220  .input_default("1.0")
221  .units( UnitSI::dimensionless() );
222 
223  *this += cross_section.name("cross_section")
224  .description("Complement dimension parameter (cross section for 1D, thickness for 2D).")
225  .input_default("1.0")
226  .units( UnitSI().m(3).md() );
227 
228  *this += conductivity.name("conductivity")
229  .description("Isotropic conductivity scalar.")
230  .input_default("1.0")
231  .units( UnitSI().m().s(-1) )
232  .set_limits(0.0);
233 
234  *this += sigma.name("sigma")
235  .description("Transition coefficient between dimensions.")
236  .input_default("1.0")
237  .units( UnitSI::dimensionless() );
238 
239  *this += water_source_density.name("water_source_density")
240  .description("Water source density.")
241  .input_default("0.0")
242  .units( UnitSI().s(-1) );
243 
244  *this += bc_type.name("bc_type")
245  .description("Boundary condition type.")
246  .input_selection( get_bc_type_selection() )
247  .input_default("\"none\"")
248  .units( UnitSI::dimensionless() );
249 
250  *this += bc_pressure
251  .disable_where(bc_type, {none, seepage} )
252  .name("bc_pressure")
253  .description("Prescribed pressure value on the boundary. Used for all values of ``bc_type`` except ``none`` and ``seepage``. "
254  "See documentation of ``bc_type`` for exact meaning of ``bc_pressure`` in individual boundary condition types.")
255  .input_default("0.0")
256  .units( UnitSI().m() );
257 
258  *this += bc_flux
259  .disable_where(bc_type, {none, dirichlet} )
260  .name("bc_flux")
261  .description("Incoming water boundary flux. Used for bc_types : ``total_flux``, ``seepage``, ``river``.")
262  .input_default("0.0")
263  .units( UnitSI().m().s(-1) );
264 
265  *this += bc_robin_sigma
266  .disable_where(bc_type, {none, dirichlet, seepage} )
267  .name("bc_robin_sigma")
268  .description("Conductivity coefficient in the ``total_flux`` or the ``river`` boundary condition type.")
269  .input_default("0.0")
270  .units( UnitSI().s(-1) );
271 
272  *this += bc_switch_pressure
273  .disable_where(bc_type, {none, dirichlet, total_flux} )
274  .name("bc_switch_pressure")
275  .description("Critical switch pressure for ``seepage`` and ``river`` boundary conditions.")
276  .input_default("0.0")
277  .units( UnitSI().m() );
278 
279 
280  //these are for unsteady
281  *this += init_pressure.name("init_pressure")
282  .description("Initial condition for pressure in time dependent problems.")
283  .input_default("0.0")
284  .units( UnitSI().m() );
285 
286  *this += storativity.name("storativity")
287  .description("Storativity (in time dependent problems).")
288  .input_default("0.0")
289  .units( UnitSI().m(-1) );
290 
291  *this += extra_storativity.name("extra_storativity")
292  .description("Storativity added from upstream equation.")
293  .units( UnitSI().m(-1) )
294  .input_default("0.0")
295  .flags( input_copy );
296 
297  *this += extra_source.name("extra_water_source_density")
298  .description("Water source density added from upstream equation.")
299  .input_default("0.0")
300  .units( UnitSI().s(-1) )
301  .flags( input_copy );
302 
303  *this += gravity_field.name("gravity")
304  .description("Gravity vector.")
305  .input_default("0.0")
306  .units( UnitSI::dimensionless() );
307 
308  *this += bc_gravity.name("bc_gravity")
309  .description("Boundary gravity vector.")
310  .input_default("0.0")
311  .units( UnitSI::dimensionless() );
312 
313  *this += init_piezo_head.name("init_piezo_head")
314  .units(UnitSI().m())
315  .input_default("0.0")
316  .description("Init piezo head.");
317 
318  *this += bc_piezo_head.name("bc_piezo_head")
319  .units(UnitSI().m())
320  .input_default("0.0")
321  .description("Boundary piezo head.");
322 
323  *this += bc_switch_piezo_head.name("bc_switch_piezo_head")
324  .units(UnitSI().m())
325  .input_default("0.0")
326  .description("Boundary switch piezo head.");
327 
328  //time_term_fields = this->subset({"storativity"});
329  //main_matrix_fields = this->subset({"anisotropy", "conductivity", "cross_section", "sigma", "bc_type", "bc_robin_sigma"});
330  //rhs_fields = this->subset({"water_source_density", "bc_pressure", "bc_flux"});
331 }
332 
334 {
335  mortar_method_=NoMortar;
336 }
337 
338 
339 
340 
341 //=============================================================================
342 // CREATE AND FILL GLOBAL MH MATRIX OF THE WATER MODEL
343 // - do it in parallel:
344 // - initial distribution of elements, edges
345 //
346 /*! @brief CREATE AND FILL GLOBAL MH MATRIX OF THE WATER MODEL
347  *
348  * Parameters {Solver,NSchurs} number of performed Schur
349  * complements (0,1,2) for water flow MH-system
350  *
351  */
352 //=============================================================================
353 DarcyMH::DarcyMH(Mesh &mesh_in, const Input::Record in_rec, TimeGovernor *tm)
354 : DarcyFlowInterface(mesh_in, in_rec),
355  output_object(nullptr),
356  data_changed_(false),
357  schur0(nullptr),
358  steady_diagonal(nullptr),
359  steady_rhs(nullptr),
360  new_diagonal(nullptr),
361  previous_solution(nullptr)
362 {
363 
364  START_TIMER("Darcy constructor");
365  {
366  auto time_rec = in_rec.val<Input::Record>("time");
367  if (tm == nullptr)
368  {
369  time_ = new TimeGovernor(time_rec);
370  }
371  else
372  {
373  TimeGovernor tm_from_rec(time_rec);
374  if (!tm_from_rec.is_default()) // is_default() == false when time record is present in input file
375  {
376  MessageOut() << "Duplicate key 'time', time in flow equation is already initialized from parent class!";
377  ASSERT_PERMANENT(false);
378  }
379  time_ = tm;
380  }
381  }
382 
383  eq_fields_ = make_shared<EqFields>();
384  eq_data_ = make_shared<EqData>();
385  this->eq_fieldset_ = eq_fields_.get();
386 
387  eq_data_->is_linear=true;
388 
389  size = mesh_->n_elements() + mesh_->n_sides() + mesh_->n_edges();
390  n_schur_compls = in_rec.val<int>("n_schurs");
391  eq_data_->mortar_method_= in_rec.val<MortarMethod>("mortar_method");
392  if (eq_data_->mortar_method_ != NoMortar) {
394  }
395 
396 
397 
398  //side_ds->view( std::cout );
399  //el_ds->view( std::cout );
400  //edge_ds->view( std::cout );
401  //rows_ds->view( std::cout );
402 
403 }
404 
406 {
407  // DebugOut() << "t = " << time_->t() << " step_end " << time_->step().end() << "\n";
409  {
410  // In steady case, the solution is computed with the data present at time t,
411  // and the steady state solution is valid until another change in data,
412  // which should correspond to time (t+dt).
413  // "The data change appears immediatly."
414  double next_t = time_->t() + time_->estimate_dt();
415  // DebugOut() << "STEADY next_t = " << next_t << "\n";
416  return next_t * (1 - 2*std::numeric_limits<double>::epsilon());
417  }
418  else
419  {
420  // In unsteady case, the solution is computed with the data present at time t,
421  // and the solution is valid at the time t+dt.
422  // "The data change does not appear immediatly, it is integrated over time interval dt."
423  // DebugOut() << "UNSTEADY\n";
424  return time_->t();
425  }
426 }
427 
429 //connecting data fields with mesh
430 {
431 
432  START_TIMER("data init");
433  eq_data_->mesh = mesh_;
434  eq_fields_->set_mesh(*mesh_);
435 
436  auto gravity_array = input_record_.val<Input::Array>("gravity");
437  std::vector<double> gvec;
438  gravity_array.copy_to(gvec);
439  gvec.push_back(0.0); // zero pressure shift
440  eq_data_->gravity_ = arma::vec(gvec);
441  eq_data_->gravity_vec_ = eq_data_->gravity_.subvec(0,2);
442 
443  FieldValue<3>::VectorFixed gvalue(eq_data_->gravity_vec_);
444  auto field_algo=std::make_shared<FieldConstant<3, FieldValue<3>::VectorFixed>>();
445  field_algo->set_value(gvalue);
446  eq_fields_->gravity_field.set(field_algo, 0.0);
447  eq_fields_->bc_gravity.set(field_algo, 0.0);
448 
449  eq_fields_->bc_pressure.add_factory(
450  std::make_shared<AddPotentialFactory<3, FieldValue<3>::Scalar> >
451  (eq_fields_->bc_gravity, eq_fields_->X(), eq_fields_->bc_piezo_head) );
452  eq_fields_->bc_switch_pressure.add_factory(
453  std::make_shared<AddPotentialFactory<3, FieldValue<3>::Scalar> >
454  (eq_fields_->bc_gravity, eq_fields_->X(), eq_fields_->bc_switch_piezo_head) );
455  eq_fields_->init_pressure.add_factory(
456  std::make_shared<AddPotentialFactory<3, FieldValue<3>::Scalar> >
457  (eq_fields_->gravity_field, eq_fields_->X(), eq_fields_->init_piezo_head) );
458 
459 
460  eq_fields_->set_input_list( this->input_record_.val<Input::Array>("input_fields"), *time_ );
461  // read optional user fields
462 // Input::Array user_fields_arr;
463 // if (input_record_.opt_val("user_fields", user_fields_arr)) {
464 // eq_fields_->init_user_fields(user_fields_arr, time_->step());
465 // }
466 
467  // Check that the time step was set for the transient simulation.
468  if (! zero_time_term(true) && time_->is_default() ) {
469  //THROW(ExcAssertMsg());
470  //THROW(ExcMissingTimeGovernor() << input_record_.ei_address());
471  MessageOut() << "Missing the key 'time', obligatory for the transient problems." << endl;
472  ASSERT_PERMANENT(false);
473  }
474 
475  eq_fields_->mark_input_times(*time_);
476 }
477 
478 
479 
481 
482  { // init DOF handler for pressure fields
483 // std::shared_ptr< FiniteElement<0> > fe0_rt = std::make_shared<FE_RT0_disc<0>>();
484  std::shared_ptr< FiniteElement<1> > fe1_rt = std::make_shared<FE_RT0_disc<1>>();
485  std::shared_ptr< FiniteElement<2> > fe2_rt = std::make_shared<FE_RT0_disc<2>>();
486  std::shared_ptr< FiniteElement<3> > fe3_rt = std::make_shared<FE_RT0_disc<3>>();
487  std::shared_ptr< FiniteElement<0> > fe0_disc = std::make_shared<FE_P_disc<0>>(0);
488  std::shared_ptr< FiniteElement<1> > fe1_disc = std::make_shared<FE_P_disc<1>>(0);
489  std::shared_ptr< FiniteElement<2> > fe2_disc = std::make_shared<FE_P_disc<2>>(0);
490  std::shared_ptr< FiniteElement<3> > fe3_disc = std::make_shared<FE_P_disc<3>>(0);
491  std::shared_ptr< FiniteElement<0> > fe0_cr = std::make_shared<FE_CR<0>>();
492  std::shared_ptr< FiniteElement<1> > fe1_cr = std::make_shared<FE_CR<1>>();
493  std::shared_ptr< FiniteElement<2> > fe2_cr = std::make_shared<FE_CR<2>>();
494  std::shared_ptr< FiniteElement<3> > fe3_cr = std::make_shared<FE_CR<3>>();
495 // static FiniteElement<0> fe0_sys = FE_P_disc<0>(0); //TODO fix and use solution with FESystem<0>( {fe0_rt, fe0_disc, fe0_cr} )
496  FESystem<0> fe0_sys( {fe0_disc, fe0_disc, fe0_cr} );
497  FESystem<1> fe1_sys( {fe1_rt, fe1_disc, fe1_cr} );
498  FESystem<2> fe2_sys( {fe2_rt, fe2_disc, fe2_cr} );
499  FESystem<3> fe3_sys( {fe3_rt, fe3_disc, fe3_cr} );
500  MixedPtr<FESystem> fe_sys( std::make_shared<FESystem<0>>(fe0_sys), std::make_shared<FESystem<1>>(fe1_sys),
501  std::make_shared<FESystem<2>>(fe2_sys), std::make_shared<FESystem<3>>(fe3_sys) );
502  std::shared_ptr<DiscreteSpace> ds = std::make_shared<EqualOrderDiscreteSpace>( mesh_, fe_sys);
503  eq_data_->dh_ = std::make_shared<DOFHandlerMultiDim>(*mesh_);
504  eq_data_->dh_->distribute_dofs(ds);
505  }
506 
507  init_eq_data();
508  this->eq_data_->multidim_assembler = AssemblyFlowBase::create< AssemblyMH >(eq_fields_, eq_data_);
510 
511  eq_fields_->add_coords_field();
512 
513  { // construct pressure, velocity and piezo head fields
514  uint rt_component = 0;
515  eq_data_->full_solution = eq_data_->dh_->create_vector();
516  auto ele_flux_ptr = create_field_fe<3, FieldValue<3>::VectorFixed>(eq_data_->dh_, &eq_data_->full_solution, rt_component);
517  eq_fields_->flux.set(ele_flux_ptr, 0.0);
518 
519  eq_fields_->field_ele_velocity.set(Model<3, FieldValue<3>::VectorFixed>::create(fn_mh_velocity(), eq_fields_->flux, eq_fields_->cross_section), 0.0);
520 
521  uint p_ele_component = 1;
522  auto ele_pressure_ptr = create_field_fe<3, FieldValue<3>::Scalar>(eq_data_->dh_, &eq_data_->full_solution, p_ele_component);
523  eq_fields_->field_ele_pressure.set(ele_pressure_ptr, 0.0);
524 
525  uint p_edge_component = 2;
526  auto edge_pressure_ptr = create_field_fe<3, FieldValue<3>::Scalar>(eq_data_->dh_, &eq_data_->full_solution, p_edge_component);
527  eq_fields_->field_edge_pressure.set(edge_pressure_ptr, 0.0);
528 
529  eq_fields_->field_ele_piezo_head.set(
530  Model<3, FieldValue<3>::Scalar>::create(fn_mh_piezohead(), eq_fields_->gravity_field, eq_fields_->X(), eq_fields_->field_ele_pressure),
531  0.0
532  );
533  }
534 
535  { // init DOF handlers represents edge DOFs
536  uint p_edge_component = 2;
537  eq_data_->dh_cr_ = std::make_shared<SubDOFHandlerMultiDim>(eq_data_->dh_,p_edge_component);
538  }
539 
540  { // init DOF handlers represents side DOFs
541  MixedPtr<FE_CR_disc> fe_cr_disc;
542  std::shared_ptr<DiscreteSpace> ds_cr_disc = std::make_shared<EqualOrderDiscreteSpace>( mesh_, fe_cr_disc);
543  eq_data_->dh_cr_disc_ = std::make_shared<DOFHandlerMultiDim>(*mesh_);
544  eq_data_->dh_cr_disc_->distribute_dofs(ds_cr_disc);
545  }
546 
547  // Initialize bc_switch_dirichlet to size of global boundary.
548  eq_data_->bc_switch_dirichlet.resize(mesh_->n_elements()+mesh_->bc_mesh()->n_elements(), 1);
549 
550 
553  .val<Input::Record>("nonlinear_solver")
554  .val<Input::AbstractRecord>("linear_solver");
555 
556  // auxiliary set_time call since allocation assembly evaluates fields as well
557  time_->step().use_fparser_ = true;
560 
561 
562 
563  // allocate time term vectors
564  VecDuplicate(schur0->get_solution(), &previous_solution);
565  VecCreateMPI(PETSC_COMM_WORLD, eq_data_->dh_->distr()->lsize(),PETSC_DETERMINE,&(steady_diagonal));
566  VecDuplicate(steady_diagonal,& new_diagonal);
567  VecZeroEntries(new_diagonal);
568  VecDuplicate(steady_diagonal, &steady_rhs);
569 
570 
571  // initialization of balance object
572  balance_ = std::make_shared<Balance>("water", mesh_);
573  balance_->init_from_input(input_record_.val<Input::Record>("balance"), time());
574  eq_data_->water_balance_idx = balance_->add_quantity("water_volume");
575  balance_->allocate(eq_data_->dh_->distr()->lsize(), 1);
576  balance_->units(UnitSI().m(3));
577 
578 
579  eq_data_->balance = balance_;
580  eq_data_->lin_sys = schur0;
581 
582 
584 }
585 
587 {}
588 
590 {
591 
592  /* TODO:
593  * - Allow solution reconstruction (pressure and velocity) from initial condition on user request.
594  * - Steady solution as an intitial condition may be forced by setting inti_time =-1, and set data for the steady solver in that time.
595  * Solver should be able to switch from and to steady case depending on the zero time term.
596  */
597 
598  time_->step().use_fparser_ = true;
600 
601  // zero_time_term means steady case
602  bool zero_time_term_from_right = zero_time_term();
603 
604 
605  if (zero_time_term_from_right) {
606  MessageOut() << "Flow zero time step - steady case\n";
607  // steady case
608  VecZeroEntries(schur0->get_solution());
609  //read_initial_condition(); // Possible solution guess for steady case.
610  use_steady_assembly_ = true;
611  solve_nonlinear(); // with right limit data
612  // eq_data_->full_solution.local_to_ghost_begin();
613  // eq_data_->full_solution.local_to_ghost_end();
614  } else {
615  MessageOut() << "Flow zero time step - unsteady case\n";
616  VecZeroEntries(schur0->get_solution());
617  VecZeroEntries(previous_solution);
618 
620  reconstruct_solution_from_schur(eq_data_->multidim_assembler);
621 
622  assembly_linear_system(); // in particular due to balance
623 // print_matlab_matrix("matrix_zero");
624  }
625  //solution_output(T,right_limit); // data for time T in any case
626 
627  eq_data_->full_solution.local_to_ghost_begin();
628  eq_data_->full_solution.local_to_ghost_end();
629  output_data();
630 }
631 
632 //=============================================================================
633 // COMPOSE and SOLVE WATER MH System possibly through Schur complements
634 //=============================================================================
636 {
637  START_TIMER("Solving MH system");
638 
639  time_->next_time();
640 
641  time_->view("DARCY"); //time governor information output
642 
643  solve_time_step();
644 
645  eq_data_->full_solution.local_to_ghost_begin();
646  eq_data_->full_solution.local_to_ghost_end();
647 }
648 
649 
650 void DarcyMH::solve_time_step(bool output)
651 {
652  time_->step().use_fparser_ = true;
654 
655  bool zero_time_term_from_left=zero_time_term();
656 
657  bool jump_time = eq_fields_->storativity.is_jump_time();
658  if (! zero_time_term_from_left) {
659  MessageOut() << "Flow time step - unsteady case\n";
660  // time term not treated as zero
661  // Unsteady solution up to the T.
662 
663  // this flag is necesssary for switching BC to avoid setting zero neumann on the whole boundary in the steady case
664  use_steady_assembly_ = false;
665 
667  solve_nonlinear(); // with left limit data
668  if (jump_time) {
669  WarningOut() << "Output of solution discontinuous in time not supported yet.\n";
670  //solution_output(T, left_limit); // output use time T- delta*dt
671  //output_data();
672  }
673  }
674 
675  if (time_->is_end()) {
676  // output for unsteady case, end_time should not be the jump time
677  // but rether check that
678  if (! zero_time_term_from_left && ! jump_time && output) output_data();
679  return;
680  }
681 
682  time_->step().use_fparser_ = true;
684  bool zero_time_term_from_right=zero_time_term();
685  if (zero_time_term_from_right) {
686  MessageOut() << "Flow time step - steady case\n";
687  // this flag is necesssary for switching BC to avoid setting zero neumann on the whole boundary in the steady case
688  use_steady_assembly_ = true;
689  solve_nonlinear(); // with right limit data
690 
691  } else if (! zero_time_term_from_left && jump_time) {
692  WarningOut() << "Discontinuous time term not supported yet.\n";
693  //solution_transfer(); // internally call set_time(T, left) and set_time(T,right) again
694  //solve_nonlinear(); // with right limit data
695  }
696 
697  //solution_output(T,right_limit); // data for time T in any case
698  if (output) output_data();
699 }
700 
701 
702 bool DarcyMH::zero_time_term(bool time_global) {
703  if (time_global) {
704  return (eq_fields_->storativity.input_list_size() == 0);
705  } else {
706  return eq_fields_->storativity.field_result(mesh_->region_db().get_region_set("BULK")) == result_zeros;
707  }
708 }
709 
710 
712 {
713 
715  double residual_norm = schur0->compute_residual();
717  MessageOut().fmt("[nonlinear solver] norm of initial residual: {}\n", residual_norm);
718 
719  // Reduce is_linear flag.
720  int is_linear_common;
721  MPI_Allreduce(&(eq_data_->is_linear), &is_linear_common,1, MPI_INT ,MPI_MIN,PETSC_COMM_WORLD);
722 
723  Input::Record nl_solver_rec = input_record_.val<Input::Record>("nonlinear_solver");
724  this->tolerance_ = nl_solver_rec.val<double>("tolerance");
725  this->max_n_it_ = nl_solver_rec.val<unsigned int>("max_it");
726  this->min_n_it_ = nl_solver_rec.val<unsigned int>("min_it");
727  if (this->min_n_it_ > this->max_n_it_) this->min_n_it_ = this->max_n_it_;
728 
729  if (! is_linear_common) {
730  // set tolerances of the linear solver unless they are set by user.
731  schur0->set_tolerances(0.1*this->tolerance_, 0.01*this->tolerance_, 100);
732  }
733  vector<double> convergence_history;
734 
735  Vec save_solution;
736  VecDuplicate(schur0->get_solution(), &save_solution);
737  while (nonlinear_iteration_ < this->min_n_it_ ||
738  (residual_norm > this->tolerance_ && nonlinear_iteration_ < this->max_n_it_ )) {
739  ASSERT_EQ( convergence_history.size(), nonlinear_iteration_ );
740  convergence_history.push_back(residual_norm);
741 
742  // stagnation test
743  if (convergence_history.size() >= 5 &&
744  convergence_history[ convergence_history.size() - 1]/convergence_history[ convergence_history.size() - 2] > 0.9 &&
745  convergence_history[ convergence_history.size() - 1]/convergence_history[ convergence_history.size() - 5] > 0.8) {
746  // stagnation
747  if (input_record_.val<Input::Record>("nonlinear_solver").val<bool>("converge_on_stagnation")) {
748  WarningOut().fmt("Accept solution on stagnation. Its: {} Residual: {}\n", nonlinear_iteration_, residual_norm);
749  break;
750  } else {
751  THROW(ExcSolverDiverge() << EI_Reason("Stagnation."));
752  }
753  }
754 
755  if (! is_linear_common)
756  VecCopy( schur0->get_solution(), save_solution);
759 
760  // hack to make BDDC work with empty compute_residual
761  if (is_linear_common){
762  // we want to print this info in linear (and steady) case
763  residual_norm = schur0->compute_residual();
764  MessageOut().fmt("[nonlinear solver] lin. it: {}, reason: {}, residual: {}\n",
765  si.n_iterations, si.converged_reason, residual_norm);
766  break;
767  }
768  data_changed_=true; // force reassembly for non-linear case
769 
770  double alpha = 1; // how much of new solution
771  VecAXPBY(schur0->get_solution(), (1-alpha), alpha, save_solution);
772 
773  //LogOut().fmt("Linear solver ended with reason: {} \n", si.converged_reason );
774  //ASSERT_PERMANENT_GE( si.converged_reason, 0).error("Linear solver failed to converge. \n");
776 
777  residual_norm = schur0->compute_residual();
778  MessageOut().fmt("[nonlinear solver] it: {} lin. it: {}, reason: {}, residual: {}\n",
779  nonlinear_iteration_, si.n_iterations, si.converged_reason, residual_norm);
780  }
781  chkerr(VecDestroy(&save_solution));
782  this -> postprocess();
783 
784  // adapt timestep
785  if (! this->zero_time_term()) {
786  double mult = 1.0;
787  if (nonlinear_iteration_ < 3) mult = 1.6;
788  if (nonlinear_iteration_ > 7) mult = 0.7;
789  time_->set_upper_constraint(time_->dt() * mult, "Darcy adaptivity.");
790  // int result = time_->set_upper_constraint(time_->dt() * mult, "Darcy adaptivity.");
791  //DebugOut().fmt("time adaptivity, res: {} it: {} m: {} dt: {} edt: {}\n", result, nonlinear_iteration_, mult, time_->dt(), time_->estimate_dt());
792  }
793 }
794 
796 {
798 }
799 
801 {
802  START_TIMER("postprocess");
803 
804  //fix velocity when mortar method is used
805  if(eq_data_->mortar_method_ != MortarMethod::NoMortar){
806  auto multidim_assembler = AssemblyFlowBase::create< AssemblyMH >(eq_fields_, eq_data_);
807  for ( DHCellAccessor dh_cell : eq_data_->dh_->own_range() ) {
808  unsigned int dim = dh_cell.dim();
809  multidim_assembler[dim-1]->fix_velocity(dh_cell);
810  }
811  }
812  //ElementAccessor<3> ele;
813 
814  // modify side fluxes in parallel
815  // for every local edge take time term on digonal and add it to the corresponding flux
816  /*
817  for (unsigned int i_loc = 0; i_loc < el_ds->lsize(); i_loc++) {
818  ele = mesh_->element_accessor(el_4_loc[i_loc]);
819  for (unsigned int i=0; i<ele->n_sides(); i++) {
820  side_rows[i] = side_row_4_id[ mh_dh.side_dof( ele_ac.side(i) ) ];
821  values[i] = -1.0 * ele_ac.measure() *
822  data.cross_section.value(ele_ac.centre(), ele_ac.element_accessor()) *
823  data.water_source_density.value(ele_ac.centre(), ele_ac.element_accessor()) /
824  ele_ac.n_sides();
825  }
826  VecSetValues(schur0->get_solution(), ele_ac.n_sides(), side_rows, values, ADD_VALUES);
827  }
828  VecAssemblyBegin(schur0->get_solution());
829  VecAssemblyEnd(schur0->get_solution());
830  */
831 }
832 
833 
835  START_TIMER("Darcy output data");
836 
837  //print_matlab_matrix("matrix_" + std::to_string(time_->step().index()));
838 
839  //time_->view("DARCY"); //time governor information output
840  this->output_object->output();
841 
842 
843  START_TIMER("Darcy balance output");
844  balance_->calculate_cumulative(eq_data_->water_balance_idx, schur0->get_solution());
845  balance_->calculate_instant(eq_data_->water_balance_idx, schur0->get_solution());
846  balance_->output();
847 }
848 
849 
851 {
852  return schur0->get_solution_precision();
853 }
854 
855 
856 // ===========================================================================================
857 //
858 // MATRIX ASSEMBLY - we use abstract assembly routine, where LS Mat/Vec SetValues
859 // are in fact pointers to allocating or filling functions - this is governed by Linsystem roitunes
860 //
861 // =======================================================================================
863 {
864  START_TIMER("DarcyFlowMHy::assembly_mh_matrix");
865 
866  // set auxiliary flag for switchting Dirichlet like BC
867  eq_data_->force_no_neumann_bc = use_steady_assembly_ && (nonlinear_iteration_ == 0);
868  eq_data_->n_schur_compls = n_schur_compls;
869 
870 
871  balance_->start_flux_assembly(eq_data_->water_balance_idx);
872 
873  // TODO: try to move this into balance, or have it in the generic assembler class, that should perform the cell loop
874  // including various pre- and post-actions
875  for ( DHCellAccessor dh_cell : eq_data_->dh_->own_range() ) {
876  unsigned int dim = dh_cell.dim();
877  assembler[dim-1]->assemble(dh_cell);
878  }
879 
880 
881  balance_->finish_flux_assembly(eq_data_->water_balance_idx);
882 
883 }
884 
885 
887 {
888  START_TIMER("DarcyFlowMH::allocate_mh_matrix");
889 
890  // set auxiliary flag for switchting Dirichlet like BC
891  eq_data_->n_schur_compls = n_schur_compls;
892  LinSys *ls = schur0;
893 
894  // to make space for second schur complement, max. 10 neighbour edges of one el.
895  double zeros[100000];
896  for(int i=0; i<100000; i++) zeros[i] = 0.0;
897 
898  std::vector<LongIdx> tmp_rows;
899  tmp_rows.reserve(200);
900 
901  std::vector<LongIdx> dofs, dofs_ngh;
902  dofs.reserve(eq_data_->dh_->max_elem_dofs());
903  dofs_ngh.reserve(eq_data_->dh_->max_elem_dofs());
904 
905  for ( DHCellAccessor dh_cell : eq_data_->dh_->own_range() ) {
906  ElementAccessor<3> ele = dh_cell.elm();
907 
908  const uint ndofs = dh_cell.n_dofs();
909  dofs.resize(ndofs);
910  dh_cell.get_dof_indices(dofs);
911 
912  // whole local MH matrix
913  ls->mat_set_values(ndofs, dofs.data(), ndofs, dofs.data(), zeros);
914 
915  tmp_rows.clear();
916 
917  // compatible neighborings rows
918  unsigned int n_neighs = ele->n_neighs_vb();
919  for ( DHCellSide neighb_side : dh_cell.neighb_sides() ) {
920  // every compatible connection adds a 2x2 matrix involving
921  // current element pressure and a connected edge pressure
922 
923  // read neighbor dofs (dh dofhandler)
924  // neighbor cell owning neighb_side
925  DHCellAccessor dh_neighb_cell = neighb_side.cell();
926 
927  const uint ndofs_ngh = dh_neighb_cell.n_dofs();
928  dofs_ngh.resize(ndofs_ngh);
929  dh_neighb_cell.get_dof_indices(dofs_ngh);
930 
931  // local index of pedge dof on neighboring cell
932  // (dim+1) is number of edges of higher dim element
933  // TODO: replace with DHCell getter when available for FESystem component
934  const unsigned int t = dh_neighb_cell.n_dofs() - (dh_neighb_cell.dim()+1) + neighb_side.side().side_idx();
935  tmp_rows.push_back(dofs_ngh[t]);
936  }
937 
938  const uint nsides = ele->n_sides();
939  LongIdx * edge_rows = dofs.data() + nsides; // pointer to start of ele
940  // allocate always also for schur 2
941  ls->mat_set_values(nsides+1, edge_rows, n_neighs, tmp_rows.data(), zeros); // (edges, ele) x (neigh edges)
942  ls->mat_set_values(n_neighs, tmp_rows.data(), nsides+1, edge_rows, zeros); // (neigh edges) x (edges, ele)
943  ls->mat_set_values(n_neighs, tmp_rows.data(), n_neighs, tmp_rows.data(), zeros); // (neigh edges) x (neigh edges)
944 
945  tmp_rows.clear();
946 
947  if (eq_data_->mortar_method_ != NoMortar) {
948  auto &isec_list = mesh_->mixed_intersections().element_intersections_[ele.idx()];
949  for(auto &isec : isec_list ) {
950  IntersectionLocalBase *local = isec.second;
951  DHCellAccessor dh_cell_slave = eq_data_->dh_->cell_accessor_from_element(local->bulk_ele_idx());
952 
953  const uint ndofs_slave = dh_cell_slave.n_dofs();
954  dofs_ngh.resize(ndofs_slave);
955  dh_cell_slave.get_dof_indices(dofs_ngh);
956 
957  //DebugOut().fmt("Alloc: {} {}", ele.idx(), local->bulk_ele_idx());
958  for(unsigned int i_side=0; i_side < dh_cell_slave.elm()->n_sides(); i_side++) {
959  // TODO: replace with DHCell getter when available for FESystem component
960  tmp_rows.push_back( dofs_ngh[(ndofs_slave+1)/2+i_side] );
961  //DebugOut() << "aedge" << print_var(tmp_rows[tmp_rows.size()-1]);
962  }
963  }
964  }
965  /*
966  for(unsigned int i_side=0; i_side < ele->n_sides(); i_side++) {
967  DebugOut() << "aedge:" << print_var(edge_rows[i_side]);
968  }*/
969 
970  edge_rows = dofs.data() + nsides +1; // pointer to start of edges
971  ls->mat_set_values(nsides, edge_rows, tmp_rows.size(), tmp_rows.data(), zeros); // master edges x neigh edges
972  ls->mat_set_values(tmp_rows.size(), tmp_rows.data(), nsides, edge_rows, zeros); // neigh edges x master edges
973  ls->mat_set_values(tmp_rows.size(), tmp_rows.data(), tmp_rows.size(), tmp_rows.data(), zeros); // neigh edges x neigh edges
974 
975  }
976 /*
977  // alloc edge diagonal entries
978  if(rank == 0)
979  for( vector<Edge>::iterator edg = mesh_->edges.begin(); edg != mesh_->edges.end(); ++edg) {
980  int edg_idx = mh_dh.row_4_edge[edg->side(0)->edge_idx()];
981 
982 // for( vector<Edge>::iterator edg2 = mesh_->edges.begin(); edg2 != mesh_->edges.end(); ++edg2){
983 // int edg_idx2 = mh_dh.row_4_edge[edg2->side(0)->edge_idx()];
984 // if(edg_idx == edg_idx2){
985 // DBGCOUT(<< "P[ " << rank << " ] " << "edg alloc: " << edg_idx << " " << edg_idx2 << "\n");
986  ls->mat_set_value(edg_idx, edg_idx, 0.0);
987 // }
988 // }
989  }
990  */
991  /*
992  if (mortar_method_ == MortarP0) {
993  P0_CouplingAssembler(*this).assembly(*ls);
994  } else if (mortar_method_ == MortarP1) {
995  P1_CouplingAssembler(*this).assembly(*ls);
996  }*/
997 }
998 
1000 {
1001  START_TIMER("assembly source term");
1002  balance_->start_source_assembly(eq_data_->water_balance_idx);
1003 
1004  std::vector<LongIdx> global_dofs(eq_data_->dh_->max_elem_dofs());
1005 
1006  for ( DHCellAccessor dh_cell : eq_data_->dh_->own_range() ) {
1007  ElementAccessor<3> ele = dh_cell.elm();
1008 
1009  const uint ndofs = dh_cell.n_dofs();
1010  global_dofs.resize(ndofs);
1011  dh_cell.get_dof_indices(global_dofs);
1012 
1013  double cs = eq_fields_->cross_section.value(ele.centre(), ele);
1014 
1015  // set sources
1016  double source = ele.measure() * cs *
1017  (eq_fields_->water_source_density.value(ele.centre(), ele)
1018  +eq_fields_->extra_source.value(ele.centre(), ele));
1019  // TODO: replace with DHCell getter when available for FESystem component
1020  schur0->rhs_set_value(global_dofs[ndofs/2], -1.0 * source );
1021 
1022  balance_->add_source_values(eq_data_->water_balance_idx, ele.region().bulk_idx(),
1023  {dh_cell.get_loc_dof_indices()[ndofs/2]}, {0}, {source});
1024  }
1025 
1026  balance_->finish_source_assembly(eq_data_->water_balance_idx);
1027 }
1028 
1029 
1030 
1031 
1032 /*******************************************************************************
1033  * COMPOSE WATER MH MATRIX WITHOUT SCHUR COMPLEMENT
1034  ******************************************************************************/
1035 
1037 
1038  START_TIMER("preallocation");
1039 
1040  if (schur0 == NULL) { // create Linear System for MH matrix
1041 
1042  if (in_rec.type() == LinSys_BDDC::get_input_type()) {
1043 #ifdef FLOW123D_HAVE_BDDCML
1044  WarningOut() << "For BDDC no Schur complements are used.";
1045  n_schur_compls = 0;
1046  LinSys_BDDC *ls = new LinSys_BDDC(&(*eq_data_->dh_->distr()),
1047  true); // swap signs of matrix and rhs to make the matrix SPD
1048  ls->set_from_input(in_rec);
1049  ls->set_solution( eq_data_->full_solution.petsc_vec() );
1050  // possible initialization particular to BDDC
1051  START_TIMER("BDDC set mesh data");
1053  schur0=ls;
1054  END_TIMER("BDDC set mesh data");
1055 #else
1056  //Exception
1057  THROW( ExcBddcmlNotSupported() );
1058 #endif // FLOW123D_HAVE_BDDCML
1059  }
1060  else if (in_rec.type() == LinSys_PETSC::get_input_type()) {
1061  // use PETSC for serial case even when user wants BDDC
1062  if (n_schur_compls > 2) {
1063  WarningOut() << "Invalid number of Schur Complements. Using 2.";
1064  n_schur_compls = 2;
1065  }
1066 
1067  LinSys_PETSC *schur1, *schur2;
1068 
1069  if (n_schur_compls == 0) {
1070  LinSys_PETSC *ls = new LinSys_PETSC( &(*eq_data_->dh_->distr()) );
1071 
1072  // temporary solution; we have to set precision also for sequantial case of BDDC
1073  // final solution should be probably call of direct solver for oneproc case
1074  if (in_rec.type() != LinSys_BDDC::get_input_type()) ls->set_from_input(in_rec);
1075  else {
1076  ls->LinSys::set_from_input(in_rec); // get only common options
1077  }
1078 
1079  ls->set_solution( eq_data_->full_solution.petsc_vec() );
1080  schur0=ls;
1081  } else {
1082  IS is;
1083  auto side_dofs_vec = get_component_indices_vec(0);
1084 
1085  ISCreateGeneral(PETSC_COMM_SELF, side_dofs_vec.size(), &(side_dofs_vec[0]), PETSC_COPY_VALUES, &is);
1086  //ISView(is, PETSC_VIEWER_STDOUT_SELF);
1087  //ASSERT_PERMANENT(err == 0).error("Error in ISCreateStride.");
1088 
1089  SchurComplement *ls = new SchurComplement(&(*eq_data_->dh_->distr()), is);
1090 
1091  // make schur1
1093  if (n_schur_compls==1) {
1094  schur1 = new LinSys_PETSC(ds);
1095  schur1->set_positive_definite();
1096  } else {
1097  IS is;
1098  auto elem_dofs_vec = get_component_indices_vec(1);
1099 
1100  const PetscInt *b_indices;
1101  ISGetIndices(ls->IsB, &b_indices);
1102  uint b_size = ls->loc_size_B;
1103  for(uint i_b=0, i_bb=0; i_b < b_size && i_bb < elem_dofs_vec.size(); i_b++) {
1104  if (b_indices[i_b] == elem_dofs_vec[i_bb])
1105  elem_dofs_vec[i_bb++] = i_b + ds->begin();
1106  }
1107  ISRestoreIndices(ls->IsB, &b_indices);
1108 
1109 
1110  ISCreateGeneral(PETSC_COMM_SELF, elem_dofs_vec.size(), &(elem_dofs_vec[0]), PETSC_COPY_VALUES, &is);
1111  //ISView(is, PETSC_VIEWER_STDOUT_SELF);
1112  //ASSERT_PERMANENT(err == 0).error("Error in ISCreateStride.");
1113  SchurComplement *ls1 = new SchurComplement(ds, is); // is is deallocated by SchurComplement
1114  ls1->set_negative_definite();
1115 
1116  // make schur2
1117  schur2 = new LinSys_PETSC( ls1->make_complement_distribution() );
1118  schur2->set_positive_definite();
1119  ls1->set_complement( schur2 );
1120  schur1 = ls1;
1121  }
1122  ls->set_complement( schur1 );
1123  ls->set_from_input(in_rec);
1124  ls->set_solution( eq_data_->full_solution.petsc_vec() );
1125  schur0=ls;
1126  }
1127 
1128  START_TIMER("PETSC PREALLOCATION");
1129  schur0->set_symmetric();
1131 
1133 
1134  VecZeroEntries(schur0->get_solution());
1135  END_TIMER("PETSC PREALLOCATION");
1136  }
1137  else {
1138  THROW( ExcUnknownSolver() );
1139  }
1140 
1141  END_TIMER("preallocation");
1142  }
1143 
1144 }
1145 
1146 
1147 
1148 
1150  START_TIMER("DarcyFlowMH_Steady::assembly_linear_system");
1151 
1152  eq_data_->is_linear=true;
1153  bool is_steady = zero_time_term();
1154  //DebugOut() << "Assembly linear system\n";
1155  if (data_changed_) {
1156  data_changed_ = false;
1157  //DebugOut() << "Data changed\n";
1158  // currently we have no optimization for cases when just time term data or RHS data are changed
1159  START_TIMER("full assembly");
1160  if (typeid(*schur0) != typeid(LinSys_BDDC)) {
1161  schur0->start_add_assembly(); // finish allocation and create matrix
1162  }
1163 
1166 
1168 
1169 
1170  assembly_mh_matrix( eq_data_->multidim_assembler ); // fill matrix
1171 
1173 // print_matlab_matrix("matrix");
1175  //MatView( *const_cast<Mat*>(schur0->get_matrix()), PETSC_VIEWER_STDOUT_WORLD );
1176  //VecView( *const_cast<Vec*>(schur0->get_rhs()), PETSC_VIEWER_STDOUT_WORLD);
1177 
1178  if (! is_steady) {
1179  START_TIMER("fix time term");
1180  //DebugOut() << "setup time term\n";
1181  // assembly time term and rhs
1182  setup_time_term();
1183  modify_system();
1184  }
1185  else
1186  {
1187  balance_->start_mass_assembly(eq_data_->water_balance_idx);
1188  balance_->finish_mass_assembly(eq_data_->water_balance_idx);
1189  }
1190  END_TIMER("full assembly");
1191  } else {
1192  START_TIMER("modify system");
1193  if (! is_steady) {
1194  modify_system();
1195  } else {
1196  //Should be replaced with exception if error will be switched on.
1197  //ASSERT_PERMANENT(false).error("Planned computation time for steady solver, but data are not changed.\n");
1198  }
1199  END_TIMER("modiffy system");
1200  }
1201 
1202 }
1203 
1204 
1205 void DarcyMH::print_matlab_matrix(std::string matlab_file)
1206 {
1207  std::string output_file;
1208 
1209  if ( typeid(*schur0) == typeid(LinSys_BDDC) ){
1210 // WarningOut() << "Can output matrix only on a single processor.";
1211 // output_file = FilePath(matlab_file + "_bddc.m", FilePath::output_file);
1212 // ofstream os( output_file );
1213 // auto bddc = static_cast<LinSys_BDDC*>(schur0);
1214 // bddc->print_matrix(os);
1215  }
1216  else {//if ( typeid(*schur0) == typeid(LinSys_PETSC) ){
1217  output_file = FilePath(matlab_file + ".m", FilePath::output_file);
1218  PetscViewer viewer;
1219  PetscViewerASCIIOpen(PETSC_COMM_WORLD, output_file.c_str(), &viewer);
1220  PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
1221  MatView( *const_cast<Mat*>(schur0->get_matrix()), viewer);
1222  VecView( *const_cast<Vec*>(schur0->get_rhs()), viewer);
1223  VecView( *const_cast<Vec*>(schur0->get_rhs()), viewer);
1224  VecView( *const_cast<Vec*>(&(schur0->get_solution())), viewer);
1225  }
1226 // else{
1227 // WarningOut() << "No matrix output available for the current solver.";
1228 // return;
1229 // }
1230 
1231  // compute h_min for different dimensions
1232  double d_max = std::numeric_limits<double>::max();
1233  double h1 = d_max, h2 = d_max, h3 = d_max;
1234  double he2 = d_max, he3 = d_max;
1235  for (auto ele : mesh_->elements_range()) {
1236  switch(ele->dim()){
1237  case 1: h1 = std::min(h1,ele.measure()); break;
1238  case 2: h2 = std::min(h2,ele.measure()); break;
1239  case 3: h3 = std::min(h3,ele.measure()); break;
1240  }
1241 
1242  for (unsigned int j=0; j<ele->n_sides(); j++) {
1243  switch(ele->dim()){
1244  case 2: he2 = std::min(he2, ele.side(j)->measure()); break;
1245  case 3: he3 = std::min(he3, ele.side(j)->measure()); break;
1246  }
1247  }
1248  }
1249  if(h1 == d_max) h1 = 0;
1250  if(h2 == d_max) h2 = 0;
1251  if(h3 == d_max) h3 = 0;
1252  if(he2 == d_max) he2 = 0;
1253  if(he3 == d_max) he3 = 0;
1254 
1255  FILE * file;
1256  file = fopen(output_file.c_str(),"a");
1257  fprintf(file, "nA = %d;\n", eq_data_->dh_cr_disc_->distr()->size());
1258  fprintf(file, "nB = %d;\n", eq_data_->dh_->mesh()->get_el_ds()->size());
1259  fprintf(file, "nBF = %d;\n", eq_data_->dh_cr_->distr()->size());
1260  fprintf(file, "h1 = %e;\nh2 = %e;\nh3 = %e;\n", h1, h2, h3);
1261  fprintf(file, "he2 = %e;\nhe3 = %e;\n", he2, he3);
1262  fclose(file);
1263 }
1264 
1265 
1266 //template <int dim>
1267 //std::vector<arma::vec3> dof_points(DHCellAccessor cell, const Mapping<dim, 3> &mapping) {
1268 //
1269 //
1270 // vector<arma::vec::fixed<dim+1>> bary_dof_points = cell->fe()->dof_points();
1271 //
1272 // std::vector<arma::vec3> points(20);
1273 // points.resize(0);
1274 //
1275 //}
1276 //
1277 
1279  START_TIMER("DarcyFlowMH_Steady::set_mesh_data_for_bddc");
1280  // prepare mesh for BDDCML
1281  // initialize arrays
1282  // auxiliary map for creating coordinates of local dofs and global-to-local numbering
1283  std::map<int, arma::vec3> localDofMap;
1284  // connectivity for the subdomain, i.e. global dof numbers on element, stored element-by-element
1285  // Indices of Nodes on Elements
1286  std::vector<int> inet;
1287  // number of degrees of freedom on elements - determines elementwise chunks of INET array
1288  // Number of Nodes on Elements
1289  std::vector<int> nnet;
1290  // Indices of Subdomain Elements in Global Numbering - for local elements, their global indices
1291  std::vector<int> isegn;
1292 
1293  // This array is currently not used in BDDCML, it was used as an interface scaling alternative to scaling
1294  // by diagonal. It corresponds to the rho-scaling.
1295  std::vector<double> element_permeability;
1296 
1297  // maximal and minimal dimension of elements
1298  uint elDimMax = 1;
1299  uint elDimMin = 3;
1300  std::vector<LongIdx> cell_dofs_global(10);
1301 
1302 
1303 
1304  for ( DHCellAccessor dh_cell : eq_data_->dh_->own_range() ) {
1305  // for each element, create local numbering of dofs as fluxes (sides), pressure (element centre), Lagrange multipliers (edges), compatible connections
1306 
1307  dh_cell.get_dof_indices(cell_dofs_global);
1308 
1309  inet.insert(inet.end(), cell_dofs_global.begin(), cell_dofs_global.end());
1310  uint n_inet = cell_dofs_global.size();
1311 
1312 
1313  uint dim = dh_cell.elm().dim();
1314  elDimMax = std::max( elDimMax, dim );
1315  elDimMin = std::min( elDimMin, dim );
1316 
1317  // TODO: this is consistent with previous implementation, but may be wrong as it use global element numbering
1318  // used in sequential mesh, do global numbering of distributed elements.
1319  isegn.push_back( dh_cell.elm_idx());
1320 
1321  // TODO: use FiniteElement::dof_points
1322  for (unsigned int si=0; si<dh_cell.elm()->n_sides(); si++) {
1323  arma::vec3 coord = dh_cell.elm().side(si)->centre();
1324  // TODO: replace with DHCell getter when available for FESystem component
1325  // flux dof points
1326  localDofMap.insert( std::make_pair( cell_dofs_global[si], coord ) );
1327  // pressure trace dof points
1328  localDofMap.insert( std::make_pair( cell_dofs_global[si+dim+2], coord ) );
1329  }
1330  // pressure dof points
1331  arma::vec3 elm_centre = dh_cell.elm().centre();
1332  localDofMap.insert( std::make_pair( cell_dofs_global[dim+1], elm_centre ) );
1333 
1334  // insert dofs related to compatible connections
1335  //const Element *ele = dh_cell.elm().element();
1336  for(DHCellSide side : dh_cell.neighb_sides()) {
1337  uint neigh_dim = side.cell().elm().dim();
1338  side.cell().get_dof_indices(cell_dofs_global);
1339  int edge_row = cell_dofs_global[neigh_dim+2+side.side_idx()];
1340  localDofMap.insert( std::make_pair( edge_row, side.centre() ) );
1341  inet.push_back( edge_row );
1342  n_inet++;
1343  }
1344  nnet.push_back(n_inet);
1345 
1346 
1347  // version for rho scaling
1348  // trace computation
1349  double conduct = eq_fields_->conductivity.value( elm_centre , dh_cell.elm() );
1350  auto aniso = eq_fields_->anisotropy.value( elm_centre , dh_cell.elm() );
1351 
1352  // compute mean on the diagonal
1353  double coef = 0.;
1354  for ( int i = 0; i < 3; i++) {
1355  coef = coef + aniso.at(i,i);
1356  }
1357  // Maybe divide by cs
1358  coef = conduct*coef / 3;
1359 
1360  ASSERT_GT( coef, 0. ).error("Zero coefficient of hydrodynamic resistance. \n");
1361  element_permeability.push_back( 1. / coef );
1362  }
1363 // uint i_inet = 0;
1364 // for(int n_dofs : nnet) {
1365 // DebugOut() << "nnet: " << n_dofs;
1366 // for(int j=0; j < n_dofs; j++, i_inet++) {
1367 // DebugOut() << "inet: " << inet[i_inet];
1368 // }
1369 // }
1370 
1371  auto distr = eq_data_->dh_->distr();
1372 // for(auto pair : localDofMap) {
1373 // DebugOut().every_proc() << "r: " << distr->myp() << " gi: " << pair.first << "xyz: " << pair.second[0];
1374 //
1375 // }
1376 
1377 
1378  //convert set of dofs to vectors
1379  // number of nodes (= dofs) on the subdomain
1380  int numNodeSub = localDofMap.size();
1381  //ASSERT_PERMANENT_EQ( (unsigned int)numNodeSub, eq:data_->dh_->lsize() );
1382  // Indices of Subdomain Nodes in Global Numbering - for local nodes, their global indices
1383  std::vector<int> isngn( numNodeSub );
1384  // pseudo-coordinates of local nodes (i.e. dofs)
1385  // they need not be exact, they are used just for some geometrical considerations in BDDCML,
1386  // such as selection of corners maximizing area of a triangle, bounding boxes fro subdomains to
1387  // find candidate neighbours etc.
1388  std::vector<double> xyz( numNodeSub * 3 ) ;
1389  int ind = 0;
1390  std::map<int,arma::vec3>::iterator itB = localDofMap.begin();
1391  for ( ; itB != localDofMap.end(); ++itB ) {
1392  isngn[ind] = itB -> first;
1393 
1394  arma::vec3 coord = itB -> second;
1395  for ( int j = 0; j < 3; j++ ) {
1396  xyz[ j*numNodeSub + ind ] = coord[j];
1397  }
1398 
1399  ind++;
1400  }
1401  localDofMap.clear();
1402 
1403  // Number of Nodal Degrees of Freedom
1404  // nndf is trivially one - dofs coincide with nodes
1405  std::vector<int> nndf( numNodeSub, 1 );
1406 
1407  // prepare auxiliary map for renumbering nodes
1408  typedef std::map<int,int> Global2LocalMap_; //! type for storage of global to local map
1409  Global2LocalMap_ global2LocalNodeMap;
1410  for ( unsigned ind = 0; ind < isngn.size(); ++ind ) {
1411  global2LocalNodeMap.insert( std::make_pair( static_cast<unsigned>( isngn[ind] ), ind ) );
1412  }
1413 
1414  // renumber nodes in the inet array to locals
1415  int indInet = 0;
1416  for ( unsigned int iEle = 0; iEle < isegn.size(); iEle++ ) {
1417  int nne = nnet[ iEle ];
1418  for ( int ien = 0; ien < nne; ien++ ) {
1419 
1420  int indGlob = inet[indInet];
1421  // map it to local node
1422  Global2LocalMap_::iterator pos = global2LocalNodeMap.find( indGlob );
1423  ASSERT( pos != global2LocalNodeMap.end() )(indGlob).error("Cannot remap node global index to local indices.\n");
1424  int indLoc = static_cast<int> ( pos -> second );
1425 
1426  // store the node
1427  inet[ indInet++ ] = indLoc;
1428  }
1429  }
1430 
1431  int numNodes = size;
1432  int numDofsInt = size;
1433  int spaceDim = 3; // TODO: what is the proper value here?
1434  int meshDim = elDimMax;
1435 
1436  /**
1437  * We need:
1438  * - local to global element map (possibly mesh->el_4_loc
1439  * - inet, nnet - local dof numbers per element, local numbering of only those dofs that are on owned elements
1440  * 1. collect DH local dof indices on elements, manage map from DH local indices to BDDC local dof indices
1441  * 2. map collected DH indices to BDDC indices using the map
1442  * - local BDDC dofs to global dofs, use DH to BDDC map with DH local to global map
1443  * - XYZ - permuted, collect in main loop into array of size of all DH local dofs, compress and rearrange latter
1444  * - element_permeability - in main loop
1445  */
1446  bddc_ls -> load_mesh( LinSys_BDDC::BDDCMatrixType::SPD_VIA_SYMMETRICGENERAL, spaceDim, numNodes, numDofsInt, inet, nnet, nndf, isegn, isngn, isngn, xyz, element_permeability, meshDim );
1447 }
1448 
1449 
1450 
1451 
1452 //=============================================================================
1453 // DESTROY WATER MH SYSTEM STRUCTURE
1454 //=============================================================================
1456  if (previous_solution != nullptr) chkerr(VecDestroy(&previous_solution));
1457  if (steady_diagonal != nullptr) chkerr(VecDestroy(&steady_diagonal));
1458  if (new_diagonal != nullptr) chkerr(VecDestroy(&new_diagonal));
1459  if (steady_rhs != nullptr) chkerr(VecDestroy(&steady_rhs));
1460 
1461 
1462  if (schur0 != NULL) {
1463  delete schur0;
1464  }
1465 
1466  if (output_object) delete output_object;
1467 
1468  if(time_ != nullptr)
1469  delete time_;
1470 
1471 }
1472 
1473 
1474 /*
1475 void mat_count_off_proc_values(Mat m, Vec v) {
1476  int n, first, last;
1477  const PetscInt *cols;
1478  Distribution distr(v);
1479 
1480  int n_off = 0;
1481  int n_on = 0;
1482  int n_off_rows = 0;
1483  MatGetOwnershipRange(m, &first, &last);
1484  for (int row = first; row < last; row++) {
1485  MatGetRow(m, row, &n, &cols, PETSC_NULL);
1486  bool exists_off = false;
1487  for (int i = 0; i < n; i++)
1488  if (distr.get_proc(cols[i]) != distr.myp())
1489  n_off++, exists_off = true;
1490  else
1491  n_on++;
1492  if (exists_off)
1493  n_off_rows++;
1494  MatRestoreRow(m, row, &n, &cols, PETSC_NULL);
1495  }
1496 }
1497 */
1498 
1499 
1500 
1501 
1502 
1503 
1504 
1505 
1506 
1507 
1508 
1509 
1511 {
1512  DebugOut().fmt("Setup with dt: {}\n", time_->dt());
1513  for ( DHCellAccessor dh_cell : eq_data_->dh_->own_range() ) {
1514  ElementAccessor<3> ele = dh_cell.elm();
1515 
1516  double init_value = eq_fields_->init_pressure.value(ele.centre(),ele);
1517 
1518  LocDofVec loc_dofs = dh_cell.get_loc_dof_indices();
1519  for (unsigned int i=0; i<ele->n_sides(); i++) {
1520  uint n_sides_of_edge = ele.side(i)->edge().n_sides();
1521  uint edge_dof = ele->n_sides()+1+i;
1522  // set initial condition
1523  eq_data_->full_solution.add(loc_dofs[edge_dof], init_value/n_sides_of_edge);
1524  }
1525  }
1526 
1527  eq_data_->full_solution.ghost_to_local_begin();
1528  eq_data_->full_solution.ghost_to_local_end();
1529  eq_data_->full_solution.local_to_ghost_begin();
1530  eq_data_->full_solution.local_to_ghost_end();
1531 }
1532 
1533 
1535 {
1536  START_TIMER("DarcyFlowMH::reconstruct_solution_from_schur");
1537 
1538  for ( DHCellAccessor dh_cell : eq_data_->dh_->own_range() ) {
1539  unsigned int dim = dh_cell.dim();
1540  assembler[dim-1]->assemble_reconstruct(dh_cell);
1541  }
1542 
1543  eq_data_->full_solution.local_to_ghost_begin();
1544  eq_data_->full_solution.local_to_ghost_end();
1545 }
1546 
1547 
1549  // save diagonal of steady matrix
1550  MatGetDiagonal(*( schur0->get_matrix() ), steady_diagonal);
1551  // save RHS
1552  VecCopy(*( schur0->get_rhs()), steady_rhs);
1553 
1554 
1555  PetscScalar *local_diagonal;
1556  VecGetArray(new_diagonal,& local_diagonal);
1557 
1558  DebugOut().fmt("Setup with dt: {}\n", time_->dt());
1559 
1560  balance_->start_mass_assembly(eq_data_->water_balance_idx);
1561 
1562  std::vector<LongIdx> dofs;
1563  dofs.reserve(eq_data_->dh_->max_elem_dofs());
1564 
1565  for ( DHCellAccessor dh_cell : eq_data_->dh_->own_range() ) {
1566  ElementAccessor<3> ele = dh_cell.elm();
1567  dofs.resize(dh_cell.n_dofs());
1568  dh_cell.get_dof_indices(dofs);
1569 
1570  // TODO: replace with DHCell getter when available for FESystem component
1571  const uint p_ele_dof = dh_cell.n_dofs() / 2;
1572  // set new diagonal
1573  double diagonal_coeff = eq_fields_->cross_section.value(ele.centre(), ele)
1574  * ( eq_fields_->storativity.value(ele.centre(), ele)
1575  +eq_fields_->extra_storativity.value(ele.centre(), ele))
1576  * ele.measure();
1577  local_diagonal[dh_cell.get_loc_dof_indices()[p_ele_dof]]= - diagonal_coeff / time_->dt();
1578 
1579  balance_->add_mass_values(eq_data_->water_balance_idx, dh_cell,
1580  {dh_cell.get_loc_dof_indices()[p_ele_dof]},
1581  {diagonal_coeff}, 0.0);
1582  }
1583  VecRestoreArray(new_diagonal,& local_diagonal);
1584  MatDiagonalSet(*( schur0->get_matrix() ), new_diagonal, ADD_VALUES);
1585 
1586  schur0->set_matrix_changed();
1587 
1588  balance_->finish_mass_assembly(eq_data_->water_balance_idx);
1589 }
1590 
1592  START_TIMER("modify system");
1593  if (time_->is_changed_dt() && time_->step().index()>0) {
1594  double scale_factor=time_->step(-2).length()/time_->step().length();
1595  if (scale_factor != 1.0) {
1596  // if time step has changed and setup_time_term not called
1597  MatDiagonalSet(*( schur0->get_matrix() ),steady_diagonal, INSERT_VALUES);
1598 
1599  VecScale(new_diagonal, time_->last_dt()/time_->dt());
1600  MatDiagonalSet(*( schur0->get_matrix() ),new_diagonal, ADD_VALUES);
1602  }
1603  }
1604 
1605  // modify RHS - add previous solution
1606  VecPointwiseMult(*( schur0->get_rhs()), new_diagonal, previous_solution);
1607 // VecPointwiseMult(*( schur0->get_rhs()), new_diagonal, schur0->get_solution());
1608  VecAXPY(*( schur0->get_rhs()), 1.0, steady_rhs);
1610 
1611  //VecSwap(previous_solution, schur0->get_solution());
1612 }
1613 
1614 
1615 /// Helper method fills range (min and max) of given component
1616 void dofs_range(unsigned int n_dofs, unsigned int &min, unsigned int &max, unsigned int component) {
1617  if (component==0) {
1618  min = 0;
1619  max = n_dofs/2;
1620  } else if (component==1) {
1621  min = n_dofs/2;
1622  max = (n_dofs+1)/2;
1623  } else {
1624  min = (n_dofs+1)/2;
1625  max = n_dofs;
1626  }
1627 }
1628 
1629 
1631  ASSERT_LT(component, 3).error("Invalid component!");
1632  unsigned int i, n_dofs, min, max;
1633  std::vector<int> dof_vec;
1634  std::vector<LongIdx> dof_indices(eq_data_->dh_->max_elem_dofs());
1635  for ( DHCellAccessor dh_cell : eq_data_->dh_->own_range() ) {
1636  n_dofs = dh_cell.get_dof_indices(dof_indices);
1637  dofs_range(n_dofs, min, max, component);
1638  for (i=min; i<max; ++i) dof_vec.push_back(dof_indices[i]);
1639  }
1640  return dof_vec;
1641 }
1642 
1643 
1644 //-----------------------------------------------------------------------------
1645 // vim: set cindent:
Functors of FieldModels used in Darcy flow module.
#define ASSERT(expr)
Definition: asserts.hh:351
#define ASSERT_PERMANENT(expr)
Allow use shorter versions of macro names if these names is not used with external library.
Definition: asserts.hh:348
#define ASSERT_LT(a, b)
Definition of comparative assert macro (Less Than) only for debug mode.
Definition: asserts.hh:301
#define ASSERT_EQ(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
Definition: asserts.hh:333
#define ASSERT_GT(a, b)
Definition of comparative assert macro (Greater Than) only for debug mode.
Definition: asserts.hh:317
static const Input::Type::Record & get_input_type()
Main balance input record type.
Definition: balance.cc:53
Cell accessor allow iterate over DOF handler cells.
const ElementAccessor< 3 > elm() const
Return ElementAccessor to element of loc_ele_idx_.
unsigned int n_dofs() const
Return number of dofs on given cell.
unsigned int get_dof_indices(std::vector< LongIdx > &indices) const
Fill vector of the global indices of dofs associated to the cell.
unsigned int dim() const
Return dimension of element appropriate to cell.
Side accessor allows to iterate over sides of DOF handler cell.
static Input::Type::Abstract & get_input_type()
MortarMethod
Type of experimental Mortar-like method for non-compatible 1d-2d interaction.
static const Input::Type::Instance & get_input_type_specific()
void output()
Calculate values for output.
static const Input::Type::Instance & get_input_type(FieldSet &eq_data, const std::string &equation_name)
EqData()
Constructor.
EqFields()
Creation of all fields.
static const Input::Type::Selection & get_bc_type_selection()
Return a Selection corresponding to enum BC_Type.
DarcyFlowMHOutput * output_object
unsigned int max_n_it_
void set_mesh_data_for_bddc(LinSys_BDDC *bddc_ls)
Vec previous_solution
virtual double solved_time() override
void prepare_new_time_step()
int n_schur_compls
void assembly_mh_matrix(MultidimAssembly &assembler)
static const Input::Type::Record & type_field_descriptor()
virtual void read_initial_condition()
void print_matlab_matrix(string matlab_file)
Print darcy flow matrix in matlab format into a file.
virtual void initialize_specific()
std::shared_ptr< Balance > balance_
double tolerance_
std::vector< int > get_component_indices_vec(unsigned int component) const
Get vector of all DOF indices of given component (0..side, 1..element, 2..edge)
void allocate_mh_matrix()
Vec new_diagonal
void initialize() override
LinSys * schur0
static std::string equation_name()
void solve_nonlinear()
Solve method common to zero_time_step and update solution.
void init_eq_data()
virtual ~DarcyMH() override
std::shared_ptr< EqData > eq_data_
bool use_steady_assembly_
void zero_time_step() override
Vec steady_rhs
std::shared_ptr< EqFields > eq_fields_
static const int registrar
Registrar of class to factory.
unsigned int min_n_it_
virtual double solution_precision() const
unsigned int nonlinear_iteration_
virtual void output_data() override
Write computed fields.
virtual void assembly_linear_system()
static const Input::Type::Selection & get_mh_mortar_selection()
Selection for enum MortarMethod.
friend class DarcyFlowMHOutput
void reconstruct_solution_from_schur(MultidimAssembly &assembler)
virtual bool zero_time_term(bool time_global=false)
Vec steady_diagonal
void modify_system()
EqFields & eq_fields()
DarcyMH(Mesh &mesh, const Input::Record in_rec, TimeGovernor *tm=nullptr)
CREATE AND FILL GLOBAL MH MATRIX OF THE WATER MODEL.
virtual void postprocess()
postprocess velocity field (add sources)
static const Input::Type::Record & get_input_type()
void update_solution() override
void create_linear_system(Input::AbstractRecord rec)
bool data_changed_
void solve_time_step(bool output=true)
Solve the problem without moving to next time and without output.
virtual void setup_time_term()
virtual void assembly_source_term()
Source term is implemented differently in LMH version.
unsigned int begin(int proc) const
get starting local index
unsigned int n_sides() const
Returns number of sides aligned with the edge.
Definition: accessors.hh:329
double measure() const
Computes the measure of the element.
SideIter side(const unsigned int loc_index)
Region region() const
Definition: accessors.hh:198
unsigned int idx() const
We need this method after replacing Region by RegionIdx, and movinf RegionDB instance into particular...
Definition: accessors.hh:215
arma::vec::fixed< spacedim > centre() const
Computes the barycenter.
unsigned int n_neighs_vb() const
Return number of neighbours.
Definition: elements.h:65
unsigned int n_sides() const
Definition: elements.h:131
Input::Record input_record_
Definition: equation.hh:242
static Input::Type::Record & record_template()
Template Record with common keys for derived equations.
Definition: equation.cc:39
std::shared_ptr< FieldSet > eq_fieldset_
Definition: equation.hh:249
TimeGovernor * time_
Definition: equation.hh:241
Mesh * mesh_
Definition: equation.hh:240
TimeGovernor & time()
Definition: equation.hh:151
Compound finite element on dim dimensional simplex.
Definition: fe_system.hh:102
static const std::string field_descriptor_record_description(const string &record_name)
Definition: field_common.cc:72
static constexpr Mask input_copy
Definition: field_flag.hh:44
static constexpr Mask equation_result
Match result fields. These are never given by input or copy of input.
Definition: field_flag.hh:55
Dedicated class for storing path to input and output files.
Definition: file_path.hh:54
@ output_file
Definition: file_path.hh:69
Accessor to the polymorphic input data of a type given by an AbstracRecord object.
Definition: accessors.hh:458
Input::Type::Record type() const
Definition: accessors.cc:273
Accessor to input data conforming to declared Array.
Definition: accessors.hh:566
Accessor to the data with type Type::Record.
Definition: accessors.hh:291
const Ret val(const string &key) const
Class for declaration of inputs sequences.
Definition: type_base.hh:339
Class for declaration of the input of type Bool.
Definition: type_base.hh:452
Class Input::Type::Default specifies default value of keys of a Input::Type::Record.
Definition: type_record.hh:61
static Default obligatory()
The factory function to make an empty default value which is obligatory.
Definition: type_record.hh:110
static Default optional()
The factory function to make an empty default value which is optional.
Definition: type_record.hh:124
Class for declaration of the input data that are floating point numbers.
Definition: type_base.hh:534
Class for declaration of the integral input data.
Definition: type_base.hh:483
Record type proxy class.
Definition: type_record.hh:182
unsigned int size() const
Returns number of keys in the Record.
Definition: type_record.hh:602
virtual Record & derive_from(Abstract &parent)
Method to derive new Record from an AbstractRecord parent.
Definition: type_record.cc:196
Record & close() const
Close the Record for further declarations of keys.
Definition: type_record.cc:304
Record & copy_keys(const Record &other)
Copy keys from other record.
Definition: type_record.cc:216
Record & declare_key(const string &key, std::shared_ptr< TypeBase > type, const Default &default_value, const string &description, TypeBase::attribute_map key_attributes=TypeBase::attribute_map())
Declares a new key of the Record.
Definition: type_record.cc:503
Template for classes storing finite set of named values.
Selection & add_value(const int value, const std::string &key, const std::string &description="", TypeBase::attribute_map attributes=TypeBase::attribute_map())
Adds one new value with name given by key to the Selection.
const Selection & close() const
Close the Selection, no more values can be added.
Common base for intersection object.
unsigned int bulk_ele_idx() const
Returns index of bulk element.
static const Input::Type::Record & get_input_type()
Definition: linsys_BDDC.cc:39
void set_from_input(const Input::Record in_rec) override
Definition: linsys_BDDC.cc:279
void set_from_input(const Input::Record in_rec) override
static const Input::Type::Record & get_input_type()
Definition: linsys_PETSC.cc:32
const Vec & get_solution()
Definition: linsys.hh:282
virtual const Vec * get_rhs()
Definition: linsys.hh:203
void set_negative_definite(bool flag=true)
Definition: linsys.hh:588
void set_solution(Vec sol_vec)
Definition: linsys.hh:290
void set_rhs_changed()
Definition: linsys.hh:218
void set_matrix_changed()
Definition: linsys.hh:212
virtual double get_solution_precision()=0
virtual void start_add_assembly()
Definition: linsys.hh:341
virtual void finish_assembly()=0
void rhs_set_value(int row, double val)
Definition: linsys.hh:381
virtual const Mat * get_matrix()
Definition: linsys.hh:187
virtual void set_tolerances(double r_tol, double a_tol, double d_tol, unsigned int max_it)=0
virtual void mat_set_values(int nrow, int *rows, int ncol, int *cols, double *vals)=0
virtual double compute_residual()=0
void set_symmetric(bool flag=true)
Definition: linsys.hh:561
virtual PetscErrorCode rhs_zero_entries()
Definition: linsys.hh:273
virtual void start_allocation()
Definition: linsys.hh:333
virtual SolveInfo solve()=0
void set_positive_definite(bool flag=true)
Definition: linsys.hh:576
virtual PetscErrorCode mat_zero_entries()
Definition: linsys.hh:264
static Input::Type::Abstract & get_input_type()
Definition: linsys.cc:29
const RegionDB & region_db() const
Definition: mesh.h:175
unsigned int n_edges() const
Definition: mesh.h:114
unsigned int n_elements() const
Definition: mesh.h:111
Range< ElementAccessor< 3 > > elements_range() const
Returns range of mesh elements.
Definition: mesh.cc:1174
Definition: mesh.h:362
BCMesh * bc_mesh() const override
Implement MeshBase::bc_mesh(), getter of boundary mesh.
Definition: mesh.h:567
MixedMeshIntersections & mixed_intersections()
Definition: mesh.cc:849
unsigned int n_sides() const
Definition: mesh.cc:308
std::vector< std::vector< ILpair > > element_intersections_
static const Input::Type::Record & get_input_type()
The specification of output stream.
Definition: output_time.cc:38
RegionSet get_region_set(const std::string &set_name) const
Definition: region.cc:328
unsigned int bulk_idx() const
Returns index of the region in the bulk set.
Definition: region.hh:90
Distribution * make_complement_distribution()
get distribution of complement object if complement is defined
Definition: schur.cc:294
void set_from_input(const Input::Record in_rec) override
Definition: schur.cc:138
void set_complement(LinSys_PETSC *ls)
Set complement LinSys object.
Definition: schur.cc:287
int loc_size_B
Definition: schur.hh:137
Edge edge() const
Returns pointer to the edge connected to the side.
Basic time management functionality for unsteady (and steady) solvers (class Equation).
double dt() const
double t() const
bool is_end() const
Returns true if the actual time is greater than or equal to the end time.
int set_upper_constraint(double upper, std::string message)
Sets upper constraint for the next time step estimating.
bool is_changed_dt() const
void view(const char *name="") const
double last_dt() const
double estimate_dt() const
Estimate choice of next time step according to actual setting of constraints.
const TimeStep & step(int index=-1) const
void next_time()
Proceed to the next time according to current estimated time step.
double length() const
unsigned int index() const
Class for representation SI units of Fields.
Definition: unit_si.hh:40
static UnitSI & dimensionless()
Returns dimensionless unit.
Definition: unit_si.cc:55
void dofs_range(unsigned int n_dofs, unsigned int &min, unsigned int &max, unsigned int component)
Helper method fills range (min and max) of given component.
Output class for darcy_flow_mh model.
Support classes for parallel programing.
Definitions of basic Lagrangean finite elements with polynomial shape functions.
@ result_zeros
#define FLOW123D_FORCE_LINK_IN_CHILD(x)
Definition: global_defs.h:104
#define THROW(whole_exception_expr)
Wrapper for throw. Saves the throwing point.
Definition: exceptions.hh:53
arma::Col< IntIdx > LocDofVec
Definition: index_types.hh:28
int LongIdx
Define type that represents indices of large arrays (elements, nodes, dofs etc.)
Definition: index_types.hh:24
Classes with algorithms for computation of intersections of meshes.
Wrappers for linear systems based on MPIAIJ and MATIS format.
Solver based on Multilevel BDDC - using corresponding class of OpenFTL package.
Solver based on the original PETSc solver using MPIAIJ matrix and succesive Schur complement construc...
#define WarningOut()
Macro defining 'warning' record of log.
Definition: logger.hh:278
#define MessageOut()
Macro defining 'message' record of log.
Definition: logger.hh:275
#define DebugOut()
Macro defining 'debug' record of log.
Definition: logger.hh:284
unsigned int uint
#define MPI_INT
Definition: mpi.h:160
#define MPI_Allreduce(sendbuf, recvbuf, count, datatype, op, comm)
Definition: mpi.h:612
#define MPI_MIN
Definition: mpi.h:198
ArmaVec< double, N > vec
Definition: armor.hh:885
FMT_FUNC int fprintf(std::ostream &os, CStringRef format, ArgList args)
Definition: ostream.cc:56
Implementation of range helper class.
Assembly explicit Schur complement for the given linear system. Provides method for resolution of the...
SchurComplement SchurComplement
Definition: linsys.hh:91
int converged_reason
Definition: linsys.hh:108
#define END_TIMER(tag)
Ends a timer with specified tag.
#define START_TIMER(tag)
Starts a timer with specified tag.
void chkerr(unsigned int ierr)
Replacement of new/delete operator in the spirit of xmalloc.
Definition: system.hh:142
Basic time management class.