Flow123d  master-1fea4ce
darcy_flow_lmh.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_lmh.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_lmh_old_.hh"
52 #include "flow/darcy_flow_lmh.hh"
54 #include "flow/assembly_lmh.hh"
55 #include "flow/assembly_models.hh"
56 
57 #include "tools/time_governor.hh"
59 #include "fields/field.hh"
60 #include "fields/field_values.hh"
62 #include "fields/field_fe.hh"
63 #include "fields/field_model.hh"
64 #include "fields/field_constant.hh"
65 
66 #include "coupling/balance.hh"
67 
70 
71 #include "fem/fe_p.hh"
72 
73 
74 FLOW123D_FORCE_LINK_IN_CHILD(darcy_flow_lmh)
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 
90 
91  const it::Record &field_descriptor =
93  .copy_keys( DarcyLMH::EqFields().make_field_descriptor_type(equation_name() + "_Data_aux") )
94  .declare_key("bc_piezo_head", FieldAlgorithmBase< 3, FieldValue<3>::Scalar >::get_input_type_instance(),
95  "Boundary piezometric head for BC types: dirichlet, robin, and river." )
96  .declare_key("bc_switch_piezo_head", FieldAlgorithmBase< 3, FieldValue<3>::Scalar >::get_input_type_instance(),
97  "Boundary switch piezometric head for BC types: seepage, river." )
98  .declare_key("init_piezo_head", FieldAlgorithmBase< 3, FieldValue<3>::Scalar >::get_input_type_instance(),
99  "Initial condition for the pressure given as the piezometric head." )
100  .close();
101  return field_descriptor;
102 }
103 
105  it::Record ns_rec = Input::Type::Record("NonlinearSolver", "Non-linear solver settings.")
106  .declare_key("linear_solver", LinSys::get_input_type(), it::Default("{}"),
107  "Linear solver for MH problem.")
108  .declare_key("tolerance", it::Double(0.0), it::Default("1E-6"),
109  "Residual tolerance.")
110  .declare_key("min_it", it::Integer(0), it::Default("1"),
111  "Minimum number of iterations (linear solutions) to use.\nThis is usefull if the convergence criteria "
112  "does not characterize your goal well enough so it converges prematurely, possibly even without a single linear solution."
113  "If greater then 'max_it' the value is set to 'max_it'.")
114  .declare_key("max_it", it::Integer(0), it::Default("100"),
115  "Maximum number of iterations (linear solutions) of the non-linear solver.")
116  .declare_key("converge_on_stagnation", it::Bool(), it::Default("false"),
117  "If a stagnation of the nonlinear solver is detected the solver stops. "
118  "A divergence is reported by default, forcing the end of the simulation. By setting this flag to 'true', the solver "
119  "ends with convergence success on stagnation, but it reports warning about it.")
120  .close();
121 
123 
124  return it::Record(equation_name(), "Lumped Mixed-Hybrid solver for saturated Darcy flow.")
128  .declare_key("gravity", it::Array(it::Double(), 3,3), it::Default("[ 0, 0, -1]"),
129  "Vector of the gravity force. Dimensionless.")
131  "Input data for Darcy flow model.")
132  .declare_key("nonlinear_solver", ns_rec, it::Default("{}"),
133  "Non-linear solver for MH problem.")
134  .declare_key("output_stream", OutputTime::get_input_type(), it::Default("{}"),
135  "Output stream settings.\n Specify file format, precision etc.")
136 
138  IT::Default("{ \"fields\": [ \"pressure_p0\", \"velocity_p0\" ] }"),
139  "Specification of output fields and output times.")
141  "Output settings specific to Darcy flow model.\n"
142  "Includes raw output and some experimental functionality.")
143  .declare_key("balance", Balance::get_input_type(), it::Default("{}"),
144  "Settings for computing mass balance.")
145  .declare_key("mortar_method", get_mh_mortar_selection(), it::Default("\"None\""),
146  "Method for coupling Darcy flow between dimensions on incompatible meshes. [Experimental]" )
147  .close();
148 }
149 
150 
151 const int DarcyLMH::registrar =
152  Input::register_class< DarcyLMH, Mesh &, const Input::Record >(equation_name()) +
154 
155 
156 
158 {
159  *this += field_ele_pressure.name("pressure_p0")
160  .units(UnitSI().m())
162  .description("Pressure solution - P0 interpolation.");
163 
164  *this += field_edge_pressure.name("pressure_edge")
165  .units(UnitSI().m())
167  .description("Pressure solution - Crouzeix-Raviart interpolation.");
168 
169  *this += field_ele_piezo_head.name("piezo_head_p0")
170  .units(UnitSI().m())
172  .description("Piezo head solution - P0 interpolation.");
173 
174  *this += field_ele_velocity.name("velocity_p0")
175  .units(UnitSI().m().s(-1))
177  .description("Velocity solution - P0 interpolation.");
178 
179  *this += flux.name("flux")
180  .units(UnitSI().m().s(-1))
182  .description("Darcy flow flux.");
183 
184  *this += anisotropy.name("anisotropy")
185  .description("Anisotropy of the conductivity tensor.")
186  .input_default("1.0")
188 
189  *this += cross_section.name("cross_section")
190  .description("Complement dimension parameter (cross section for 1D, thickness for 2D).")
191  .input_default("1.0")
192  .units( UnitSI().m(3).md() );
193 
194  *this += conductivity.name("conductivity")
195  .description("Isotropic conductivity scalar.")
196  .input_default("1.0")
197  .units( UnitSI().m().s(-1) )
198  .set_limits(0.0);
199 
200  *this += sigma.name("sigma")
201  .description("Transition coefficient between dimensions.")
202  .input_default("1.0")
204 
205  *this += water_source_density.name("water_source_density")
206  .description("Water source density.")
207  .input_default("0.0")
208  .units( UnitSI().s(-1) );
209 
210  *this += bc_type.name("bc_type")
211  .description("Boundary condition type.")
213  .input_default("\"none\"")
215 
216  *this += bc_pressure
218  .name("bc_pressure")
219  .description("Prescribed pressure value on the boundary. Used for all values of ``bc_type`` except ``none`` and ``seepage``. "
220  "See documentation of ``bc_type`` for exact meaning of ``bc_pressure`` in individual boundary condition types.")
221  .input_default("0.0")
222  .units( UnitSI().m() );
223 
224  *this += bc_flux
226  .name("bc_flux")
227  .description("Incoming water boundary flux. Used for bc_types : ``total_flux``, ``seepage``, ``river``.")
228  .input_default("0.0")
229  .units( UnitSI().m().s(-1) );
230 
231  *this += bc_robin_sigma
233  .name("bc_robin_sigma")
234  .description("Conductivity coefficient in the ``total_flux`` or the ``river`` boundary condition type.")
235  .input_default("0.0")
236  .units( UnitSI().s(-1) );
237 
238  *this += bc_switch_pressure
240  .name("bc_switch_pressure")
241  .description("Critical switch pressure for ``seepage`` and ``river`` boundary conditions.")
242  .input_default("0.0")
243  .units( UnitSI().m() );
244 
245 
246  //these are for unsteady
247  *this += init_pressure.name("init_pressure")
248  .description("Initial condition for pressure in time dependent problems.")
249  .input_default("0.0")
250  .units( UnitSI().m() );
251 
252  *this += storativity.name("storativity")
253  .description("Storativity (in time dependent problems).")
254  .input_default("0.0")
255  .units( UnitSI().m(-1) );
256 
257  *this += extra_storativity.name("extra_storativity")
258  .description("Storativity added from upstream equation.")
259  .units( UnitSI().m(-1) )
260  .input_default("0.0")
261  .flags( input_copy );
262 
263  *this += extra_source.name("extra_water_source_density")
264  .description("Water source density added from upstream equation.")
265  .input_default("0.0")
266  .units( UnitSI().s(-1) )
267  .flags( input_copy );
268 
269  *this += gravity_field.name("gravity")
270  .description("Gravity vector.")
271  .input_default("0.0")
273 
274  *this += bc_gravity.name("bc_gravity")
275  .description("Boundary gravity vector.")
276  .input_default("0.0")
278 
279  *this += init_piezo_head.name("init_piezo_head")
280  .units(UnitSI().m())
281  .input_default("0.0")
282  .description("Init piezo head.");
283 
284  *this += bc_piezo_head.name("bc_piezo_head")
285  .units(UnitSI().m())
286  .input_default("0.0")
287  .description("Boundary piezo head.");
288 
289  *this += bc_switch_piezo_head.name("bc_switch_piezo_head")
290  .units(UnitSI().m())
291  .input_default("0.0")
292  .description("Boundary switch piezo head.");
293 
294  *this += ref_pressure.name("ref_pressure")
295  .units(UnitSI().m())
296  .input_default("0.0")
298  .description("Precomputed pressure of l2 difference output.");
299 
300  *this += ref_velocity.name("ref_velocity")
301  .units(UnitSI().m().s(-1))
302  .input_default("0.0")
304  .description("Precomputed velocity of l2 difference output.");
305 
306  *this += ref_divergence.name("ref_divergence")
307  .units(UnitSI().m())
308  .input_default("0.0")
310  .description("Precomputed divergence of l2 difference output.");
311 
312  this->set_default_fieldset();
313  //time_term_fields = this->subset({"storativity"});
314  //main_matrix_fields = this->subset({"anisotropy", "conductivity", "cross_section", "sigma", "bc_type", "bc_robin_sigma"});
315  //rhs_fields = this->subset({"water_source_density", "bc_pressure", "bc_flux"});
316 }
317 
318 
319 
321  return it::Selection("Flow_Darcy_BC_Type")
322  .add_value(none, "none",
323  "Homogeneous Neumann boundary condition\n(zero normal flux over the boundary).")
324  .add_value(dirichlet, "dirichlet",
325  "Dirichlet boundary condition. "
326  "Specify the pressure head through the ``bc_pressure`` field "
327  "or the piezometric head through the ``bc_piezo_head`` field.")
328  .add_value(total_flux, "total_flux", "Flux boundary condition (combines Neumann and Robin type). "
329  "Water inflow equal to (($ \\delta_d(q_d^N + \\sigma_d (h_d^R - h_d) )$)). "
330  "Specify the water inflow by the ``bc_flux`` field, the transition coefficient by ``bc_robin_sigma`` "
331  "and the reference pressure head or piezometric head through ``bc_pressure`` or ``bc_piezo_head`` respectively.")
332  .add_value(seepage, "seepage",
333  "Seepage face boundary condition. Pressure and inflow bounded from above. Boundary with potential seepage flow "
334  "is described by the pair of inequalities: "
335  "(($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. "
336  "Caution: setting (($q_d^N$)) strictly negative "
337  "may lead to an ill posed problem since a positive outflow is enforced. "
338  "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."
339  )
340  .add_value(river, "river",
341  "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: "
342  "(( $ \\delta_d(q_d^N + \\sigma_d(H_d^D - H_d) )$)). For the water level under the bedrock, constant infiltration is used: "
343  "(( $ \\delta_d(q_d^N + \\sigma_d(H_d^D - H_d^S) )$)). Parameters: ``bc_pressure``, ``bc_switch_pressure``, "
344  " ``bc_sigma``, ``bc_flux``."
345  )
346  .close();
347 }
348 
349 
350 
352 {
353  mortar_method_=NoMortar;
354 }
355 
356 
358 {
359  auto size = dh_p_->get_local_to_global_map().size();
360  save_local_system_.resize(size);
361  bc_fluxes_reconstruted.resize(size);
362  loc_system_.resize(size);
363  postprocess_solution_.resize(size);
364 }
365 
366 
368 {
369  std::fill(save_local_system_.begin(), save_local_system_.end(), false);
370  std::fill(bc_fluxes_reconstruted.begin(), bc_fluxes_reconstruted.end(), false);
371 }
372 
373 
374 
375 
376 
377 
378 //=============================================================================
379 // CREATE AND FILL GLOBAL MH MATRIX OF THE WATER MODEL
380 // - do it in parallel:
381 // - initial distribution of elements, edges
382 //
383 /*! @brief CREATE AND FILL GLOBAL MH MATRIX OF THE WATER MODEL
384  *
385  * Parameters {Solver,NSchurs} number of performed Schur
386  * complements (0,1,2) for water flow MH-system
387  *
388  */
389 //=============================================================================
390 DarcyLMH::DarcyLMH(Mesh &mesh_in, const Input::Record in_rec, TimeGovernor *tm)
391 : DarcyFlowInterface(mesh_in, in_rec),
392  output_object(nullptr),
393  data_changed_(false),
394  read_init_cond_assembly_(nullptr),
395  mh_matrix_assembly_(nullptr),
397 {
398 
399  START_TIMER("Darcy constructor");
400  {
401  auto time_record = input_record_.val<Input::Record>("time");
402  if (tm == nullptr)
403  {
404  time_ = new TimeGovernor(time_record);
405  }
406  else
407  {
408  TimeGovernor tm_from_rec(time_record);
409  if (!tm_from_rec.is_default()) // is_default() == false when time record is present in input file
410  {
411  MessageOut() << "Duplicate key 'time', time in flow equation is already initialized from parent class!";
412  ASSERT_PERMANENT(false);
413  }
414  time_ = tm;
415  }
416  }
417 
418  eq_fields_ = make_shared<EqFields>();
419  eq_data_ = make_shared<EqData>();
420  this->eq_fieldset_ = eq_fields_;
421 
422  eq_fields_->set_mesh(*mesh_);
423 
424  eq_data_->is_linear=true;
425 
426  size = mesh_->n_elements() + mesh_->n_sides() + mesh_->n_edges();
427  eq_data_->mortar_method_= in_rec.val<MortarMethod>("mortar_method");
428  if (eq_data_->mortar_method_ != NoMortar) {
430  }
431 
432 
433  //side_ds->view( std::cout );
434  //el_ds->view( std::cout );
435  //edge_ds->view( std::cout );
436  //rows_ds->view( std::cout );
437 
438 }
439 
441 {
442  // DebugOut() << "t = " << time_->t() << " step_end " << time_->step().end() << "\n";
443  if(eq_data_->use_steady_assembly_)
444  {
445  // In steady case, the solution is computed with the data present at time t,
446  // and the steady state solution is valid until another change in data,
447  // which should correspond to time (t+dt).
448  // "The data change appears immediatly."
449  double next_t = time_->t() + time_->estimate_dt();
450  // DebugOut() << "STEADY next_t = " << next_t << "\n";
451  return next_t * (1 - 2*std::numeric_limits<double>::epsilon());
452  }
453  else
454  {
455  // In unsteady case, the solution is computed with the data present at time t,
456  // and the solution is valid at the time t+dt.
457  // "The data change does not appear immediatly, it is integrated over time interval dt."
458  // DebugOut() << "UNSTEADY\n";
459  return time_->t();
460  }
461 }
462 
464 //connecting data fields with mesh
465 {
466 
467  START_TIMER("Darcy data init");
468  eq_data_->mesh = mesh_;
469 
470  auto gravity_array = input_record_.val<Input::Array>("gravity");
471  std::vector<double> gvec;
472  gravity_array.copy_to(gvec);
473  gvec.push_back(0.0); // zero pressure shift
474  eq_data_->gravity_ = arma::vec(gvec);
475  eq_data_->gravity_vec_ = eq_data_->gravity_.subvec(0,2);
476 
477  FieldValue<3>::VectorFixed gvalue(eq_data_->gravity_vec_);
478  auto field_algo=std::make_shared<FieldConstant<3, FieldValue<3>::VectorFixed>>();
479  field_algo->set_value(gvalue);
480  eq_fields_->gravity_field.set(field_algo, 0.0);
481  eq_fields_->bc_gravity.set(field_algo, 0.0);
482 
483  eq_fields_->bc_pressure.add_factory(
484  std::make_shared<AddPotentialFactory<3, FieldValue<3>::Scalar> >
485  (eq_fields_->bc_gravity, eq_fields_->X(), eq_fields_->bc_piezo_head) );
486  eq_fields_->bc_switch_pressure.add_factory(
487  std::make_shared<AddPotentialFactory<3, FieldValue<3>::Scalar> >
488  (eq_fields_->bc_gravity, eq_fields_->X(), eq_fields_->bc_switch_piezo_head) );
489  eq_fields_->init_pressure.add_factory(
490  std::make_shared<AddPotentialFactory<3, FieldValue<3>::Scalar> >
491  (eq_fields_->gravity_field, eq_fields_->X(), eq_fields_->init_piezo_head) );
492 
493 
494  eq_fields_->set_input_list( this->input_record_.val<Input::Array>("input_fields"), *time_ );
495 
496  // Check that the time step was set for the transient simulation.
497  if (! zero_time_term(true) && time_->is_default() ) {
498  //THROW(ExcAssertMsg());
499  //THROW(ExcMissingTimeGovernor() << input_record_.ei_address());
500  MessageOut() << "Missing the key 'time', obligatory for the transient problems." << endl;
501  ASSERT_PERMANENT(false);
502  }
503 
504  eq_fields_->mark_input_times(*time_);
505 }
506 
508 
509  { // init DOF handler for pressure fields
510 // std::shared_ptr< FiniteElement<0> > fe0_rt = std::make_shared<FE_RT0_disc<0>>();
511  std::shared_ptr< FiniteElement<1> > fe1_rt = std::make_shared<FE_RT0_disc<1>>();
512  std::shared_ptr< FiniteElement<2> > fe2_rt = std::make_shared<FE_RT0_disc<2>>();
513  std::shared_ptr< FiniteElement<3> > fe3_rt = std::make_shared<FE_RT0_disc<3>>();
514  std::shared_ptr< FiniteElement<0> > fe0_disc = std::make_shared<FE_P_disc<0>>(0);
515  std::shared_ptr< FiniteElement<1> > fe1_disc = std::make_shared<FE_P_disc<1>>(0);
516  std::shared_ptr< FiniteElement<2> > fe2_disc = std::make_shared<FE_P_disc<2>>(0);
517  std::shared_ptr< FiniteElement<3> > fe3_disc = std::make_shared<FE_P_disc<3>>(0);
518  std::shared_ptr< FiniteElement<0> > fe0_cr = std::make_shared<FE_CR<0>>();
519  std::shared_ptr< FiniteElement<1> > fe1_cr = std::make_shared<FE_CR<1>>();
520  std::shared_ptr< FiniteElement<2> > fe2_cr = std::make_shared<FE_CR<2>>();
521  std::shared_ptr< FiniteElement<3> > fe3_cr = std::make_shared<FE_CR<3>>();
522 // static FiniteElement<0> fe0_sys = FE_P_disc<0>(0); //TODO fix and use solution with FESystem<0>( {fe0_rt, fe0_disc, fe0_cr} )
523  FESystem<0> fe0_sys( {fe0_disc, fe0_disc, fe0_cr} );
524  FESystem<1> fe1_sys( {fe1_rt, fe1_disc, fe1_cr} );
525  FESystem<2> fe2_sys( {fe2_rt, fe2_disc, fe2_cr} );
526  FESystem<3> fe3_sys( {fe3_rt, fe3_disc, fe3_cr} );
527  MixedPtr<FESystem> fe_sys( std::make_shared<FESystem<0>>(fe0_sys), std::make_shared<FESystem<1>>(fe1_sys),
528  std::make_shared<FESystem<2>>(fe2_sys), std::make_shared<FESystem<3>>(fe3_sys) );
529  std::shared_ptr<DiscreteSpace> ds = std::make_shared<EqualOrderDiscreteSpace>( mesh_, fe_sys);
530  eq_data_->dh_ = std::make_shared<DOFHandlerMultiDim>(*mesh_);
531  eq_data_->dh_->distribute_dofs(ds);
532  }
533 
534  init_eq_data();
536 
537  eq_fields_->add_coords_field();
538 
539  { // construct pressure, velocity and piezo head fields
540  uint rt_component = 0;
541  eq_data_->full_solution = eq_data_->dh_->create_vector();
542  auto ele_flux_ptr = create_field_fe<3, FieldValue<3>::VectorFixed>(eq_data_->dh_, &eq_data_->full_solution, rt_component);
543  eq_fields_->flux.set(ele_flux_ptr, 0.0);
544 
545  eq_fields_->field_ele_velocity.set(Model<3, FieldValue<3>::VectorFixed>::create(fn_mh_velocity(), eq_fields_->flux, eq_fields_->cross_section), 0.0);
546 
547  uint p_ele_component = 1;
548  auto ele_pressure_ptr = create_field_fe<3, FieldValue<3>::Scalar>(eq_data_->dh_, &eq_data_->full_solution, p_ele_component);
549  eq_fields_->field_ele_pressure.set(ele_pressure_ptr, 0.0);
550 
551  uint p_edge_component = 2;
552  auto edge_pressure_ptr = create_field_fe<3, FieldValue<3>::Scalar>(eq_data_->dh_, &eq_data_->full_solution, p_edge_component);
553  eq_fields_->field_edge_pressure.set(edge_pressure_ptr, 0.0);
554 
555  eq_fields_->field_ele_piezo_head.set(
556  Model<3, FieldValue<3>::Scalar>::create(fn_mh_piezohead(), eq_fields_->gravity_field, eq_fields_->X(), eq_fields_->field_ele_pressure),
557  0.0
558  );
559  }
560 
561  { // init DOF handlers represents element pressure DOFs
562  uint p_element_component = 1;
563  eq_data_->dh_p_ = std::make_shared<SubDOFHandlerMultiDim>(eq_data_->dh_,p_element_component);
564  }
565 
566  { // init DOF handlers represents edge DOFs
567  uint p_edge_component = 2;
568  eq_data_->dh_cr_ = std::make_shared<SubDOFHandlerMultiDim>(eq_data_->dh_,p_edge_component);
569  }
570 
571  { // init DOF handlers represents side DOFs
572  MixedPtr<FE_CR_disc> fe_cr_disc;
573  std::shared_ptr<DiscreteSpace> ds_cr_disc = std::make_shared<EqualOrderDiscreteSpace>( mesh_, fe_cr_disc);
574  eq_data_->dh_cr_disc_ = std::make_shared<DOFHandlerMultiDim>(*mesh_);
575  eq_data_->dh_cr_disc_->distribute_dofs(ds_cr_disc);
576  }
577 
578  eq_data_->init();
579 
580  // create solution vector for 2. Schur complement linear system
581 // p_edge_solution = new VectorMPI(eq_data_->dh_cr_->distr()->lsize());
582 // full_solution = new VectorMPI(eq_data_->dh_->distr()->lsize());
583  // this creates mpi vector from DoFHandler, including ghost values
584  eq_data_->p_edge_solution = eq_data_->dh_cr_->create_vector();
585  eq_data_->p_edge_solution_previous = eq_data_->dh_cr_->create_vector();
586  eq_data_->p_edge_solution_previous_time = eq_data_->dh_cr_->create_vector();
587 
588  // Initialize bc_switch_dirichlet to size of global boundary.
589  eq_data_->bc_switch_dirichlet.resize(mesh_->n_elements()+mesh_->bc_mesh()->n_elements(), 1);
590 
591 
592  eq_data_->nonlinear_iteration_=0;
594  .val<Input::Record>("nonlinear_solver")
595  .val<Input::AbstractRecord>("linear_solver");
596 
598 
599  // auxiliary set_time call since allocation assembly evaluates fields as well
602 
603 
604  // initialization of balance object
605  balance_ = std::make_shared<Balance>("water", mesh_);
606  balance_->init_from_input(input_record_.val<Input::Record>("balance"), time());
607  eq_data_->water_balance_idx = balance_->add_quantity("water_volume");
608  balance_->allocate(eq_data_->dh_, 1);
609  balance_->units(UnitSI().m(3));
610 
611  eq_data_->balance_ = this->balance_;
612 
613  this->initialize_asm();
614 }
615 
617 {
618  //eq_data_->multidim_assembler = AssemblyFlowBase::create< AssemblyLMH >(eq_fields_, eq_data_);
619 }
620 
621 //void DarcyLMH::read_initial_condition()
622 //{
623 // DebugOut().fmt("Read initial condition\n");
624 //
625 // for ( DHCellAccessor dh_cell : eq_data_->dh_->own_range() ) {
626 //
627 // LocDofVec p_indices = dh_cell.cell_with_other_dh(eq_data_->dh_p_.get()).get_loc_dof_indices();
628 // ASSERT_DBG(p_indices.n_elem == 1);
629 // LocDofVec l_indices = dh_cell.cell_with_other_dh(eq_data_->dh_cr_.get()).get_loc_dof_indices();
630 // ElementAccessor<3> ele = dh_cell.elm();
631 //
632 // // set initial condition
633 // double init_value = eq_fields_->init_pressure.value(ele.centre(),ele);
634 // unsigned int p_idx = eq_data_->dh_p_->parent_indices()[p_indices[0]];
635 // eq_data_->full_solution.set(p_idx, init_value);
636 //
637 // for (unsigned int i=0; i<ele->n_sides(); i++) {
638 // uint n_sides_of_edge = ele.side(i)->edge().n_sides();
639 // unsigned int l_idx = eq_data_->dh_cr_->parent_indices()[l_indices[i]];
640 // eq_data_->full_solution.add(l_idx, init_value/n_sides_of_edge);
641 //
642 // eq_data_->p_edge_solution.add(l_indices[i], init_value/n_sides_of_edge);
643 // }
644 // }
645 //
646 // initial_condition_postprocess();
647 //
648 // eq_data_->full_solution.ghost_to_local_begin();
649 // eq_data_->full_solution.ghost_to_local_end();
650 //
651 // eq_data_->p_edge_solution.ghost_to_local_begin();
652 // eq_data_->p_edge_solution.ghost_to_local_end();
653 // eq_data_->p_edge_solution_previous_time.copy_from(eq_data_->p_edge_solution);
654 //}
655 //
656 //void DarcyLMH::initial_condition_postprocess()
657 //{}
658 
660 {
661  START_TIMER("Darcy zero time step");
662 
663  /* TODO:
664  * - Allow solution reconstruction (pressure and velocity) from initial condition on user request.
665  * - Steady solution as an intitial condition may be forced by setting inti_time =-1, and set data for the steady solver in that time.
666  * Solver should be able to switch from and to steady case depending on the zero time term.
667  */
668 
670 
671  // zero_time_term means steady case
672  eq_data_->use_steady_assembly_ = zero_time_term();
673 
674  eq_data_->p_edge_solution.zero_entries();
675 
676  if (eq_data_->use_steady_assembly_) { // steady case
677  MessageOut() << "Flow zero time step - steady case\n";
678  //read_initial_condition(); // Possible solution guess for steady case.
679  solve_nonlinear(); // with right limit data
680  } else {
681  MessageOut() << "Flow zero time step - unsteady case\n";
682  eq_data_->time_step_ = time_->dt();
684  this->read_init_cond_asm();
685  accept_time_step(); // accept zero time step, i.e. initial condition
686 
687 
688  // we reconstruct the initial solution here
689  // during the reconstruction assembly:
690  // - the balance objects are actually allocated
691  // - the full solution vector is computed
693  }
694  //solution_output(T,right_limit); // data for time T in any case
695  output_data();
696 
697  END_TIMER("Darcy zero time step");
698 }
699 
700 //=============================================================================
701 // COMPOSE and SOLVE WATER MH System possibly through Schur complements
702 //=============================================================================
704 {
705  START_TIMER("Darcy solve system");
706 
707  time_->next_time();
708 
709  time_->view("DARCY"); //time governor information output
710 
711  solve_time_step();
712 
713  eq_data_->full_solution.local_to_ghost_begin();
714  eq_data_->full_solution.local_to_ghost_end();
715 }
716 
717 void DarcyLMH::solve_time_step(bool output)
718 {
720  bool zero_time_term_from_left=zero_time_term();
721 
722  bool jump_time = eq_fields_->storativity.is_jump_time();
723  if (! zero_time_term_from_left) {
724  MessageOut() << "Flow time step - unsteady case\n";
725  // time term not treated as zero
726  // Unsteady solution up to the T.
727 
728  // this flag is necesssary for switching BC to avoid setting zero neumann on the whole boundary in the steady case
729  eq_data_->use_steady_assembly_ = false;
730 
731  solve_nonlinear(); // with left limit data
732  if(output)
734  if (jump_time) {
735  WarningOut() << "Output of solution discontinuous in time not supported yet.\n";
736  //solution_output(T, left_limit); // output use time T- delta*dt
737  //output_data();
738  }
739  }
740 
741  if (time_->is_end()) {
742  // output for unsteady case, end_time should not be the jump time
743  // but rether check that
744  if (! zero_time_term_from_left && ! jump_time && output)
745  output_data();
746  return;
747  }
748 
750  bool zero_time_term_from_right=zero_time_term();
751  if (zero_time_term_from_right) {
752  MessageOut() << "Flow time step - steady case\n";
753  // this flag is necesssary for switching BC to avoid setting zero neumann on the whole boundary in the steady case
754  eq_data_->use_steady_assembly_ = true;
755  solve_nonlinear(); // with right limit data
756  if(output)
758 
759  } else if (! zero_time_term_from_left && jump_time) {
760  WarningOut() << "Discontinuous time term not supported yet.\n";
761  //solution_transfer(); // internally call set_time(T, left) and set_time(T,right) again
762  //solve_nonlinear(); // with right limit data
763  }
764  //solution_output(T,right_limit); // data for time T in any case
765  if (output)
766  output_data();
767 }
768 
769 bool DarcyLMH::zero_time_term(bool time_global) {
770  if (time_global) {
771  return (eq_fields_->storativity.input_list_size() == 0);
772  } else {
773  return eq_fields_->storativity.field_result(mesh_->region_db().get_region_set("BULK")) == result_zeros;
774  }
775 }
776 
777 
779 {
780  START_TIMER("Darcy solve_nonlinear");
782  double residual_norm = lin_sys_schur().compute_residual();
783  eq_data_->nonlinear_iteration_ = 0;
784  MessageOut().fmt("[nonlinear solver] norm of initial residual: {}\n", residual_norm);
785 
786  // Reduce is_linear flag.
787  int is_linear_common;
788  MPI_Allreduce(&(eq_data_->is_linear), &is_linear_common,1, MPI_INT ,MPI_MIN,PETSC_COMM_WORLD);
789 
790  Input::Record nl_solver_rec = input_record_.val<Input::Record>("nonlinear_solver");
791  this->tolerance_ = nl_solver_rec.val<double>("tolerance");
792  this->max_n_it_ = nl_solver_rec.val<unsigned int>("max_it");
793  this->min_n_it_ = nl_solver_rec.val<unsigned int>("min_it");
794  if (this->min_n_it_ > this->max_n_it_) this->min_n_it_ = this->max_n_it_;
795 
796  if (! is_linear_common) {
797  // set tolerances of the linear solver unless they are set by user.
798  lin_sys_schur().set_tolerances(0.1*this->tolerance_, 0.01*this->tolerance_, 10000, 100);
799  }
800  vector<double> convergence_history;
801 
802  while (eq_data_->nonlinear_iteration_ < this->min_n_it_ ||
803  (residual_norm > this->tolerance_ && eq_data_->nonlinear_iteration_ < this->max_n_it_ )) {
804  ASSERT_EQ( convergence_history.size(), eq_data_->nonlinear_iteration_ );
805  convergence_history.push_back(residual_norm);
806 
807  // print_matlab_matrix("matrix_" + std::to_string(time_->step().index()) + "_it_" + std::to_string(nonlinear_iteration_));
808  // stagnation test
809  if (convergence_history.size() >= 5 &&
810  convergence_history[ convergence_history.size() - 1]/convergence_history[ convergence_history.size() - 2] > 0.9 &&
811  convergence_history[ convergence_history.size() - 1]/convergence_history[ convergence_history.size() - 5] > 0.8) {
812  // stagnation
813  if (input_record_.val<Input::Record>("nonlinear_solver").val<bool>("converge_on_stagnation")) {
814  WarningOut().fmt("Accept solution on stagnation. Its: {} Residual: {}\n", eq_data_->nonlinear_iteration_, residual_norm);
815  break;
816  } else {
817  THROW(ExcSolverDiverge() << EI_Reason("Stagnation."));
818  }
819  }
820 
821  if (! is_linear_common){
822  eq_data_->p_edge_solution_previous.copy_from(eq_data_->p_edge_solution);
823  eq_data_->p_edge_solution_previous.local_to_ghost_begin();
824  eq_data_->p_edge_solution_previous.local_to_ghost_end();
825  }
826 
828  MessageOut().fmt("[schur solver] lin. it: {}, reason: {}, residual: {}\n",
829  si.n_iterations, si.converged_reason, lin_sys_schur().compute_residual());
830 
831  eq_data_->nonlinear_iteration_++;
832 
833  // hack to make BDDC work with empty compute_residual
834  if (is_linear_common){
835  // we want to print this info in linear (and steady) case
836  residual_norm = lin_sys_schur().compute_residual();
837  MessageOut().fmt("[nonlinear solver] lin. it: {}, reason: {}, residual: {}\n",
838  si.n_iterations, si.converged_reason, residual_norm);
839  break;
840  }
841  data_changed_=true; // force reassembly for non-linear case
842 
843  double alpha = 1; // how much of new solution
844  VecAXPBY(eq_data_->p_edge_solution.petsc_vec(), (1-alpha), alpha, eq_data_->p_edge_solution_previous.petsc_vec());
845 
846  //LogOut().fmt("Linear solver ended with reason: {} \n", si.converged_reason );
847  //ASSERT_PERMANENT_GE( si.converged_reason, 0).error("Linear solver failed to converge.\n");
849 
850  residual_norm = lin_sys_schur().compute_residual();
851  MessageOut().fmt("[nonlinear solver] it: {} lin. it: {}, reason: {}, residual: {}\n",
852  eq_data_->nonlinear_iteration_, si.n_iterations, si.converged_reason, residual_norm);
853  }
854 
855 // reconstruct_solution_from_schur(eq_data_->multidim_assembler);
857 
858  // adapt timestep
859  if (! this->zero_time_term()) {
860  double mult = 1.0;
861  if (eq_data_->nonlinear_iteration_ < 3) mult = 1.6;
862  if (eq_data_->nonlinear_iteration_ > 7) mult = 0.7;
863  time_->set_upper_constraint(time_->dt() * mult, "Darcy adaptivity.");
864  // int result = time_->set_upper_constraint(time_->dt() * mult, "Darcy adaptivity.");
865  //DebugOut().fmt("time adaptivity, res: {} it: {} m: {} dt: {} edt: {}\n", result, nonlinear_iteration_, mult, time_->dt(), time_->estimate_dt());
866  }
867 }
868 
869 
871 {
872  eq_data_->p_edge_solution_previous_time.copy_from(eq_data_->p_edge_solution);
873  eq_data_->p_edge_solution_previous_time.local_to_ghost_begin();
874  eq_data_->p_edge_solution_previous_time.local_to_ghost_end();
875 }
876 
877 
879  START_TIMER("Darcy output data");
880 
881  // print_matlab_matrix("matrix_" + std::to_string(time_->step().index()));
882 
883  //time_->view("DARCY"); //time governor information output
884  this->output_object->output();
885 
886 
887  START_TIMER("Darcy balance output");
888  balance_->calculate_cumulative(eq_data_->water_balance_idx, eq_data_->full_solution.petsc_vec());
889  balance_->calculate_instant(eq_data_->water_balance_idx, eq_data_->full_solution.petsc_vec());
890  balance_->output();
891 }
892 
893 
894 //double DarcyLMH::solution_precision() const
895 //{
896 // return eq_data_->lin_sys_schur->get_solution_precision();
897 //}
898 
899 
900 // ===========================================================================================
901 //
902 // MATRIX ASSEMBLY - we use abstract assembly routine, where LS Mat/Vec SetValues
903 // are in fact pointers to allocating or filling functions - this is governed by Linsystem roitunes
904 //
905 // =======================================================================================
906 //void DarcyLMH::assembly_mh_matrix(FMT_UNUSED MultidimAssembly& assembler)
907 //{
908 // START_TIMER("DarcyLMH::assembly_steady_mh_matrix");
909 //
910 // // DebugOut() << "assembly_mh_matrix \n";
911 // // set auxiliary flag for switchting Dirichlet like BC
912 // eq_data_->force_no_neumann_bc = eq_data_->use_steady_assembly_ && (eq_data_->nonlinear_iteration_ == 0);
913 //
914 // balance_->start_flux_assembly(eq_data_->water_balance_idx);
915 // balance_->start_source_assembly(eq_data_->water_balance_idx);
916 // balance_->start_mass_assembly(eq_data_->water_balance_idx);
917 //
918 // // TODO: try to move this into balance, or have it in the generic assembler class, that should perform the cell loop
919 // // including various pre- and post-actions
920 //// for ( DHCellAccessor dh_cell : eq_data_->dh_->own_range() ) {
921 //// unsigned int dim = dh_cell.dim();
922 //// assembler[dim-1]->assemble(dh_cell);
923 //// }
924 // this->mh_matrix_assembly_->assemble(eq_data_->dh_);
925 //
926 //
927 // balance_->finish_mass_assembly(eq_data_->water_balance_idx);
928 // balance_->finish_source_assembly(eq_data_->water_balance_idx);
929 // balance_->finish_flux_assembly(eq_data_->water_balance_idx);
930 //
931 //}
932 
933 
935 {
936  START_TIMER("Darcy allocate_mh_matrix");
937 
938  // to make space for second schur complement, max. 10 neighbour edges of one el.
939  double zeros[100000];
940  for(int i=0; i<100000; i++) zeros[i] = 0.0;
941 
942  std::vector<LongIdx> tmp_rows;
943  tmp_rows.reserve(200);
944 
945  std::vector<LongIdx> dofs, dofs_ngh;
946  dofs.reserve(eq_data_->dh_cr_->max_elem_dofs());
947  dofs_ngh.reserve(eq_data_->dh_cr_->max_elem_dofs());
948 
949  // DebugOut() << "Allocate new schur\n";
950  for ( DHCellAccessor dh_cell : eq_data_->dh_cr_->own_range() ) {
951  ElementAccessor<3> ele = dh_cell.elm();
952 
953  const uint ndofs = dh_cell.n_dofs();
954  dofs.resize(dh_cell.n_dofs());
955  dh_cell.get_dof_indices(dofs);
956 
957  int* dofs_ptr = dofs.data();
958  lin_sys_schur().mat_set_values(ndofs, dofs_ptr, ndofs, dofs_ptr, zeros);
959 
960  tmp_rows.clear();
961 
962  // compatible neighborings rows
963  unsigned int n_neighs = ele->n_neighs_vb();
964  for ( DHCellSide neighb_side : dh_cell.neighb_sides() ) {
965  // every compatible connection adds a 2x2 matrix involving
966  // current element pressure and a connected edge pressure
967 
968  // read neighbor dofs (dh_cr dofhandler)
969  // neighbor cell owning neighb_side
970  DHCellAccessor dh_neighb_cell = neighb_side.cell();
971 
972  const uint ndofs_ngh = dh_neighb_cell.n_dofs();
973  dofs_ngh.resize(ndofs_ngh);
974  dh_neighb_cell.get_dof_indices(dofs_ngh);
975 
976  // local index of pedge dof on neighboring cell
977  tmp_rows.push_back(dofs_ngh[neighb_side.side().side_idx()]);
978  }
979 
980  lin_sys_schur().mat_set_values(ndofs, dofs_ptr, n_neighs, tmp_rows.data(), zeros); // (edges) x (neigh edges)
981  lin_sys_schur().mat_set_values(n_neighs, tmp_rows.data(), ndofs, dofs_ptr, zeros); // (neigh edges) x (edges)
982  lin_sys_schur().mat_set_values(n_neighs, tmp_rows.data(), n_neighs, tmp_rows.data(), zeros); // (neigh edges) x (neigh edges)
983 
984  tmp_rows.clear();
985 // if (eq_data_->mortar_method_ != NoMortar) {
986 // auto &isec_list = mesh_->mixed_intersections().element_intersections_[ele.idx()];
987 // for(auto &isec : isec_list ) {
988 // IntersectionLocalBase *local = isec.second;
989 // DHCellAccessor dh_cell_slave = eq_data_->dh_cr_->cell_accessor_from_element(local->bulk_ele_idx());
990 //
991 // const uint ndofs_slave = dh_cell_slave.n_dofs();
992 // dofs_ngh.resize(ndofs_slave);
993 // dh_cell_slave.get_dof_indices(dofs_ngh);
994 //
995 // //DebugOut().fmt("Alloc: {} {}", ele.idx(), local->bulk_ele_idx());
996 // for(unsigned int i_side=0; i_side < dh_cell_slave.elm()->n_sides(); i_side++) {
997 // tmp_rows.push_back( dofs_ngh[i_side] );
998 // //DebugOut() << "aedge" << print_var(tmp_rows[tmp_rows.size()-1]);
999 // }
1000 // }
1001 // }
1002 
1003  lin_sys_schur().mat_set_values(ndofs, dofs_ptr, tmp_rows.size(), tmp_rows.data(), zeros); // master edges x slave edges
1004  lin_sys_schur().mat_set_values(tmp_rows.size(), tmp_rows.data(), ndofs, dofs_ptr, zeros); // slave edges x master edges
1005  lin_sys_schur().mat_set_values(tmp_rows.size(), tmp_rows.data(), tmp_rows.size(), tmp_rows.data(), zeros); // slave edges x slave edges
1006  }
1007  // DebugOut() << "end Allocate new schur\n";
1008 
1009  // int local_dofs[10];
1010  // unsigned int nsides;
1011  // for ( DHCellAccessor dh_cell : eq_data_->dh_->own_range() ) {
1012  // LocalElementAccessorBase<3> ele_ac(dh_cell);
1013  // nsides = ele_ac.n_sides();
1014 
1015  // //allocate at once matrix [sides,ele,edges]x[sides,ele,edges]
1016  // loc_size = 1 + 2*nsides;
1017  // unsigned int i_side = 0;
1018 
1019  // for (; i_side < nsides; i_side++) {
1020  // local_dofs[i_side] = ele_ac.side_row(i_side);
1021  // local_dofs[i_side+nsides] = ele_ac.edge_row(i_side);
1022  // }
1023  // local_dofs[i_side+nsides] = ele_ac.ele_row();
1024  // int * edge_rows = local_dofs + nsides;
1025  // //int ele_row = local_dofs[0];
1026 
1027  // // whole local MH matrix
1028  // ls->mat_set_values(loc_size, local_dofs, loc_size, local_dofs, zeros);
1029 
1030 
1031  // // compatible neighborings rows
1032  // unsigned int n_neighs = ele_ac.element_accessor()->n_neighs_vb();
1033  // unsigned int i=0;
1034  // for ( DHCellSide neighb_side : dh_cell.neighb_sides() ) {
1035  // //for (unsigned int i = 0; i < n_neighs; i++) {
1036  // // every compatible connection adds a 2x2 matrix involving
1037  // // current element pressure and a connected edge pressure
1038  // Neighbour *ngh = ele_ac.element_accessor()->neigh_vb[i];
1039  // DHCellAccessor cell_higher_dim = eq_data_->dh_->cell_accessor_from_element(neighb_side.elem_idx());
1040  // LocalElementAccessorBase<3> acc_higher_dim( cell_higher_dim );
1041  // for (unsigned int j = 0; j < neighb_side.element().dim()+1; j++)
1042  // if (neighb_side.element()->edge_idx(j) == ngh->edge_idx()) {
1043  // int neigh_edge_row = acc_higher_dim.edge_row(j);
1044  // tmp_rows.push_back(neigh_edge_row);
1045  // break;
1046  // }
1047  // //DebugOut() << "CC" << print_var(tmp_rows[i]);
1048  // ++i;
1049  // }
1050 
1051  // // allocate always also for schur 2
1052  // ls->mat_set_values(nsides+1, edge_rows, n_neighs, tmp_rows.data(), zeros); // (edges, ele) x (neigh edges)
1053  // ls->mat_set_values(n_neighs, tmp_rows.data(), nsides+1, edge_rows, zeros); // (neigh edges) x (edges, ele)
1054  // ls->mat_set_values(n_neighs, tmp_rows.data(), n_neighs, tmp_rows.data(), zeros); // (neigh edges) x (neigh edges)
1055 
1056  // tmp_rows.clear();
1057 
1058  // if (eq_data_->mortar_method_ != NoMortar) {
1059  // auto &isec_list = mesh_->mixed_intersections().element_intersections_[ele_ac.ele_global_idx()];
1060  // for(auto &isec : isec_list ) {
1061  // IntersectionLocalBase *local = isec.second;
1062  // LocalElementAccessorBase<3> slave_acc( eq_data_->dh_->cell_accessor_from_element(local->bulk_ele_idx()) );
1063  // //DebugOut().fmt("Alloc: {} {}", ele_ac.ele_global_idx(), local->bulk_ele_idx());
1064  // for(unsigned int i_side=0; i_side < slave_acc.dim()+1; i_side++) {
1065  // tmp_rows.push_back( slave_acc.edge_row(i_side) );
1066  // //DebugOut() << "aedge" << print_var(tmp_rows[tmp_rows.size()-1]);
1067  // }
1068  // }
1069  // }
1070  // /*
1071  // for(unsigned int i_side=0; i_side < ele_ac.element_accessor()->n_sides(); i_side++) {
1072  // DebugOut() << "aedge:" << print_var(edge_rows[i_side]);
1073  // }*/
1074 
1075  // ls->mat_set_values(nsides, edge_rows, tmp_rows.size(), tmp_rows.data(), zeros); // master edges x neigh edges
1076  // ls->mat_set_values(tmp_rows.size(), tmp_rows.data(), nsides, edge_rows, zeros); // neigh edges x master edges
1077  // ls->mat_set_values(tmp_rows.size(), tmp_rows.data(), tmp_rows.size(), tmp_rows.data(), zeros); // neigh edges x neigh edges
1078 
1079  // }
1080 /*
1081  // alloc edge diagonal entries
1082  if(rank == 0)
1083  for( vector<Edge>::iterator edg = mesh_->edges.begin(); edg != mesh_->edges.end(); ++edg) {
1084  int edg_idx = mh_dh.row_4_edge[edg->side(0)->edge_idx()];
1085 
1086 // for( vector<Edge>::iterator edg2 = mesh_->edges.begin(); edg2 != mesh_->edges.end(); ++edg2){
1087 // int edg_idx2 = mh_dh.row_4_edge[edg2->side(0)->edge_idx()];
1088 // if(edg_idx == edg_idx2){
1089 // DBGCOUT(<< "P[ " << rank << " ] " << "edg alloc: " << edg_idx << " " << edg_idx2 << "\n");
1090  ls->mat_set_value(edg_idx, edg_idx, 0.0);
1091 // }
1092 // }
1093  }
1094  */
1095  /*
1096  if (mortar_method_ == MortarP0) {
1097  P0_CouplingAssembler(*this).assembly(*ls);
1098  } else if (mortar_method_ == MortarP1) {
1099  P1_CouplingAssembler(*this).assembly(*ls);
1100  }*/
1101 }
1102 
1103 
1104 
1105 /*******************************************************************************
1106  * COMPOSE WATER MH MATRIX WITHOUT SCHUR COMPLEMENT
1107  ******************************************************************************/
1108 
1110 
1111  START_TIMER("Darcy preallocation");
1112 
1113  // if (schur0 == NULL) { // create Linear System for MH matrix
1114 
1115 // if (in_rec.type() == LinSys_BDDC::get_input_type()) {
1116 // #ifdef FLOW123D_HAVE_BDDCML
1117 // WarningOut() << "For BDDC no Schur complements are used.";
1118 // n_schur_compls = 0;
1119 // LinSys_BDDC *ls = new LinSys_BDDC(&(*eq_data_->dh_->distr()),
1120 // true); // swap signs of matrix and rhs to make the matrix SPD
1121 // ls->set_from_input(in_rec);
1122 // ls->set_solution( eq_data_->full_solution.petsc_vec() );
1123 // // possible initialization particular to BDDC
1124 // START_TIMER("BDDC set mesh data");
1125 // set_mesh_data_for_bddc(ls);
1126 // schur0=ls;
1127 // END_TIMER("BDDC set mesh data");
1128 // #else
1129 // Exception
1130 // THROW( ExcBddcmlNotSupported() );
1131 // #endif // FLOW123D_HAVE_BDDCML
1132 // }
1133 // else
1134  if (in_rec.type() == LinSys_PETSC::get_input_type()) {
1135  // use PETSC for serial case even when user wants BDDC
1136 
1137  eq_data_->lin_sys_schur = std::make_shared<LinSys_PETSC>( &(*eq_data_->dh_cr_->distr()) );
1138  lin_sys_schur().set_from_input(in_rec);
1140  lin_sys_schur().set_solution( eq_data_->p_edge_solution.petsc_vec() );
1142  ((LinSys_PETSC *)&lin_sys_schur())->set_initial_guess_nonzero(true);
1143 
1144 // LinSys_PETSC *schur1, *schur2;
1145 
1146 // if (n_schur_compls == 0) {
1147 // LinSys_PETSC *ls = new LinSys_PETSC( &(*eq_data_->dh_->distr()) );
1148 
1149 // // temporary solution; we have to set precision also for sequantial case of BDDC
1150 // // final solution should be probably call of direct solver for oneproc case
1151 // // if (in_rec.type() != LinSys_BDDC::get_input_type()) ls->set_from_input(in_rec);
1152 // // else {
1153 // // ls->LinSys::set_from_input(in_rec); // get only common options
1154 // // }
1155 // ls->set_from_input(in_rec);
1156 
1157 // // ls->set_solution( eq_data_->full_solution.petsc_vec() );
1158 // schur0=ls;
1159 // } else {
1160 // IS is;
1161 // auto side_dofs_vec = get_component_indices_vec(0);
1162 
1163 // ISCreateGeneral(PETSC_COMM_SELF, side_dofs_vec.size(), &(side_dofs_vec[0]), PETSC_COPY_VALUES, &is);
1164 // //ISView(is, PETSC_VIEWER_STDOUT_SELF);
1165 // //ASSERT_PERMANENT(err == 0).error("Error in ISCreateStride.");
1166 
1167 // SchurComplement *ls = new SchurComplement(&(*eq_data_->dh_->distr()), is);
1168 
1169 // // make schur1
1170 // Distribution *ds = ls->make_complement_distribution();
1171 // if (n_schur_compls==1) {
1172 // schur1 = new LinSys_PETSC(ds);
1173 // schur1->set_positive_definite();
1174 // } else {
1175 // IS is;
1176 // auto elem_dofs_vec = get_component_indices_vec(1);
1177 
1178 // const PetscInt *b_indices;
1179 // ISGetIndices(ls->IsB, &b_indices);
1180 // uint b_size = ls->loc_size_B;
1181 // for(uint i_b=0, i_bb=0; i_b < b_size && i_bb < elem_dofs_vec.size(); i_b++) {
1182 // if (b_indices[i_b] == elem_dofs_vec[i_bb])
1183 // elem_dofs_vec[i_bb++] = i_b + ds->begin();
1184 // }
1185 // ISRestoreIndices(ls->IsB, &b_indices);
1186 
1187 
1188 // ISCreateGeneral(PETSC_COMM_SELF, elem_dofs_vec.size(), &(elem_dofs_vec[0]), PETSC_COPY_VALUES, &is);
1189 // //ISView(is, PETSC_VIEWER_STDOUT_SELF);
1190 // //ASSERT_PERMANENT(err == 0).error("Error in ISCreateStride.");
1191 // SchurComplement *ls1 = new SchurComplement(ds, is); // is is deallocated by SchurComplement
1192 // ls1->set_negative_definite();
1193 
1194 // // make schur2
1195 // schur2 = new LinSys_PETSC( ls1->make_complement_distribution() );
1196 // schur2->set_positive_definite();
1197 // ls1->set_complement( schur2 );
1198 // schur1 = ls1;
1199 // }
1200 // ls->set_complement( schur1 );
1201 // ls->set_from_input(in_rec);
1202 // // ls->set_solution( eq_data_->full_solution.petsc_vec() );
1203 // schur0=ls;
1204  // }
1205 
1206  START_TIMER("PETSc preallocation");
1208 
1210 
1211  eq_data_->full_solution.zero_entries();
1212  eq_data_->p_edge_solution.zero_entries();
1213  END_TIMER("PETSc preallocation");
1214  }
1215  else {
1216  THROW( ExcUnknownSolver() );
1217  }
1218 
1219  END_TIMER("Darcy preallocation");
1220 }
1221 
1223 {}
1224 
1225 //void DarcyLMH::reconstruct_solution_from_schur(MultidimAssembly& assembler)
1226 //{
1227 // START_TIMER("DarcyFlowMH::reconstruct_solution_from_schur");
1228 //
1229 // eq_data_->full_solution.zero_entries();
1230 // eq_data_->p_edge_solution.local_to_ghost_begin();
1231 // eq_data_->p_edge_solution.local_to_ghost_end();
1232 //
1233 // balance_->start_flux_assembly(eq_data_->water_balance_idx);
1234 // balance_->start_source_assembly(eq_data_->water_balance_idx);
1235 // balance_->start_mass_assembly(eq_data_->water_balance_idx);
1236 //
1237 // for ( DHCellAccessor dh_cell : eq_data_->dh_->own_range() ) {
1238 // unsigned int dim = dh_cell.dim();
1239 // assembler[dim-1]->assemble_reconstruct(dh_cell);
1240 // }
1241 //
1242 // eq_data_->full_solution.local_to_ghost_begin();
1243 // eq_data_->full_solution.local_to_ghost_end();
1244 //
1245 // balance_->finish_mass_assembly(eq_data_->water_balance_idx);
1246 // balance_->finish_source_assembly(eq_data_->water_balance_idx);
1247 // balance_->finish_flux_assembly(eq_data_->water_balance_idx);
1248 //}
1249 
1251  START_TIMER("Darcy assembly_linear_system");
1252 // DebugOut() << "DarcyLMH::assembly_linear_system\n";
1253 
1254  eq_data_->p_edge_solution.local_to_ghost_begin();
1255  eq_data_->p_edge_solution.local_to_ghost_end();
1256 
1257  eq_data_->is_linear=true;
1258  //DebugOut() << "Assembly linear system\n";
1259 // if (data_changed_) {
1260 // data_changed_ = false;
1261  {
1262  //DebugOut() << "Data changed\n";
1263  // currently we have no optimization for cases when just time term data or RHS data are changed
1264 // if (typeid(*schur0) != typeid(LinSys_BDDC)) {
1265 // schur0->start_add_assembly(); // finish allocation and create matrix
1266 // schur_compl->start_add_assembly();
1267 // }
1268 
1270 
1273 
1274  eq_data_->time_step_ = time_->dt();
1275 
1276  this->mh_matrix_assembly_->assemble(eq_data_->dh_);; // fill matrix
1277 // assembly_mh_matrix( eq_data_->multidim_assembler ); // fill matrix
1278 
1281 
1282  // print_matlab_matrix("matrix");
1283  }
1284 }
1285 
1286 
1287 void DarcyLMH::print_matlab_matrix(std::string matlab_file)
1288 {
1289  std::string output_file;
1290 
1291  // compute h_min for different dimensions
1292  double d_max = std::numeric_limits<double>::max();
1293  double h1 = d_max, h2 = d_max, h3 = d_max;
1294  double he2 = d_max, he3 = d_max;
1295  for (auto ele : mesh_->elements_range()) {
1296  switch(ele->dim()){
1297  case 1: h1 = std::min(h1,ele.measure()); break;
1298  case 2: h2 = std::min(h2,ele.measure()); break;
1299  case 3: h3 = std::min(h3,ele.measure()); break;
1300  }
1301 
1302  for (unsigned int j=0; j<ele->n_sides(); j++) {
1303  switch(ele->dim()){
1304  case 2: he2 = std::min(he2, ele.side(j)->measure()); break;
1305  case 3: he3 = std::min(he3, ele.side(j)->measure()); break;
1306  }
1307  }
1308  }
1309  if(h1 == d_max) h1 = 0;
1310  if(h2 == d_max) h2 = 0;
1311  if(h3 == d_max) h3 = 0;
1312  if(he2 == d_max) he2 = 0;
1313  if(he3 == d_max) he3 = 0;
1314 
1315  FILE * file;
1316  file = fopen(output_file.c_str(),"a");
1317  fprintf(file, "nA = %d;\n", eq_data_->dh_cr_disc_->distr()->size());
1318  fprintf(file, "nB = %d;\n", eq_data_->dh_->mesh()->get_el_ds()->size());
1319  fprintf(file, "nBF = %d;\n", eq_data_->dh_cr_->distr()->size());
1320  fprintf(file, "h1 = %e;\nh2 = %e;\nh3 = %e;\n", h1, h2, h3);
1321  fprintf(file, "he2 = %e;\nhe3 = %e;\n", he2, he3);
1322  fclose(file);
1323 
1324  {
1325  output_file = FilePath(matlab_file + "_sch_new.m", FilePath::output_file);
1326  PetscViewer viewer;
1327  PetscViewerASCIIOpen(PETSC_COMM_WORLD, output_file.c_str(), &viewer);
1328  PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
1329  MatView( *const_cast<Mat*>(lin_sys_schur().get_matrix()), viewer);
1330  VecView( *const_cast<Vec*>(lin_sys_schur().get_rhs()), viewer);
1331  VecView( *const_cast<Vec*>(&(lin_sys_schur().get_solution())), viewer);
1332  VecView( *const_cast<Vec*>(&(eq_data_->full_solution.petsc_vec())), viewer);
1333  }
1334 }
1335 
1336 
1337 //template <int dim>
1338 //std::vector<arma::vec3> dof_points(DHCellAccessor cell, const Mapping<dim, 3> &mapping) {
1339 //
1340 //
1341 // vector<arma::vec::fixed<dim+1>> bary_dof_points = cell->fe()->dof_points();
1342 //
1343 // std::vector<arma::vec3> points(20);
1344 // points.resize(0);
1345 //
1346 //}
1347 //
1348 
1349 // void DarcyLMH::set_mesh_data_for_bddc(LinSys_BDDC * bddc_ls) {
1350 // START_TIMER("DarcyFlowMH_Steady::set_mesh_data_for_bddc");
1351 // // prepare mesh for BDDCML
1352 // // initialize arrays
1353 // // auxiliary map for creating coordinates of local dofs and global-to-local numbering
1354 // std::map<int, arma::vec3> localDofMap;
1355 // // connectivity for the subdomain, i.e. global dof numbers on element, stored element-by-element
1356 // // Indices of Nodes on Elements
1357 // std::vector<int> inet;
1358 // // number of degrees of freedom on elements - determines elementwise chunks of INET array
1359 // // Number of Nodes on Elements
1360 // std::vector<int> nnet;
1361 // // Indices of Subdomain Elements in Global Numbering - for local elements, their global indices
1362 // std::vector<int> isegn;
1363 //
1364 // // This array is currently not used in BDDCML, it was used as an interface scaling alternative to scaling
1365 // // by diagonal. It corresponds to the rho-scaling.
1366 // std::vector<double> element_permeability;
1367 //
1368 // // maximal and minimal dimension of elements
1369 // uint elDimMax = 1;
1370 // uint elDimMin = 3;
1371 // std::vector<LongIdx> cell_dofs_global(10);
1372 //
1373 //
1374 //
1375 // for ( DHCellAccessor dh_cell : eq_data_->dh_->own_range() ) {
1376 // // LocalElementAccessorBase<3> ele_ac(dh_cell);
1377 // // for each element, create local numbering of dofs as fluxes (sides), pressure (element centre), Lagrange multipliers (edges), compatible connections
1378 //
1379 // dh_cell.get_dof_indices(cell_dofs_global);
1380 //
1381 // inet.insert(inet.end(), cell_dofs_global.begin(), cell_dofs_global.end());
1382 // uint n_inet = cell_dofs_global.size();
1383 //
1384 //
1385 // uint dim = dh_cell.elm().dim();
1386 // elDimMax = std::max( elDimMax, dim );
1387 // elDimMin = std::min( elDimMin, dim );
1388 //
1389 // // TODO: this is consistent with previous implementation, but may be wrong as it use global element numbering
1390 // // used in sequential mesh, do global numbering of distributed elements.
1391 // isegn.push_back( dh_cell.elm_idx());
1392 //
1393 // // TODO: use FiniteElement::dof_points
1394 // for (unsigned int si=0; si<dh_cell.elm()->n_sides(); si++) {
1395 // arma::vec3 coord = dh_cell.elm().side(si)->centre();
1396 // // flux dof points
1397 // localDofMap.insert( std::make_pair( cell_dofs_global[si], coord ) );
1398 // // pressure trace dof points
1399 // localDofMap.insert( std::make_pair( cell_dofs_global[si+dim+2], coord ) );
1400 // }
1401 // // pressure dof points
1402 // arma::vec3 elm_centre = dh_cell.elm().centre();
1403 // localDofMap.insert( std::make_pair( cell_dofs_global[dim+1], elm_centre ) );
1404 //
1405 // // insert dofs related to compatible connections
1406 // //const Element *ele = dh_cell.elm().element();
1407 // for(DHCellSide side : dh_cell.neighb_sides()) {
1408 // uint neigh_dim = side.cell().elm().dim();
1409 // side.cell().get_dof_indices(cell_dofs_global);
1410 // int edge_row = cell_dofs_global[neigh_dim+2+side.side_idx()];
1411 // localDofMap.insert( std::make_pair( edge_row, side.centre() ) );
1412 // inet.push_back( edge_row );
1413 // n_inet++;
1414 // }
1415 // nnet.push_back(n_inet);
1416 //
1417 //
1418 // // version for rho scaling
1419 // // trace computation
1420 // double conduct = eq_fields_->conductivity.value( elm_centre , dh_cell.elm() );
1421 // auto aniso = eq_fields_->anisotropy.value( elm_centre , dh_cell.elm() );
1422 //
1423 // // compute mean on the diagonal
1424 // double coef = 0.;
1425 // for ( int i = 0; i < 3; i++) {
1426 // coef = coef + aniso.at(i,i);
1427 // }
1428 // // Maybe divide by cs
1429 // coef = conduct*coef / 3;
1430 //
1431 // ASSERT_PERMANENT_GT(coef, 0).error("Zero coefficient of hydrodynamic resistance.\n");
1432 // element_permeability.push_back( 1. / coef );
1433 // }
1434 // // uint i_inet = 0;
1435 // // for(int n_dofs : nnet) {
1436 // // DebugOut() << "nnet: " << n_dofs;
1437 // // for(int j=0; j < n_dofs; j++, i_inet++) {
1438 // // DebugOut() << "inet: " << inet[i_inet];
1439 // // }
1440 // // }
1441 //
1442 // auto distr = eq_data_->dh_->distr();
1443 // // for(auto pair : localDofMap) {
1444 // // DebugOut().every_proc() << "r: " << distr->myp() << " gi: " << pair.first << "xyz: " << pair.second[0];
1445 // //
1446 // // }
1447 //
1448 //
1449 // //convert set of dofs to vectors
1450 // // number of nodes (= dofs) on the subdomain
1451 // int numNodeSub = localDofMap.size();
1452 // //ASSERT_PERMANENT_EQ( (unsigned int)numNodeSub, eq_data_->dh_->lsize() );
1453 // // Indices of Subdomain Nodes in Global Numbering - for local nodes, their global indices
1454 // std::vector<int> isngn( numNodeSub );
1455 // // pseudo-coordinates of local nodes (i.e. dofs)
1456 // // they need not be exact, they are used just for some geometrical considerations in BDDCML,
1457 // // such as selection of corners maximizing area of a triangle, bounding boxes fro subdomains to
1458 // // find candidate neighbours etc.
1459 // std::vector<double> xyz( numNodeSub * 3 ) ;
1460 // int ind = 0;
1461 // std::map<int,arma::vec3>::iterator itB = localDofMap.begin();
1462 // for ( ; itB != localDofMap.end(); ++itB ) {
1463 // isngn[ind] = itB -> first;
1464 //
1465 // arma::vec3 coord = itB -> second;
1466 // for ( int j = 0; j < 3; j++ ) {
1467 // xyz[ j*numNodeSub + ind ] = coord[j];
1468 // }
1469 //
1470 // ind++;
1471 // }
1472 // localDofMap.clear();
1473 //
1474 // // Number of Nodal Degrees of Freedom
1475 // // nndf is trivially one - dofs coincide with nodes
1476 // std::vector<int> nndf( numNodeSub, 1 );
1477 //
1478 // // prepare auxiliary map for renumbering nodes
1479 // typedef std::map<int,int> Global2LocalMap_; //! type for storage of global to local map
1480 // Global2LocalMap_ global2LocalNodeMap;
1481 // for ( unsigned ind = 0; ind < isngn.size(); ++ind ) {
1482 // global2LocalNodeMap.insert( std::make_pair( static_cast<unsigned>( isngn[ind] ), ind ) );
1483 // }
1484 //
1485 // // renumber nodes in the inet array to locals
1486 // int indInet = 0;
1487 // for ( unsigned int iEle = 0; iEle < isegn.size(); iEle++ ) {
1488 // int nne = nnet[ iEle ];
1489 // for ( int ien = 0; ien < nne; ien++ ) {
1490 //
1491 // int indGlob = inet[indInet];
1492 // // map it to local node
1493 // Global2LocalMap_::iterator pos = global2LocalNodeMap.find( indGlob );
1494 // ASSERT_PERMANENT( pos != global2LocalNodeMap.end())(indGlob).error("Cannot remap node index to local indices. \n " );
1495 // int indLoc = static_cast<int> ( pos -> second );
1496 //
1497 // // store the node
1498 // inet[ indInet++ ] = indLoc;
1499 // }
1500 // }
1501 //
1502 // int numNodes = size;
1503 // int numDofsInt = size;
1504 // int spaceDim = 3; // TODO: what is the proper value here?
1505 // int meshDim = elDimMax;
1506 //
1507 // /**
1508 // * We need:
1509 // * - local to global element map (possibly mesh->el_4_loc
1510 // * - inet, nnet - local dof numbers per element, local numbering of only those dofs that are on owned elements
1511 // * 1. collect DH local dof indices on elements, manage map from DH local indices to BDDC local dof indices
1512 // * 2. map collected DH indices to BDDC indices using the map
1513 // * - local BDDC dofs to global dofs, use DH to BDDC map with DH local to global map
1514 // * - XYZ - permuted, collect in main loop into array of size of all DH local dofs, compress and rearrange latter
1515 // * - element_permeability - in main loop
1516 // */
1517 // bddc_ls -> load_mesh( LinSys_BDDC::BDDCMatrixType::SPD_VIA_SYMMETRICGENERAL, spaceDim, numNodes, numDofsInt, inet, nnet, nndf, isegn, isngn, isngn, xyz, element_permeability, meshDim );
1518 // }
1519 
1520 
1521 
1522 
1523 //=============================================================================
1524 // DESTROY WATER MH SYSTEM STRUCTURE
1525 //=============================================================================
1527  if (output_object) delete output_object;
1528 
1529  if(time_ != nullptr)
1530  delete time_;
1531 
1532  if (read_init_cond_assembly_!=nullptr) {
1533  delete read_init_cond_assembly_;
1534  read_init_cond_assembly_ = nullptr;
1535  }
1536  if (mh_matrix_assembly_!=nullptr) {
1537  delete mh_matrix_assembly_;
1538  mh_matrix_assembly_ = nullptr;
1539  }
1540  if (reconstruct_schur_assembly_!=nullptr) {
1542  reconstruct_schur_assembly_ = nullptr;
1543  }
1544 }
1545 
1546 
1547 /// Helper method fills range (min and max) of given component
1548 void dofs_range(unsigned int n_dofs, unsigned int &min, unsigned int &max, unsigned int component) {
1549  if (component==0) {
1550  min = 0;
1551  max = n_dofs/2;
1552  } else if (component==1) {
1553  min = n_dofs/2;
1554  max = (n_dofs+1)/2;
1555  } else {
1556  min = (n_dofs+1)/2;
1557  max = n_dofs;
1558  }
1559 }
1560 
1561 
1563  ASSERT_LT(component, 3).error("Invalid component!");
1564  unsigned int i, n_dofs, min, max;
1565  std::vector<int> dof_vec;
1566  std::vector<LongIdx> dof_indices(eq_data_->dh_->max_elem_dofs());
1567  for ( DHCellAccessor dh_cell : eq_data_->dh_->own_range() ) {
1568  n_dofs = dh_cell.get_dof_indices(dof_indices);
1569  dofs_range(n_dofs, min, max, component);
1570  for (i=min; i<max; ++i) dof_vec.push_back(dof_indices[i]);
1571  }
1572  return dof_vec;
1573 }
1574 
1575 
1580 }
1581 
1582 
1584  this->read_init_cond_assembly_->assemble(eq_data_->dh_cr_);
1585 }
1586 
1587 
1588 //-----------------------------------------------------------------------------
1589 // vim: set cindent:
Functors of FieldModels used in Darcy flow module.
#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
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.
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.
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)
void reset()
Reset data members.
void init()
Initialize vectors, ...
Field< 3, FieldValue< 3 >::Scalar > water_source_density
Field< 3, FieldValue< 3 >::Scalar > ref_divergence
Field< 3, FieldValue< 3 >::Scalar > extra_storativity
Field< 3, FieldValue< 3 >::Scalar > field_ele_pressure
Externally added water source.
BCField< 3, FieldValue< 3 >::Enum > bc_type
Field< 3, FieldValue< 3 >::VectorFixed > gravity_field
Field< 3, FieldValue< 3 >::Scalar > sigma
Field< 3, FieldValue< 3 >::Scalar > init_pressure
Field< 3, FieldValue< 3 >::Scalar > storativity
BCField< 3, FieldValue< 3 >::Scalar > bc_pressure
BCField< 3, FieldValue< 3 >::Scalar > bc_flux
Field< 3, FieldValue< 3 >::TensorFixed > anisotropy
Field< 3, FieldValue< 3 >::Scalar > cross_section
Field< 3, FieldValue< 3 >::Scalar > conductivity
Field< 3, FieldValue< 3 >::VectorFixed > field_ele_velocity
Field< 3, FieldValue< 3 >::Scalar > field_edge_pressure
BCField< 3, FieldValue< 3 >::Scalar > bc_switch_piezo_head
Field< 3, FieldValue< 3 >::Scalar > field_ele_piezo_head
Field< 3, FieldValue< 3 >::Scalar > init_piezo_head
Same as previous but used in boundary fields.
BCField< 3, FieldValue< 3 >::Scalar > bc_robin_sigma
static const Input::Type::Selection & get_bc_type_selection()
Return a Selection corresponding to enum BC_Type.
Field< 3, FieldValue< 3 >::Scalar > extra_source
Externally added storativity.
Field< 3, FieldValue< 3 >::VectorFixed > ref_velocity
Precompute l2 difference outputs.
Field< 3, FieldValue< 3 >::VectorFixed > flux
EqFields()
Creation of all fields.
BCField< 3, FieldValue< 3 >::Scalar > bc_switch_pressure
BCField< 3, FieldValue< 3 >::VectorFixed > bc_gravity
Holds gravity vector acceptable in FieldModel.
BCField< 3, FieldValue< 3 >::Scalar > bc_piezo_head
Field< 3, FieldValue< 3 >::Scalar > ref_pressure
GenericAssemblyBase * mh_matrix_assembly_
void initialize() override
void print_matlab_matrix(string matlab_file)
Print darcy flow matrix in matlab format into a file.
virtual double solved_time() override
void zero_time_step() override
virtual void postprocess()
EqFields & eq_fields()
void solve_nonlinear()
Solve method common to zero_time_step and update solution.
static std::string equation_name()
virtual void initialize_specific()
static const int registrar
Registrar of class to factory.
DarcyFlowMHOutput * output_object
std::shared_ptr< EqData > eq_data_
unsigned int max_n_it_
GenericAssembly< ReadInitCondAssemblyLMH > * read_init_cond_assembly_
general assembly objects, hold assembly objects of appropriate dimension
bool data_changed_
GenericAssemblyBase * reconstruct_schur_assembly_
double tolerance_
virtual void initialize_asm()
Create and initialize assembly objects.
void update_solution() override
void allocate_mh_matrix()
unsigned int min_n_it_
void solve_time_step(bool output=true)
Solve the problem without moving to next time and without output.
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)
static const Input::Type::Record & type_field_descriptor()
void create_linear_system(Input::AbstractRecord rec)
static const Input::Type::Record & get_input_type()
DarcyLMH(Mesh &mesh, const Input::Record in_rec, TimeGovernor *tm=nullptr)
CREATE AND FILL GLOBAL MH MATRIX OF THE WATER MODEL.
friend class DarcyFlowMHOutput
virtual bool zero_time_term(bool time_global=false)
std::shared_ptr< EqFields > eq_fields_
static const Input::Type::Selection & get_mh_mortar_selection()
Selection for enum MortarMethod.
std::shared_ptr< Balance > balance_
void init_eq_data()
virtual void accept_time_step()
postprocess velocity field (add sources)
virtual ~DarcyLMH() override
virtual void output_data() override
Write computed fields.
LinSys & lin_sys_schur()
Getter for the linear system of the 2. Schur complement.
virtual void assembly_linear_system()
virtual void read_init_cond_asm()
Call assemble of read_init_cond_assembly_.
unsigned int n_neighs_vb() const
Return number of neighbours.
Definition: elements.h:65
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
static Input::Type::Record & user_fields_template(std::string equation_name)
Template Record with common key user_fields for derived equations.
Definition: equation.cc:46
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
FieldCommon & input_selection(Input::Type::Selection element_selection)
FieldCommon & description(const string &description)
FieldCommon & flags(FieldFlag::Flags::Mask mask)
FieldCommon & name(const string &name)
FieldCommon & set_limits(double min, double max=std::numeric_limits< double >::max())
FieldCommon & units(const UnitSI &units)
Set basic units of the field.
FieldCommon & input_default(const string &input_default)
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
void set_default_fieldset()
Definition: field_set.hh:408
auto disable_where(const Field< spacedim, typename FieldValue< spacedim >::Enum > &control_field, const vector< FieldEnum > &value_list) -> Field &
Definition: field.impl.hh:193
Dedicated class for storing path to input and output files.
Definition: file_path.hh:54
@ output_file
Definition: file_path.hh:69
virtual void assemble(std::shared_ptr< DOFHandlerMultiDim > dh)=0
void assemble(std::shared_ptr< DOFHandlerMultiDim > dh) override
General assemble methods.
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.
static const Input::Type::Record & get_input_type()
Definition: linsys_PETSC.cc:32
void set_solution(Vec sol_vec)
Definition: linsys.hh:290
void set_matrix_changed()
Definition: linsys.hh:212
virtual void set_from_input(const Input::Record in_rec)
Definition: linsys.hh:641
virtual void start_add_assembly()
Definition: linsys.hh:341
virtual void finish_assembly()=0
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
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
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.
void view(const char *name="") 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.
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.
Lumped mixed-hybrid model of linear Darcy flow, possibly unsteady.
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
Classes with algorithms for computation of intersections of meshes.
Wrappers for linear systems based on MPIAIJ and MATIS format.
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
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:933
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...
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.
Basic time management class.