Flow123d  master-f44eb46
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("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 
662  /* TODO:
663  * - Allow solution reconstruction (pressure and velocity) from initial condition on user request.
664  * - Steady solution as an intitial condition may be forced by setting inti_time =-1, and set data for the steady solver in that time.
665  * Solver should be able to switch from and to steady case depending on the zero time term.
666  */
667 
669 
670  // zero_time_term means steady case
671  eq_data_->use_steady_assembly_ = zero_time_term();
672 
673  eq_data_->p_edge_solution.zero_entries();
674 
675  if (eq_data_->use_steady_assembly_) { // steady case
676  MessageOut() << "Flow zero time step - steady case\n";
677  //read_initial_condition(); // Possible solution guess for steady case.
678  solve_nonlinear(); // with right limit data
679  } else {
680  MessageOut() << "Flow zero time step - unsteady case\n";
681  eq_data_->time_step_ = time_->dt();
683  this->read_init_cond_asm();
684  accept_time_step(); // accept zero time step, i.e. initial condition
685 
686 
687  // we reconstruct the initial solution here
688  // during the reconstruction assembly:
689  // - the balance objects are actually allocated
690  // - the full solution vector is computed
691  START_TIMER("DarcyFlowMH::reconstruct_solution_from_schur");
693  END_TIMER("DarcyFlowMH::reconstruct_solution_from_schur");
694  }
695  //solution_output(T,right_limit); // data for time T in any case
696  output_data();
697 }
698 
699 //=============================================================================
700 // COMPOSE and SOLVE WATER MH System possibly through Schur complements
701 //=============================================================================
703 {
704  START_TIMER("Solving MH system");
705 
706  time_->next_time();
707 
708  time_->view("DARCY"); //time governor information output
709 
710  solve_time_step();
711 
712  eq_data_->full_solution.local_to_ghost_begin();
713  eq_data_->full_solution.local_to_ghost_end();
714 }
715 
716 void DarcyLMH::solve_time_step(bool output)
717 {
719  bool zero_time_term_from_left=zero_time_term();
720 
721  bool jump_time = eq_fields_->storativity.is_jump_time();
722  if (! zero_time_term_from_left) {
723  MessageOut() << "Flow time step - unsteady case\n";
724  // time term not treated as zero
725  // Unsteady solution up to the T.
726 
727  // this flag is necesssary for switching BC to avoid setting zero neumann on the whole boundary in the steady case
728  eq_data_->use_steady_assembly_ = false;
729 
730  solve_nonlinear(); // with left limit data
731  if(output)
733  if (jump_time) {
734  WarningOut() << "Output of solution discontinuous in time not supported yet.\n";
735  //solution_output(T, left_limit); // output use time T- delta*dt
736  //output_data();
737  }
738  }
739 
740  if (time_->is_end()) {
741  // output for unsteady case, end_time should not be the jump time
742  // but rether check that
743  if (! zero_time_term_from_left && ! jump_time && output)
744  output_data();
745  return;
746  }
747 
749  bool zero_time_term_from_right=zero_time_term();
750  if (zero_time_term_from_right) {
751  MessageOut() << "Flow time step - steady case\n";
752  // this flag is necesssary for switching BC to avoid setting zero neumann on the whole boundary in the steady case
753  eq_data_->use_steady_assembly_ = true;
754  solve_nonlinear(); // with right limit data
755  if(output)
757 
758  } else if (! zero_time_term_from_left && jump_time) {
759  WarningOut() << "Discontinuous time term not supported yet.\n";
760  //solution_transfer(); // internally call set_time(T, left) and set_time(T,right) again
761  //solve_nonlinear(); // with right limit data
762  }
763  //solution_output(T,right_limit); // data for time T in any case
764  if (output)
765  output_data();
766 }
767 
768 bool DarcyLMH::zero_time_term(bool time_global) {
769  if (time_global) {
770  return (eq_fields_->storativity.input_list_size() == 0);
771  } else {
772  return eq_fields_->storativity.field_result(mesh_->region_db().get_region_set("BULK")) == result_zeros;
773  }
774 }
775 
776 
778 {
779 
781  double residual_norm = lin_sys_schur().compute_residual();
782  eq_data_->nonlinear_iteration_ = 0;
783  MessageOut().fmt("[nonlinear solver] norm of initial residual: {}\n", residual_norm);
784 
785  // Reduce is_linear flag.
786  int is_linear_common;
787  MPI_Allreduce(&(eq_data_->is_linear), &is_linear_common,1, MPI_INT ,MPI_MIN,PETSC_COMM_WORLD);
788 
789  Input::Record nl_solver_rec = input_record_.val<Input::Record>("nonlinear_solver");
790  this->tolerance_ = nl_solver_rec.val<double>("tolerance");
791  this->max_n_it_ = nl_solver_rec.val<unsigned int>("max_it");
792  this->min_n_it_ = nl_solver_rec.val<unsigned int>("min_it");
793  if (this->min_n_it_ > this->max_n_it_) this->min_n_it_ = this->max_n_it_;
794 
795  if (! is_linear_common) {
796  // set tolerances of the linear solver unless they are set by user.
797  lin_sys_schur().set_tolerances(0.1*this->tolerance_, 0.01*this->tolerance_, 10000, 100);
798  }
799  vector<double> convergence_history;
800 
801  while (eq_data_->nonlinear_iteration_ < this->min_n_it_ ||
802  (residual_norm > this->tolerance_ && eq_data_->nonlinear_iteration_ < this->max_n_it_ )) {
803  ASSERT_EQ( convergence_history.size(), eq_data_->nonlinear_iteration_ );
804  convergence_history.push_back(residual_norm);
805 
806  // print_matlab_matrix("matrix_" + std::to_string(time_->step().index()) + "_it_" + std::to_string(nonlinear_iteration_));
807  // stagnation test
808  if (convergence_history.size() >= 5 &&
809  convergence_history[ convergence_history.size() - 1]/convergence_history[ convergence_history.size() - 2] > 0.9 &&
810  convergence_history[ convergence_history.size() - 1]/convergence_history[ convergence_history.size() - 5] > 0.8) {
811  // stagnation
812  if (input_record_.val<Input::Record>("nonlinear_solver").val<bool>("converge_on_stagnation")) {
813  WarningOut().fmt("Accept solution on stagnation. Its: {} Residual: {}\n", eq_data_->nonlinear_iteration_, residual_norm);
814  break;
815  } else {
816  THROW(ExcSolverDiverge() << EI_Reason("Stagnation."));
817  }
818  }
819 
820  if (! is_linear_common){
821  eq_data_->p_edge_solution_previous.copy_from(eq_data_->p_edge_solution);
822  eq_data_->p_edge_solution_previous.local_to_ghost_begin();
823  eq_data_->p_edge_solution_previous.local_to_ghost_end();
824  }
825 
827  MessageOut().fmt("[schur solver] lin. it: {}, reason: {}, residual: {}\n",
828  si.n_iterations, si.converged_reason, lin_sys_schur().compute_residual());
829 
830  eq_data_->nonlinear_iteration_++;
831 
832  // hack to make BDDC work with empty compute_residual
833  if (is_linear_common){
834  // we want to print this info in linear (and steady) case
835  residual_norm = lin_sys_schur().compute_residual();
836  MessageOut().fmt("[nonlinear solver] lin. it: {}, reason: {}, residual: {}\n",
837  si.n_iterations, si.converged_reason, residual_norm);
838  break;
839  }
840  data_changed_=true; // force reassembly for non-linear case
841 
842  double alpha = 1; // how much of new solution
843  VecAXPBY(eq_data_->p_edge_solution.petsc_vec(), (1-alpha), alpha, eq_data_->p_edge_solution_previous.petsc_vec());
844 
845  //LogOut().fmt("Linear solver ended with reason: {} \n", si.converged_reason );
846  //ASSERT_PERMANENT_GE( si.converged_reason, 0).error("Linear solver failed to converge.\n");
848 
849  residual_norm = lin_sys_schur().compute_residual();
850  MessageOut().fmt("[nonlinear solver] it: {} lin. it: {}, reason: {}, residual: {}\n",
851  eq_data_->nonlinear_iteration_, si.n_iterations, si.converged_reason, residual_norm);
852  }
853 
854 // reconstruct_solution_from_schur(eq_data_->multidim_assembler);
855  START_TIMER("DarcyFlowMH::reconstruct_solution_from_schur");
857  END_TIMER("DarcyFlowMH::reconstruct_solution_from_schur");
858 
859  // adapt timestep
860  if (! this->zero_time_term()) {
861  double mult = 1.0;
862  if (eq_data_->nonlinear_iteration_ < 3) mult = 1.6;
863  if (eq_data_->nonlinear_iteration_ > 7) mult = 0.7;
864  time_->set_upper_constraint(time_->dt() * mult, "Darcy adaptivity.");
865  // int result = time_->set_upper_constraint(time_->dt() * mult, "Darcy adaptivity.");
866  //DebugOut().fmt("time adaptivity, res: {} it: {} m: {} dt: {} edt: {}\n", result, nonlinear_iteration_, mult, time_->dt(), time_->estimate_dt());
867  }
868 }
869 
870 
872 {
873  eq_data_->p_edge_solution_previous_time.copy_from(eq_data_->p_edge_solution);
874  eq_data_->p_edge_solution_previous_time.local_to_ghost_begin();
875  eq_data_->p_edge_solution_previous_time.local_to_ghost_end();
876 }
877 
878 
880  START_TIMER("Darcy output data");
881 
882  // print_matlab_matrix("matrix_" + std::to_string(time_->step().index()));
883 
884  //time_->view("DARCY"); //time governor information output
885  this->output_object->output();
886 
887 
888  START_TIMER("Darcy balance output");
889  balance_->calculate_cumulative(eq_data_->water_balance_idx, eq_data_->full_solution.petsc_vec());
890  balance_->calculate_instant(eq_data_->water_balance_idx, eq_data_->full_solution.petsc_vec());
891  balance_->output();
892 }
893 
894 
895 //double DarcyLMH::solution_precision() const
896 //{
897 // return eq_data_->lin_sys_schur->get_solution_precision();
898 //}
899 
900 
901 // ===========================================================================================
902 //
903 // MATRIX ASSEMBLY - we use abstract assembly routine, where LS Mat/Vec SetValues
904 // are in fact pointers to allocating or filling functions - this is governed by Linsystem roitunes
905 //
906 // =======================================================================================
907 //void DarcyLMH::assembly_mh_matrix(FMT_UNUSED MultidimAssembly& assembler)
908 //{
909 // START_TIMER("DarcyLMH::assembly_steady_mh_matrix");
910 //
911 // // DebugOut() << "assembly_mh_matrix \n";
912 // // set auxiliary flag for switchting Dirichlet like BC
913 // eq_data_->force_no_neumann_bc = eq_data_->use_steady_assembly_ && (eq_data_->nonlinear_iteration_ == 0);
914 //
915 // balance_->start_flux_assembly(eq_data_->water_balance_idx);
916 // balance_->start_source_assembly(eq_data_->water_balance_idx);
917 // balance_->start_mass_assembly(eq_data_->water_balance_idx);
918 //
919 // // TODO: try to move this into balance, or have it in the generic assembler class, that should perform the cell loop
920 // // including various pre- and post-actions
921 //// for ( DHCellAccessor dh_cell : eq_data_->dh_->own_range() ) {
922 //// unsigned int dim = dh_cell.dim();
923 //// assembler[dim-1]->assemble(dh_cell);
924 //// }
925 // this->mh_matrix_assembly_->assemble(eq_data_->dh_);
926 //
927 //
928 // balance_->finish_mass_assembly(eq_data_->water_balance_idx);
929 // balance_->finish_source_assembly(eq_data_->water_balance_idx);
930 // balance_->finish_flux_assembly(eq_data_->water_balance_idx);
931 //
932 //}
933 
934 
936 {
937  START_TIMER("DarcyLMH::allocate_mh_matrix");
938 
939  // to make space for second schur complement, max. 10 neighbour edges of one el.
940  double zeros[100000];
941  for(int i=0; i<100000; i++) zeros[i] = 0.0;
942 
943  std::vector<LongIdx> tmp_rows;
944  tmp_rows.reserve(200);
945 
946  std::vector<LongIdx> dofs, dofs_ngh;
947  dofs.reserve(eq_data_->dh_cr_->max_elem_dofs());
948  dofs_ngh.reserve(eq_data_->dh_cr_->max_elem_dofs());
949 
950  // DebugOut() << "Allocate new schur\n";
951  for ( DHCellAccessor dh_cell : eq_data_->dh_cr_->own_range() ) {
952  ElementAccessor<3> ele = dh_cell.elm();
953 
954  const uint ndofs = dh_cell.n_dofs();
955  dofs.resize(dh_cell.n_dofs());
956  dh_cell.get_dof_indices(dofs);
957 
958  int* dofs_ptr = dofs.data();
959  lin_sys_schur().mat_set_values(ndofs, dofs_ptr, ndofs, dofs_ptr, zeros);
960 
961  tmp_rows.clear();
962 
963  // compatible neighborings rows
964  unsigned int n_neighs = ele->n_neighs_vb();
965  for ( DHCellSide neighb_side : dh_cell.neighb_sides() ) {
966  // every compatible connection adds a 2x2 matrix involving
967  // current element pressure and a connected edge pressure
968 
969  // read neighbor dofs (dh_cr dofhandler)
970  // neighbor cell owning neighb_side
971  DHCellAccessor dh_neighb_cell = neighb_side.cell();
972 
973  const uint ndofs_ngh = dh_neighb_cell.n_dofs();
974  dofs_ngh.resize(ndofs_ngh);
975  dh_neighb_cell.get_dof_indices(dofs_ngh);
976 
977  // local index of pedge dof on neighboring cell
978  tmp_rows.push_back(dofs_ngh[neighb_side.side().side_idx()]);
979  }
980 
981  lin_sys_schur().mat_set_values(ndofs, dofs_ptr, n_neighs, tmp_rows.data(), zeros); // (edges) x (neigh edges)
982  lin_sys_schur().mat_set_values(n_neighs, tmp_rows.data(), ndofs, dofs_ptr, zeros); // (neigh edges) x (edges)
983  lin_sys_schur().mat_set_values(n_neighs, tmp_rows.data(), n_neighs, tmp_rows.data(), zeros); // (neigh edges) x (neigh edges)
984 
985  tmp_rows.clear();
986 // if (eq_data_->mortar_method_ != NoMortar) {
987 // auto &isec_list = mesh_->mixed_intersections().element_intersections_[ele.idx()];
988 // for(auto &isec : isec_list ) {
989 // IntersectionLocalBase *local = isec.second;
990 // DHCellAccessor dh_cell_slave = eq_data_->dh_cr_->cell_accessor_from_element(local->bulk_ele_idx());
991 //
992 // const uint ndofs_slave = dh_cell_slave.n_dofs();
993 // dofs_ngh.resize(ndofs_slave);
994 // dh_cell_slave.get_dof_indices(dofs_ngh);
995 //
996 // //DebugOut().fmt("Alloc: {} {}", ele.idx(), local->bulk_ele_idx());
997 // for(unsigned int i_side=0; i_side < dh_cell_slave.elm()->n_sides(); i_side++) {
998 // tmp_rows.push_back( dofs_ngh[i_side] );
999 // //DebugOut() << "aedge" << print_var(tmp_rows[tmp_rows.size()-1]);
1000 // }
1001 // }
1002 // }
1003 
1004  lin_sys_schur().mat_set_values(ndofs, dofs_ptr, tmp_rows.size(), tmp_rows.data(), zeros); // master edges x slave edges
1005  lin_sys_schur().mat_set_values(tmp_rows.size(), tmp_rows.data(), ndofs, dofs_ptr, zeros); // slave edges x master edges
1006  lin_sys_schur().mat_set_values(tmp_rows.size(), tmp_rows.data(), tmp_rows.size(), tmp_rows.data(), zeros); // slave edges x slave edges
1007  }
1008  // DebugOut() << "end Allocate new schur\n";
1009 
1010  // int local_dofs[10];
1011  // unsigned int nsides;
1012  // for ( DHCellAccessor dh_cell : eq_data_->dh_->own_range() ) {
1013  // LocalElementAccessorBase<3> ele_ac(dh_cell);
1014  // nsides = ele_ac.n_sides();
1015 
1016  // //allocate at once matrix [sides,ele,edges]x[sides,ele,edges]
1017  // loc_size = 1 + 2*nsides;
1018  // unsigned int i_side = 0;
1019 
1020  // for (; i_side < nsides; i_side++) {
1021  // local_dofs[i_side] = ele_ac.side_row(i_side);
1022  // local_dofs[i_side+nsides] = ele_ac.edge_row(i_side);
1023  // }
1024  // local_dofs[i_side+nsides] = ele_ac.ele_row();
1025  // int * edge_rows = local_dofs + nsides;
1026  // //int ele_row = local_dofs[0];
1027 
1028  // // whole local MH matrix
1029  // ls->mat_set_values(loc_size, local_dofs, loc_size, local_dofs, zeros);
1030 
1031 
1032  // // compatible neighborings rows
1033  // unsigned int n_neighs = ele_ac.element_accessor()->n_neighs_vb();
1034  // unsigned int i=0;
1035  // for ( DHCellSide neighb_side : dh_cell.neighb_sides() ) {
1036  // //for (unsigned int i = 0; i < n_neighs; i++) {
1037  // // every compatible connection adds a 2x2 matrix involving
1038  // // current element pressure and a connected edge pressure
1039  // Neighbour *ngh = ele_ac.element_accessor()->neigh_vb[i];
1040  // DHCellAccessor cell_higher_dim = eq_data_->dh_->cell_accessor_from_element(neighb_side.elem_idx());
1041  // LocalElementAccessorBase<3> acc_higher_dim( cell_higher_dim );
1042  // for (unsigned int j = 0; j < neighb_side.element().dim()+1; j++)
1043  // if (neighb_side.element()->edge_idx(j) == ngh->edge_idx()) {
1044  // int neigh_edge_row = acc_higher_dim.edge_row(j);
1045  // tmp_rows.push_back(neigh_edge_row);
1046  // break;
1047  // }
1048  // //DebugOut() << "CC" << print_var(tmp_rows[i]);
1049  // ++i;
1050  // }
1051 
1052  // // allocate always also for schur 2
1053  // ls->mat_set_values(nsides+1, edge_rows, n_neighs, tmp_rows.data(), zeros); // (edges, ele) x (neigh edges)
1054  // ls->mat_set_values(n_neighs, tmp_rows.data(), nsides+1, edge_rows, zeros); // (neigh edges) x (edges, ele)
1055  // ls->mat_set_values(n_neighs, tmp_rows.data(), n_neighs, tmp_rows.data(), zeros); // (neigh edges) x (neigh edges)
1056 
1057  // tmp_rows.clear();
1058 
1059  // if (eq_data_->mortar_method_ != NoMortar) {
1060  // auto &isec_list = mesh_->mixed_intersections().element_intersections_[ele_ac.ele_global_idx()];
1061  // for(auto &isec : isec_list ) {
1062  // IntersectionLocalBase *local = isec.second;
1063  // LocalElementAccessorBase<3> slave_acc( eq_data_->dh_->cell_accessor_from_element(local->bulk_ele_idx()) );
1064  // //DebugOut().fmt("Alloc: {} {}", ele_ac.ele_global_idx(), local->bulk_ele_idx());
1065  // for(unsigned int i_side=0; i_side < slave_acc.dim()+1; i_side++) {
1066  // tmp_rows.push_back( slave_acc.edge_row(i_side) );
1067  // //DebugOut() << "aedge" << print_var(tmp_rows[tmp_rows.size()-1]);
1068  // }
1069  // }
1070  // }
1071  // /*
1072  // for(unsigned int i_side=0; i_side < ele_ac.element_accessor()->n_sides(); i_side++) {
1073  // DebugOut() << "aedge:" << print_var(edge_rows[i_side]);
1074  // }*/
1075 
1076  // ls->mat_set_values(nsides, edge_rows, tmp_rows.size(), tmp_rows.data(), zeros); // master edges x neigh edges
1077  // ls->mat_set_values(tmp_rows.size(), tmp_rows.data(), nsides, edge_rows, zeros); // neigh edges x master edges
1078  // ls->mat_set_values(tmp_rows.size(), tmp_rows.data(), tmp_rows.size(), tmp_rows.data(), zeros); // neigh edges x neigh edges
1079 
1080  // }
1081 /*
1082  // alloc edge diagonal entries
1083  if(rank == 0)
1084  for( vector<Edge>::iterator edg = mesh_->edges.begin(); edg != mesh_->edges.end(); ++edg) {
1085  int edg_idx = mh_dh.row_4_edge[edg->side(0)->edge_idx()];
1086 
1087 // for( vector<Edge>::iterator edg2 = mesh_->edges.begin(); edg2 != mesh_->edges.end(); ++edg2){
1088 // int edg_idx2 = mh_dh.row_4_edge[edg2->side(0)->edge_idx()];
1089 // if(edg_idx == edg_idx2){
1090 // DBGCOUT(<< "P[ " << rank << " ] " << "edg alloc: " << edg_idx << " " << edg_idx2 << "\n");
1091  ls->mat_set_value(edg_idx, edg_idx, 0.0);
1092 // }
1093 // }
1094  }
1095  */
1096  /*
1097  if (mortar_method_ == MortarP0) {
1098  P0_CouplingAssembler(*this).assembly(*ls);
1099  } else if (mortar_method_ == MortarP1) {
1100  P1_CouplingAssembler(*this).assembly(*ls);
1101  }*/
1102 }
1103 
1104 
1105 
1106 /*******************************************************************************
1107  * COMPOSE WATER MH MATRIX WITHOUT SCHUR COMPLEMENT
1108  ******************************************************************************/
1109 
1111 
1112  START_TIMER("preallocation");
1113 
1114  // if (schur0 == NULL) { // create Linear System for MH matrix
1115 
1116 // if (in_rec.type() == LinSys_BDDC::get_input_type()) {
1117 // #ifdef FLOW123D_HAVE_BDDCML
1118 // WarningOut() << "For BDDC no Schur complements are used.";
1119 // n_schur_compls = 0;
1120 // LinSys_BDDC *ls = new LinSys_BDDC(&(*eq_data_->dh_->distr()),
1121 // true); // swap signs of matrix and rhs to make the matrix SPD
1122 // ls->set_from_input(in_rec);
1123 // ls->set_solution( eq_data_->full_solution.petsc_vec() );
1124 // // possible initialization particular to BDDC
1125 // START_TIMER("BDDC set mesh data");
1126 // set_mesh_data_for_bddc(ls);
1127 // schur0=ls;
1128 // END_TIMER("BDDC set mesh data");
1129 // #else
1130 // Exception
1131 // THROW( ExcBddcmlNotSupported() );
1132 // #endif // FLOW123D_HAVE_BDDCML
1133 // }
1134 // else
1135  if (in_rec.type() == LinSys_PETSC::get_input_type()) {
1136  // use PETSC for serial case even when user wants BDDC
1137 
1138  eq_data_->lin_sys_schur = std::make_shared<LinSys_PETSC>( &(*eq_data_->dh_cr_->distr()) );
1139  lin_sys_schur().set_from_input(in_rec);
1141  lin_sys_schur().set_solution( eq_data_->p_edge_solution.petsc_vec() );
1143  ((LinSys_PETSC *)&lin_sys_schur())->set_initial_guess_nonzero(true);
1144 
1145 // LinSys_PETSC *schur1, *schur2;
1146 
1147 // if (n_schur_compls == 0) {
1148 // LinSys_PETSC *ls = new LinSys_PETSC( &(*eq_data_->dh_->distr()) );
1149 
1150 // // temporary solution; we have to set precision also for sequantial case of BDDC
1151 // // final solution should be probably call of direct solver for oneproc case
1152 // // if (in_rec.type() != LinSys_BDDC::get_input_type()) ls->set_from_input(in_rec);
1153 // // else {
1154 // // ls->LinSys::set_from_input(in_rec); // get only common options
1155 // // }
1156 // ls->set_from_input(in_rec);
1157 
1158 // // ls->set_solution( eq_data_->full_solution.petsc_vec() );
1159 // schur0=ls;
1160 // } else {
1161 // IS is;
1162 // auto side_dofs_vec = get_component_indices_vec(0);
1163 
1164 // ISCreateGeneral(PETSC_COMM_SELF, side_dofs_vec.size(), &(side_dofs_vec[0]), PETSC_COPY_VALUES, &is);
1165 // //ISView(is, PETSC_VIEWER_STDOUT_SELF);
1166 // //ASSERT_PERMANENT(err == 0).error("Error in ISCreateStride.");
1167 
1168 // SchurComplement *ls = new SchurComplement(&(*eq_data_->dh_->distr()), is);
1169 
1170 // // make schur1
1171 // Distribution *ds = ls->make_complement_distribution();
1172 // if (n_schur_compls==1) {
1173 // schur1 = new LinSys_PETSC(ds);
1174 // schur1->set_positive_definite();
1175 // } else {
1176 // IS is;
1177 // auto elem_dofs_vec = get_component_indices_vec(1);
1178 
1179 // const PetscInt *b_indices;
1180 // ISGetIndices(ls->IsB, &b_indices);
1181 // uint b_size = ls->loc_size_B;
1182 // for(uint i_b=0, i_bb=0; i_b < b_size && i_bb < elem_dofs_vec.size(); i_b++) {
1183 // if (b_indices[i_b] == elem_dofs_vec[i_bb])
1184 // elem_dofs_vec[i_bb++] = i_b + ds->begin();
1185 // }
1186 // ISRestoreIndices(ls->IsB, &b_indices);
1187 
1188 
1189 // ISCreateGeneral(PETSC_COMM_SELF, elem_dofs_vec.size(), &(elem_dofs_vec[0]), PETSC_COPY_VALUES, &is);
1190 // //ISView(is, PETSC_VIEWER_STDOUT_SELF);
1191 // //ASSERT_PERMANENT(err == 0).error("Error in ISCreateStride.");
1192 // SchurComplement *ls1 = new SchurComplement(ds, is); // is is deallocated by SchurComplement
1193 // ls1->set_negative_definite();
1194 
1195 // // make schur2
1196 // schur2 = new LinSys_PETSC( ls1->make_complement_distribution() );
1197 // schur2->set_positive_definite();
1198 // ls1->set_complement( schur2 );
1199 // schur1 = ls1;
1200 // }
1201 // ls->set_complement( schur1 );
1202 // ls->set_from_input(in_rec);
1203 // // ls->set_solution( eq_data_->full_solution.petsc_vec() );
1204 // schur0=ls;
1205  // }
1206 
1207  START_TIMER("PETSC PREALLOCATION");
1209 
1211 
1212  eq_data_->full_solution.zero_entries();
1213  eq_data_->p_edge_solution.zero_entries();
1214  END_TIMER("PETSC PREALLOCATION");
1215  }
1216  else {
1217  THROW( ExcUnknownSolver() );
1218  }
1219 
1220  END_TIMER("preallocation");
1221 }
1222 
1224 {}
1225 
1226 //void DarcyLMH::reconstruct_solution_from_schur(MultidimAssembly& assembler)
1227 //{
1228 // START_TIMER("DarcyFlowMH::reconstruct_solution_from_schur");
1229 //
1230 // eq_data_->full_solution.zero_entries();
1231 // eq_data_->p_edge_solution.local_to_ghost_begin();
1232 // eq_data_->p_edge_solution.local_to_ghost_end();
1233 //
1234 // balance_->start_flux_assembly(eq_data_->water_balance_idx);
1235 // balance_->start_source_assembly(eq_data_->water_balance_idx);
1236 // balance_->start_mass_assembly(eq_data_->water_balance_idx);
1237 //
1238 // for ( DHCellAccessor dh_cell : eq_data_->dh_->own_range() ) {
1239 // unsigned int dim = dh_cell.dim();
1240 // assembler[dim-1]->assemble_reconstruct(dh_cell);
1241 // }
1242 //
1243 // eq_data_->full_solution.local_to_ghost_begin();
1244 // eq_data_->full_solution.local_to_ghost_end();
1245 //
1246 // balance_->finish_mass_assembly(eq_data_->water_balance_idx);
1247 // balance_->finish_source_assembly(eq_data_->water_balance_idx);
1248 // balance_->finish_flux_assembly(eq_data_->water_balance_idx);
1249 //}
1250 
1252  START_TIMER("DarcyFlowMH::assembly_linear_system");
1253 // DebugOut() << "DarcyLMH::assembly_linear_system\n";
1254 
1255  eq_data_->p_edge_solution.local_to_ghost_begin();
1256  eq_data_->p_edge_solution.local_to_ghost_end();
1257 
1258  eq_data_->is_linear=true;
1259  //DebugOut() << "Assembly linear system\n";
1260 // if (data_changed_) {
1261 // data_changed_ = false;
1262  {
1263  //DebugOut() << "Data changed\n";
1264  // currently we have no optimization for cases when just time term data or RHS data are changed
1265  START_TIMER("full assembly");
1266 // if (typeid(*schur0) != typeid(LinSys_BDDC)) {
1267 // schur0->start_add_assembly(); // finish allocation and create matrix
1268 // schur_compl->start_add_assembly();
1269 // }
1270 
1272 
1275 
1276  eq_data_->time_step_ = time_->dt();
1277 
1278  START_TIMER("DarcyLMH::assembly_steady_mh_matrix");
1279  this->mh_matrix_assembly_->assemble(eq_data_->dh_);; // fill matrix
1280  END_TIMER("DarcyLMH::assembly_steady_mh_matrix");
1281 // assembly_mh_matrix( eq_data_->multidim_assembler ); // fill matrix
1282 
1285 
1286  // print_matlab_matrix("matrix");
1287  }
1288 }
1289 
1290 
1291 void DarcyLMH::print_matlab_matrix(std::string matlab_file)
1292 {
1293  std::string output_file;
1294 
1295  // compute h_min for different dimensions
1296  double d_max = std::numeric_limits<double>::max();
1297  double h1 = d_max, h2 = d_max, h3 = d_max;
1298  double he2 = d_max, he3 = d_max;
1299  for (auto ele : mesh_->elements_range()) {
1300  switch(ele->dim()){
1301  case 1: h1 = std::min(h1,ele.measure()); break;
1302  case 2: h2 = std::min(h2,ele.measure()); break;
1303  case 3: h3 = std::min(h3,ele.measure()); break;
1304  }
1305 
1306  for (unsigned int j=0; j<ele->n_sides(); j++) {
1307  switch(ele->dim()){
1308  case 2: he2 = std::min(he2, ele.side(j)->measure()); break;
1309  case 3: he3 = std::min(he3, ele.side(j)->measure()); break;
1310  }
1311  }
1312  }
1313  if(h1 == d_max) h1 = 0;
1314  if(h2 == d_max) h2 = 0;
1315  if(h3 == d_max) h3 = 0;
1316  if(he2 == d_max) he2 = 0;
1317  if(he3 == d_max) he3 = 0;
1318 
1319  FILE * file;
1320  file = fopen(output_file.c_str(),"a");
1321  fprintf(file, "nA = %d;\n", eq_data_->dh_cr_disc_->distr()->size());
1322  fprintf(file, "nB = %d;\n", eq_data_->dh_->mesh()->get_el_ds()->size());
1323  fprintf(file, "nBF = %d;\n", eq_data_->dh_cr_->distr()->size());
1324  fprintf(file, "h1 = %e;\nh2 = %e;\nh3 = %e;\n", h1, h2, h3);
1325  fprintf(file, "he2 = %e;\nhe3 = %e;\n", he2, he3);
1326  fclose(file);
1327 
1328  {
1329  output_file = FilePath(matlab_file + "_sch_new.m", FilePath::output_file);
1330  PetscViewer viewer;
1331  PetscViewerASCIIOpen(PETSC_COMM_WORLD, output_file.c_str(), &viewer);
1332  PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
1333  MatView( *const_cast<Mat*>(lin_sys_schur().get_matrix()), viewer);
1334  VecView( *const_cast<Vec*>(lin_sys_schur().get_rhs()), viewer);
1335  VecView( *const_cast<Vec*>(&(lin_sys_schur().get_solution())), viewer);
1336  VecView( *const_cast<Vec*>(&(eq_data_->full_solution.petsc_vec())), viewer);
1337  }
1338 }
1339 
1340 
1341 //template <int dim>
1342 //std::vector<arma::vec3> dof_points(DHCellAccessor cell, const Mapping<dim, 3> &mapping) {
1343 //
1344 //
1345 // vector<arma::vec::fixed<dim+1>> bary_dof_points = cell->fe()->dof_points();
1346 //
1347 // std::vector<arma::vec3> points(20);
1348 // points.resize(0);
1349 //
1350 //}
1351 //
1352 
1353 // void DarcyLMH::set_mesh_data_for_bddc(LinSys_BDDC * bddc_ls) {
1354 // START_TIMER("DarcyFlowMH_Steady::set_mesh_data_for_bddc");
1355 // // prepare mesh for BDDCML
1356 // // initialize arrays
1357 // // auxiliary map for creating coordinates of local dofs and global-to-local numbering
1358 // std::map<int, arma::vec3> localDofMap;
1359 // // connectivity for the subdomain, i.e. global dof numbers on element, stored element-by-element
1360 // // Indices of Nodes on Elements
1361 // std::vector<int> inet;
1362 // // number of degrees of freedom on elements - determines elementwise chunks of INET array
1363 // // Number of Nodes on Elements
1364 // std::vector<int> nnet;
1365 // // Indices of Subdomain Elements in Global Numbering - for local elements, their global indices
1366 // std::vector<int> isegn;
1367 //
1368 // // This array is currently not used in BDDCML, it was used as an interface scaling alternative to scaling
1369 // // by diagonal. It corresponds to the rho-scaling.
1370 // std::vector<double> element_permeability;
1371 //
1372 // // maximal and minimal dimension of elements
1373 // uint elDimMax = 1;
1374 // uint elDimMin = 3;
1375 // std::vector<LongIdx> cell_dofs_global(10);
1376 //
1377 //
1378 //
1379 // for ( DHCellAccessor dh_cell : eq_data_->dh_->own_range() ) {
1380 // // LocalElementAccessorBase<3> ele_ac(dh_cell);
1381 // // for each element, create local numbering of dofs as fluxes (sides), pressure (element centre), Lagrange multipliers (edges), compatible connections
1382 //
1383 // dh_cell.get_dof_indices(cell_dofs_global);
1384 //
1385 // inet.insert(inet.end(), cell_dofs_global.begin(), cell_dofs_global.end());
1386 // uint n_inet = cell_dofs_global.size();
1387 //
1388 //
1389 // uint dim = dh_cell.elm().dim();
1390 // elDimMax = std::max( elDimMax, dim );
1391 // elDimMin = std::min( elDimMin, dim );
1392 //
1393 // // TODO: this is consistent with previous implementation, but may be wrong as it use global element numbering
1394 // // used in sequential mesh, do global numbering of distributed elements.
1395 // isegn.push_back( dh_cell.elm_idx());
1396 //
1397 // // TODO: use FiniteElement::dof_points
1398 // for (unsigned int si=0; si<dh_cell.elm()->n_sides(); si++) {
1399 // arma::vec3 coord = dh_cell.elm().side(si)->centre();
1400 // // flux dof points
1401 // localDofMap.insert( std::make_pair( cell_dofs_global[si], coord ) );
1402 // // pressure trace dof points
1403 // localDofMap.insert( std::make_pair( cell_dofs_global[si+dim+2], coord ) );
1404 // }
1405 // // pressure dof points
1406 // arma::vec3 elm_centre = dh_cell.elm().centre();
1407 // localDofMap.insert( std::make_pair( cell_dofs_global[dim+1], elm_centre ) );
1408 //
1409 // // insert dofs related to compatible connections
1410 // //const Element *ele = dh_cell.elm().element();
1411 // for(DHCellSide side : dh_cell.neighb_sides()) {
1412 // uint neigh_dim = side.cell().elm().dim();
1413 // side.cell().get_dof_indices(cell_dofs_global);
1414 // int edge_row = cell_dofs_global[neigh_dim+2+side.side_idx()];
1415 // localDofMap.insert( std::make_pair( edge_row, side.centre() ) );
1416 // inet.push_back( edge_row );
1417 // n_inet++;
1418 // }
1419 // nnet.push_back(n_inet);
1420 //
1421 //
1422 // // version for rho scaling
1423 // // trace computation
1424 // double conduct = eq_fields_->conductivity.value( elm_centre , dh_cell.elm() );
1425 // auto aniso = eq_fields_->anisotropy.value( elm_centre , dh_cell.elm() );
1426 //
1427 // // compute mean on the diagonal
1428 // double coef = 0.;
1429 // for ( int i = 0; i < 3; i++) {
1430 // coef = coef + aniso.at(i,i);
1431 // }
1432 // // Maybe divide by cs
1433 // coef = conduct*coef / 3;
1434 //
1435 // ASSERT_PERMANENT_GT(coef, 0).error("Zero coefficient of hydrodynamic resistance.\n");
1436 // element_permeability.push_back( 1. / coef );
1437 // }
1438 // // uint i_inet = 0;
1439 // // for(int n_dofs : nnet) {
1440 // // DebugOut() << "nnet: " << n_dofs;
1441 // // for(int j=0; j < n_dofs; j++, i_inet++) {
1442 // // DebugOut() << "inet: " << inet[i_inet];
1443 // // }
1444 // // }
1445 //
1446 // auto distr = eq_data_->dh_->distr();
1447 // // for(auto pair : localDofMap) {
1448 // // DebugOut().every_proc() << "r: " << distr->myp() << " gi: " << pair.first << "xyz: " << pair.second[0];
1449 // //
1450 // // }
1451 //
1452 //
1453 // //convert set of dofs to vectors
1454 // // number of nodes (= dofs) on the subdomain
1455 // int numNodeSub = localDofMap.size();
1456 // //ASSERT_PERMANENT_EQ( (unsigned int)numNodeSub, eq_data_->dh_->lsize() );
1457 // // Indices of Subdomain Nodes in Global Numbering - for local nodes, their global indices
1458 // std::vector<int> isngn( numNodeSub );
1459 // // pseudo-coordinates of local nodes (i.e. dofs)
1460 // // they need not be exact, they are used just for some geometrical considerations in BDDCML,
1461 // // such as selection of corners maximizing area of a triangle, bounding boxes fro subdomains to
1462 // // find candidate neighbours etc.
1463 // std::vector<double> xyz( numNodeSub * 3 ) ;
1464 // int ind = 0;
1465 // std::map<int,arma::vec3>::iterator itB = localDofMap.begin();
1466 // for ( ; itB != localDofMap.end(); ++itB ) {
1467 // isngn[ind] = itB -> first;
1468 //
1469 // arma::vec3 coord = itB -> second;
1470 // for ( int j = 0; j < 3; j++ ) {
1471 // xyz[ j*numNodeSub + ind ] = coord[j];
1472 // }
1473 //
1474 // ind++;
1475 // }
1476 // localDofMap.clear();
1477 //
1478 // // Number of Nodal Degrees of Freedom
1479 // // nndf is trivially one - dofs coincide with nodes
1480 // std::vector<int> nndf( numNodeSub, 1 );
1481 //
1482 // // prepare auxiliary map for renumbering nodes
1483 // typedef std::map<int,int> Global2LocalMap_; //! type for storage of global to local map
1484 // Global2LocalMap_ global2LocalNodeMap;
1485 // for ( unsigned ind = 0; ind < isngn.size(); ++ind ) {
1486 // global2LocalNodeMap.insert( std::make_pair( static_cast<unsigned>( isngn[ind] ), ind ) );
1487 // }
1488 //
1489 // // renumber nodes in the inet array to locals
1490 // int indInet = 0;
1491 // for ( unsigned int iEle = 0; iEle < isegn.size(); iEle++ ) {
1492 // int nne = nnet[ iEle ];
1493 // for ( int ien = 0; ien < nne; ien++ ) {
1494 //
1495 // int indGlob = inet[indInet];
1496 // // map it to local node
1497 // Global2LocalMap_::iterator pos = global2LocalNodeMap.find( indGlob );
1498 // ASSERT_PERMANENT( pos != global2LocalNodeMap.end())(indGlob).error("Cannot remap node index to local indices. \n " );
1499 // int indLoc = static_cast<int> ( pos -> second );
1500 //
1501 // // store the node
1502 // inet[ indInet++ ] = indLoc;
1503 // }
1504 // }
1505 //
1506 // int numNodes = size;
1507 // int numDofsInt = size;
1508 // int spaceDim = 3; // TODO: what is the proper value here?
1509 // int meshDim = elDimMax;
1510 //
1511 // /**
1512 // * We need:
1513 // * - local to global element map (possibly mesh->el_4_loc
1514 // * - inet, nnet - local dof numbers per element, local numbering of only those dofs that are on owned elements
1515 // * 1. collect DH local dof indices on elements, manage map from DH local indices to BDDC local dof indices
1516 // * 2. map collected DH indices to BDDC indices using the map
1517 // * - local BDDC dofs to global dofs, use DH to BDDC map with DH local to global map
1518 // * - XYZ - permuted, collect in main loop into array of size of all DH local dofs, compress and rearrange latter
1519 // * - element_permeability - in main loop
1520 // */
1521 // bddc_ls -> load_mesh( LinSys_BDDC::BDDCMatrixType::SPD_VIA_SYMMETRICGENERAL, spaceDim, numNodes, numDofsInt, inet, nnet, nndf, isegn, isngn, isngn, xyz, element_permeability, meshDim );
1522 // }
1523 
1524 
1525 
1526 
1527 //=============================================================================
1528 // DESTROY WATER MH SYSTEM STRUCTURE
1529 //=============================================================================
1531  if (output_object) delete output_object;
1532 
1533  if(time_ != nullptr)
1534  delete time_;
1535 
1536  if (read_init_cond_assembly_!=nullptr) {
1537  delete read_init_cond_assembly_;
1538  read_init_cond_assembly_ = nullptr;
1539  }
1540  if (mh_matrix_assembly_!=nullptr) {
1541  delete mh_matrix_assembly_;
1542  mh_matrix_assembly_ = nullptr;
1543  }
1544  if (reconstruct_schur_assembly_!=nullptr) {
1546  reconstruct_schur_assembly_ = nullptr;
1547  }
1548 }
1549 
1550 
1551 /// Helper method fills range (min and max) of given component
1552 void dofs_range(unsigned int n_dofs, unsigned int &min, unsigned int &max, unsigned int component) {
1553  if (component==0) {
1554  min = 0;
1555  max = n_dofs/2;
1556  } else if (component==1) {
1557  min = n_dofs/2;
1558  max = (n_dofs+1)/2;
1559  } else {
1560  min = (n_dofs+1)/2;
1561  max = n_dofs;
1562  }
1563 }
1564 
1565 
1567  ASSERT_LT(component, 3).error("Invalid component!");
1568  unsigned int i, n_dofs, min, max;
1569  std::vector<int> dof_vec;
1570  std::vector<LongIdx> dof_indices(eq_data_->dh_->max_elem_dofs());
1571  for ( DHCellAccessor dh_cell : eq_data_->dh_->own_range() ) {
1572  n_dofs = dh_cell.get_dof_indices(dof_indices);
1573  dofs_range(n_dofs, min, max, component);
1574  for (i=min; i<max; ++i) dof_vec.push_back(dof_indices[i]);
1575  }
1576  return dof_vec;
1577 }
1578 
1579 
1584 }
1585 
1586 
1588  this->read_init_cond_assembly_->assemble(eq_data_->dh_cr_);
1589 }
1590 
1591 
1592 //-----------------------------------------------------------------------------
1593 // 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: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...
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.