28 #include "petscviewer.h"
29 #include "petscerror.h"
51 #include "flow/assembly_mh_old.hh"
52 #include "flow/darcy_flow_mh.hh"
85 .
add_value(
MortarP1,
"P1",
"Mortar space: P1 on intersections, using non-conforming pressures.")
93 "Homogeneous Neumann boundary condition\n(zero normal flux over the boundary).")
95 "Dirichlet boundary condition. "
96 "Specify the pressure head through the ``bc_pressure`` field "
97 "or the piezometric head through the ``bc_piezo_head`` field.")
98 .
add_value(
total_flux,
"total_flux",
"Flux boundary condition (combines Neumann and Robin type). "
99 "Water inflow equal to (($ \\delta_d(q_d^N + \\sigma_d (h_d^R - h_d) )$)). "
100 "Specify the water inflow by the ``bc_flux`` field, the transition coefficient by ``bc_robin_sigma`` "
101 "and the reference pressure head or piezometric head through ``bc_pressure`` or ``bc_piezo_head`` respectively.")
103 "Seepage face boundary condition. Pressure and inflow bounded from above. Boundary with potential seepage flow "
104 "is described by the pair of inequalities: "
105 "(($h_d \\le h_d^D$)) and (($ -\\boldsymbol q_d\\cdot\\boldsymbol n \\le \\delta q_d^N$)), where the equality holds in at least one of them. "
106 "Caution: setting (($q_d^N$)) strictly negative "
107 "may lead to an ill posed problem since a positive outflow is enforced. "
108 "Parameters (($h_d^D$)) and (($q_d^N$)) are given by the fields ``bc_switch_pressure`` (or ``bc_switch_piezo_head``) and ``bc_flux`` respectively."
111 "River boundary condition. For the water level above the bedrock, (($H_d > H_d^S$)), the Robin boundary condition is used with the inflow given by: "
112 "(( $ \\delta_d(q_d^N + \\sigma_d(H_d^D - H_d) )$)). For the water level under the bedrock, constant infiltration is used: "
113 "(( $ \\delta_d(q_d^N + \\sigma_d(H_d^D - H_d^S) )$)). Parameters: ``bc_pressure``, ``bc_switch_pressure``, "
114 " ``bc_sigma``, ``bc_flux``."
125 "Boundary piezometric head for BC types: dirichlet, robin, and river." )
127 "Boundary switch piezometric head for BC types: seepage, river." )
129 "Initial condition for the pressure given as the piezometric head." )
131 return field_descriptor;
137 "Linear solver for MH problem.")
139 "Residual tolerance.")
141 "Minimum number of iterations (linear solutions) to use.\nThis is usefull if the convergence criteria "
142 "does not characterize your goal well enough so it converges prematurely, possibly even without a single linear solution."
143 "If greater then 'max_it' the value is set to 'max_it'.")
145 "Maximum number of iterations (linear solutions) of the non-linear solver.")
147 "If a stagnation of the nonlinear solver is detected the solver stops. "
148 "A divergence is reported by default, forcing the end of the simulation. By setting this flag to 'true', the solver "
149 "ends with convergence success on stagnation, but it reports warning about it.")
158 "Vector of the gravity force. Dimensionless.")
163 "Input data for Darcy flow model.")
165 "Non-linear solver for MH problem.")
167 "Output stream settings.\n Specify file format, precision etc.")
170 IT::Default(
"{ \"fields\": [ \"pressure_p0\", \"velocity_p0\" ] }"),
171 "Specification of output fields and output times.")
173 "Output settings specific to Darcy flow model.\n"
174 "Includes raw output and some experimental functionality.")
176 "Settings for computing mass balance.")
178 "Number of Schur complements to perform when solving MH system.")
180 "Method for coupling Darcy flow between dimensions on incompatible meshes. [Experimental]" )
186 Input::register_class< DarcyMH, Mesh &, const Input::Record >(
equation_name()) +
193 *
this += field_ele_pressure.name(
"pressure_p0")
196 .description(
"Pressure solution - P0 interpolation.");
198 *
this += field_edge_pressure.name(
"pressure_edge")
201 .description(
"Pressure solution - Crouzeix-Raviart interpolation.");
203 *
this += field_ele_piezo_head.name(
"piezo_head_p0")
206 .description(
"Piezo head solution - P0 interpolation.");
208 *
this += field_ele_velocity.name(
"velocity_p0")
209 .units(
UnitSI().m().s(-1))
211 .description(
"Velocity solution - P0 interpolation.");
213 *
this += flux.name(
"flux")
214 .units(
UnitSI().m().s(-1))
216 .description(
"Darcy flow flux.");
218 *
this += anisotropy.name(
"anisotropy")
219 .description(
"Anisotropy of the conductivity tensor.")
220 .input_default(
"1.0")
223 *
this += cross_section.name(
"cross_section")
224 .description(
"Complement dimension parameter (cross section for 1D, thickness for 2D).")
225 .input_default(
"1.0")
226 .units(
UnitSI().m(3).md() );
228 *
this += conductivity.name(
"conductivity")
229 .description(
"Isotropic conductivity scalar.")
230 .input_default(
"1.0")
231 .units(
UnitSI().m().s(-1) )
234 *
this += sigma.name(
"sigma")
235 .description(
"Transition coefficient between dimensions.")
236 .input_default(
"1.0")
239 *
this += water_source_density.name(
"water_source_density")
240 .description(
"Water source density.")
241 .input_default(
"0.0")
244 *
this += bc_type.name(
"bc_type")
245 .description(
"Boundary condition type.")
246 .input_selection( get_bc_type_selection() )
247 .input_default(
"\"none\"")
251 .disable_where(bc_type, {
none, seepage} )
253 .description(
"Prescribed pressure value on the boundary. Used for all values of ``bc_type`` except ``none`` and ``seepage``. "
254 "See documentation of ``bc_type`` for exact meaning of ``bc_pressure`` in individual boundary condition types.")
255 .input_default(
"0.0")
259 .disable_where(bc_type, {
none, dirichlet} )
261 .description(
"Incoming water boundary flux. Used for bc_types : ``total_flux``, ``seepage``, ``river``.")
262 .input_default(
"0.0")
263 .units(
UnitSI().m().s(-1) );
265 *
this += bc_robin_sigma
266 .disable_where(bc_type, {
none, dirichlet, seepage} )
267 .name(
"bc_robin_sigma")
268 .description(
"Conductivity coefficient in the ``total_flux`` or the ``river`` boundary condition type.")
269 .input_default(
"0.0")
272 *
this += bc_switch_pressure
273 .disable_where(bc_type, {
none, dirichlet, total_flux} )
274 .name(
"bc_switch_pressure")
275 .description(
"Critical switch pressure for ``seepage`` and ``river`` boundary conditions.")
276 .input_default(
"0.0")
281 *
this += init_pressure.name(
"init_pressure")
282 .description(
"Initial condition for pressure in time dependent problems.")
283 .input_default(
"0.0")
286 *
this += storativity.name(
"storativity")
287 .description(
"Storativity (in time dependent problems).")
288 .input_default(
"0.0")
291 *
this += extra_storativity.name(
"extra_storativity")
292 .description(
"Storativity added from upstream equation.")
294 .input_default(
"0.0")
295 .flags( input_copy );
297 *
this += extra_source.name(
"extra_water_source_density")
298 .description(
"Water source density added from upstream equation.")
299 .input_default(
"0.0")
301 .flags( input_copy );
303 *
this += gravity_field.name(
"gravity")
304 .description(
"Gravity vector.")
305 .input_default(
"0.0")
308 *
this += bc_gravity.name(
"bc_gravity")
309 .description(
"Boundary gravity vector.")
310 .input_default(
"0.0")
313 *
this += init_piezo_head.name(
"init_piezo_head")
315 .input_default(
"0.0")
316 .description(
"Init piezo head.");
318 *
this += bc_piezo_head.name(
"bc_piezo_head")
320 .input_default(
"0.0")
321 .description(
"Boundary piezo head.");
323 *
this += bc_switch_piezo_head.name(
"bc_switch_piezo_head")
325 .input_default(
"0.0")
326 .description(
"Boundary switch piezo head.");
376 MessageOut() <<
"Duplicate key 'time', time in flow equation is already initialized from parent class!";
416 return next_t * (1 - 2*std::numeric_limits<double>::epsilon());
438 gravity_array.copy_to(gvec);
444 auto field_algo=std::make_shared<FieldConstant<3, FieldValue<3>::VectorFixed>>();
445 field_algo->set_value(gvalue);
446 eq_fields_->gravity_field.set(field_algo, 0.0);
471 MessageOut() <<
"Missing the key 'time', obligatory for the transient problems." << endl;
484 std::shared_ptr< FiniteElement<1> > fe1_rt = std::make_shared<FE_RT0_disc<1>>();
485 std::shared_ptr< FiniteElement<2> > fe2_rt = std::make_shared<FE_RT0_disc<2>>();
486 std::shared_ptr< FiniteElement<3> > fe3_rt = std::make_shared<FE_RT0_disc<3>>();
487 std::shared_ptr< FiniteElement<0> > fe0_disc = std::make_shared<FE_P_disc<0>>(0);
488 std::shared_ptr< FiniteElement<1> > fe1_disc = std::make_shared<FE_P_disc<1>>(0);
489 std::shared_ptr< FiniteElement<2> > fe2_disc = std::make_shared<FE_P_disc<2>>(0);
490 std::shared_ptr< FiniteElement<3> > fe3_disc = std::make_shared<FE_P_disc<3>>(0);
491 std::shared_ptr< FiniteElement<0> > fe0_cr = std::make_shared<FE_CR<0>>();
492 std::shared_ptr< FiniteElement<1> > fe1_cr = std::make_shared<FE_CR<1>>();
493 std::shared_ptr< FiniteElement<2> > fe2_cr = std::make_shared<FE_CR<2>>();
494 std::shared_ptr< FiniteElement<3> > fe3_cr = std::make_shared<FE_CR<3>>();
496 FESystem<0> fe0_sys( {fe0_disc, fe0_disc, fe0_cr} );
502 std::shared_ptr<DiscreteSpace> ds = std::make_shared<EqualOrderDiscreteSpace>(
mesh_, fe_sys);
503 eq_data_->dh_ = std::make_shared<DOFHandlerMultiDim>(*
mesh_);
514 uint rt_component = 0;
516 auto ele_flux_ptr = create_field_fe<3, FieldValue<3>::VectorFixed>(
eq_data_->dh_, &
eq_data_->full_solution, rt_component);
521 uint p_ele_component = 1;
522 auto ele_pressure_ptr = create_field_fe<3, FieldValue<3>::Scalar>(
eq_data_->dh_, &
eq_data_->full_solution, p_ele_component);
523 eq_fields_->field_ele_pressure.set(ele_pressure_ptr, 0.0);
525 uint p_edge_component = 2;
526 auto edge_pressure_ptr = create_field_fe<3, FieldValue<3>::Scalar>(
eq_data_->dh_, &
eq_data_->full_solution, p_edge_component);
527 eq_fields_->field_edge_pressure.set(edge_pressure_ptr, 0.0);
536 uint p_edge_component = 2;
537 eq_data_->dh_cr_ = std::make_shared<SubDOFHandlerMultiDim>(
eq_data_->dh_,p_edge_component);
542 std::shared_ptr<DiscreteSpace> ds_cr_disc = std::make_shared<EqualOrderDiscreteSpace>(
mesh_, fe_cr_disc);
543 eq_data_->dh_cr_disc_ = std::make_shared<DOFHandlerMultiDim>(*
mesh_);
544 eq_data_->dh_cr_disc_->distribute_dofs(ds_cr_disc);
554 .val<Input::AbstractRecord>(
"linear_solver");
605 if (zero_time_term_from_right) {
606 MessageOut() <<
"Flow zero time step - steady case\n";
615 MessageOut() <<
"Flow zero time step - unsteady case\n";
627 eq_data_->full_solution.local_to_ghost_begin();
628 eq_data_->full_solution.local_to_ghost_end();
645 eq_data_->full_solution.local_to_ghost_begin();
646 eq_data_->full_solution.local_to_ghost_end();
657 bool jump_time =
eq_fields_->storativity.is_jump_time();
658 if (! zero_time_term_from_left) {
659 MessageOut() <<
"Flow time step - unsteady case\n";
669 WarningOut() <<
"Output of solution discontinuous in time not supported yet.\n";
678 if (! zero_time_term_from_left && ! jump_time && output)
output_data();
685 if (zero_time_term_from_right) {
686 MessageOut() <<
"Flow time step - steady case\n";
691 }
else if (! zero_time_term_from_left && jump_time) {
692 WarningOut() <<
"Discontinuous time term not supported yet.\n";
704 return (
eq_fields_->storativity.input_list_size() == 0);
717 MessageOut().fmt(
"[nonlinear solver] norm of initial residual: {}\n", residual_norm);
720 int is_linear_common;
725 this->
max_n_it_ = nl_solver_rec.
val<
unsigned int>(
"max_it");
726 this->
min_n_it_ = nl_solver_rec.
val<
unsigned int>(
"min_it");
727 if (this->min_n_it_ > this->max_n_it_) this->min_n_it_ = this->
max_n_it_;
729 if (! is_linear_common) {
737 while (nonlinear_iteration_ < this->
min_n_it_ ||
738 (residual_norm > this->tolerance_ && nonlinear_iteration_ < this->
max_n_it_ )) {
740 convergence_history.push_back(residual_norm);
743 if (convergence_history.size() >= 5 &&
744 convergence_history[ convergence_history.size() - 1]/convergence_history[ convergence_history.size() - 2] > 0.9 &&
745 convergence_history[ convergence_history.size() - 1]/convergence_history[ convergence_history.size() - 5] > 0.8) {
751 THROW(ExcSolverDiverge() << EI_Reason(
"Stagnation."));
755 if (! is_linear_common)
761 if (is_linear_common){
764 MessageOut().fmt(
"[nonlinear solver] lin. it: {}, reason: {}, residual: {}\n",
778 MessageOut().fmt(
"[nonlinear solver] it: {} lin. it: {}, reason: {}, residual: {}\n",
781 chkerr(VecDestroy(&save_solution));
805 if(
eq_data_->mortar_method_ != MortarMethod::NoMortar){
806 auto multidim_assembler = AssemblyFlowBase::create< AssemblyMH >(
eq_fields_,
eq_data_);
808 unsigned int dim = dh_cell.dim();
809 multidim_assembler[dim-1]->fix_velocity(dh_cell);
876 unsigned int dim = dh_cell.dim();
877 assembler[dim-1]->assemble(dh_cell);
895 double zeros[100000];
896 for(
int i=0; i<100000; i++) zeros[i] = 0.0;
899 tmp_rows.reserve(200);
902 dofs.reserve(
eq_data_->dh_->max_elem_dofs());
903 dofs_ngh.reserve(
eq_data_->dh_->max_elem_dofs());
908 const uint ndofs = dh_cell.n_dofs();
910 dh_cell.get_dof_indices(dofs);
919 for (
DHCellSide neighb_side : dh_cell.neighb_sides() ) {
927 const uint ndofs_ngh = dh_neighb_cell.
n_dofs();
928 dofs_ngh.resize(ndofs_ngh);
934 const unsigned int t = dh_neighb_cell.
n_dofs() - (dh_neighb_cell.
dim()+1) + neighb_side.side().side_idx();
935 tmp_rows.push_back(dofs_ngh[t]);
949 for(
auto &isec : isec_list ) {
953 const uint ndofs_slave = dh_cell_slave.
n_dofs();
954 dofs_ngh.resize(ndofs_slave);
958 for(
unsigned int i_side=0; i_side < dh_cell_slave.
elm()->
n_sides(); i_side++) {
960 tmp_rows.push_back( dofs_ngh[(ndofs_slave+1)/2+i_side] );
970 edge_rows = dofs.
data() + nsides +1;
1009 const uint ndofs = dh_cell.n_dofs();
1010 global_dofs.resize(ndofs);
1011 dh_cell.get_dof_indices(global_dofs);
1016 double source = ele.
measure() * cs *
1023 {dh_cell.get_loc_dof_indices()[ndofs/2]}, {0}, {source});
1043 #ifdef FLOW123D_HAVE_BDDCML
1044 WarningOut() <<
"For BDDC no Schur complements are used.";
1057 THROW( ExcBddcmlNotSupported() );
1063 WarningOut() <<
"Invalid number of Schur Complements. Using 2.";
1076 ls->LinSys::set_from_input(in_rec);
1085 ISCreateGeneral(PETSC_COMM_SELF, side_dofs_vec.size(), &(side_dofs_vec[0]), PETSC_COPY_VALUES, &is);
1100 const PetscInt *b_indices;
1101 ISGetIndices(ls->
IsB, &b_indices);
1103 for(
uint i_b=0, i_bb=0; i_b < b_size && i_bb < elem_dofs_vec.size(); i_b++) {
1104 if (b_indices[i_b] == elem_dofs_vec[i_bb])
1105 elem_dofs_vec[i_bb++] = i_b + ds->
begin();
1107 ISRestoreIndices(ls->
IsB, &b_indices);
1110 ISCreateGeneral(PETSC_COMM_SELF, elem_dofs_vec.size(), &(elem_dofs_vec[0]), PETSC_COPY_VALUES, &is);
1138 THROW( ExcUnknownSolver() );
1150 START_TIMER(
"DarcyFlowMH_Steady::assembly_linear_system");
1207 std::string output_file;
1219 PetscViewerASCIIOpen(PETSC_COMM_WORLD, output_file.c_str(), &viewer);
1220 PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
1232 double d_max = std::numeric_limits<double>::max();
1233 double h1 = d_max, h2 = d_max, h3 = d_max;
1234 double he2 = d_max, he3 = d_max;
1237 case 1: h1 = std::min(h1,ele.measure());
break;
1238 case 2: h2 = std::min(h2,ele.measure());
break;
1239 case 3: h3 = std::min(h3,ele.measure());
break;
1242 for (
unsigned int j=0; j<ele->n_sides(); j++) {
1244 case 2: he2 = std::min(he2, ele.side(j)->measure());
break;
1245 case 3: he3 = std::min(he3, ele.side(j)->measure());
break;
1249 if(h1 == d_max) h1 = 0;
1250 if(h2 == d_max) h2 = 0;
1251 if(h3 == d_max) h3 = 0;
1252 if(he2 == d_max) he2 = 0;
1253 if(he3 == d_max) he3 = 0;
1256 file = fopen(output_file.c_str(),
"a");
1258 fprintf(file,
"nB = %d;\n",
eq_data_->dh_->mesh()->get_el_ds()->size());
1260 fprintf(file,
"h1 = %e;\nh2 = %e;\nh3 = %e;\n", h1, h2, h3);
1261 fprintf(file,
"he2 = %e;\nhe3 = %e;\n", he2, he3);
1279 START_TIMER(
"DarcyFlowMH_Steady::set_mesh_data_for_bddc");
1307 dh_cell.get_dof_indices(cell_dofs_global);
1309 inet.insert(inet.end(), cell_dofs_global.begin(), cell_dofs_global.end());
1310 uint n_inet = cell_dofs_global.size();
1313 uint dim = dh_cell.elm().dim();
1314 elDimMax = std::max( elDimMax, dim );
1315 elDimMin = std::min( elDimMin, dim );
1319 isegn.push_back( dh_cell.elm_idx());
1322 for (
unsigned int si=0; si<dh_cell.elm()->n_sides(); si++) {
1323 arma::vec3 coord = dh_cell.elm().side(si)->centre();
1326 localDofMap.insert( std::make_pair( cell_dofs_global[si], coord ) );
1328 localDofMap.insert( std::make_pair( cell_dofs_global[si+dim+2], coord ) );
1331 arma::vec3 elm_centre = dh_cell.elm().centre();
1332 localDofMap.insert( std::make_pair( cell_dofs_global[dim+1], elm_centre ) );
1336 for(
DHCellSide side : dh_cell.neighb_sides()) {
1337 uint neigh_dim = side.cell().elm().dim();
1338 side.cell().get_dof_indices(cell_dofs_global);
1339 int edge_row = cell_dofs_global[neigh_dim+2+side.side_idx()];
1340 localDofMap.insert( std::make_pair( edge_row, side.centre() ) );
1341 inet.push_back( edge_row );
1344 nnet.push_back(n_inet);
1349 double conduct =
eq_fields_->conductivity.value( elm_centre , dh_cell.elm() );
1350 auto aniso =
eq_fields_->anisotropy.value( elm_centre , dh_cell.elm() );
1354 for (
int i = 0; i < 3; i++) {
1355 coef = coef + aniso.at(i,i);
1358 coef = conduct*coef / 3;
1360 ASSERT_GT( coef, 0. ).error(
"Zero coefficient of hydrodynamic resistance. \n");
1361 element_permeability.push_back( 1. / coef );
1371 auto distr =
eq_data_->dh_->distr();
1380 int numNodeSub = localDofMap.size();
1391 for ( ; itB != localDofMap.end(); ++itB ) {
1392 isngn[ind] = itB -> first;
1395 for (
int j = 0; j < 3; j++ ) {
1396 xyz[ j*numNodeSub + ind ] = coord[j];
1401 localDofMap.clear();
1409 Global2LocalMap_ global2LocalNodeMap;
1410 for (
unsigned ind = 0; ind < isngn.size(); ++ind ) {
1411 global2LocalNodeMap.insert( std::make_pair(
static_cast<unsigned>( isngn[ind] ), ind ) );
1416 for (
unsigned int iEle = 0; iEle < isegn.size(); iEle++ ) {
1417 int nne = nnet[ iEle ];
1418 for (
int ien = 0; ien < nne; ien++ ) {
1420 int indGlob = inet[indInet];
1422 Global2LocalMap_::iterator pos = global2LocalNodeMap.find( indGlob );
1423 ASSERT( pos != global2LocalNodeMap.end() )(indGlob).error(
"Cannot remap node global index to local indices.\n");
1424 int indLoc =
static_cast<int> ( pos -> second );
1427 inet[ indInet++ ] = indLoc;
1431 int numNodes =
size;
1432 int numDofsInt =
size;
1434 int meshDim = elDimMax;
1446 bddc_ls -> load_mesh( LinSys_BDDC::BDDCMatrixType::SPD_VIA_SYMMETRICGENERAL, spaceDim, numNodes, numDofsInt, inet, nnet, nndf, isegn, isngn, isngn, xyz, element_permeability, meshDim );
1468 if(
time_ !=
nullptr)
1518 LocDofVec loc_dofs = dh_cell.get_loc_dof_indices();
1519 for (
unsigned int i=0; i<ele->
n_sides(); i++) {
1523 eq_data_->full_solution.add(loc_dofs[edge_dof], init_value/n_sides_of_edge);
1527 eq_data_->full_solution.ghost_to_local_begin();
1528 eq_data_->full_solution.ghost_to_local_end();
1529 eq_data_->full_solution.local_to_ghost_begin();
1530 eq_data_->full_solution.local_to_ghost_end();
1536 START_TIMER(
"DarcyFlowMH::reconstruct_solution_from_schur");
1539 unsigned int dim = dh_cell.dim();
1540 assembler[dim-1]->assemble_reconstruct(dh_cell);
1543 eq_data_->full_solution.local_to_ghost_begin();
1544 eq_data_->full_solution.local_to_ghost_end();
1555 PetscScalar *local_diagonal;
1563 dofs.reserve(
eq_data_->dh_->max_elem_dofs());
1567 dofs.resize(dh_cell.n_dofs());
1568 dh_cell.get_dof_indices(dofs);
1571 const uint p_ele_dof = dh_cell.n_dofs() / 2;
1577 local_diagonal[dh_cell.get_loc_dof_indices()[p_ele_dof]]= - diagonal_coeff /
time_->
dt();
1580 {dh_cell.get_loc_dof_indices()[p_ele_dof]},
1581 {diagonal_coeff}, 0.0);
1583 VecRestoreArray(new_diagonal,& local_diagonal);
1584 MatDiagonalSet(*( schur0->get_matrix() ), new_diagonal, ADD_VALUES);
1586 schur0->set_matrix_changed();
1588 balance_->finish_mass_assembly(eq_data_->water_balance_idx);
1595 if (scale_factor != 1.0) {
1616 void dofs_range(
unsigned int n_dofs,
unsigned int &min,
unsigned int &max,
unsigned int component) {
1620 }
else if (component==1) {
1631 ASSERT_LT(component, 3).error(
"Invalid component!");
1632 unsigned int i, n_dofs, min, max;
1636 n_dofs = dh_cell.get_dof_indices(dof_indices);
1638 for (i=min; i<max; ++i) dof_vec.push_back(dof_indices[i]);
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.
#define ASSERT_LT(a, b)
Definition of comparative assert macro (Less Than) only for debug mode.
#define ASSERT_EQ(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
#define ASSERT_GT(a, b)
Definition of comparative assert macro (Greater Than) only for debug mode.
static const Input::Type::Record & get_input_type()
Main balance input record type.
Cell accessor allow iterate over DOF handler cells.
const ElementAccessor< 3 > elm() const
Return ElementAccessor to element of loc_ele_idx_.
unsigned int n_dofs() const
Return number of dofs on given cell.
unsigned int get_dof_indices(std::vector< LongIdx > &indices) const
Fill vector of the global indices of dofs associated to the cell.
unsigned int dim() const
Return dimension of element appropriate to cell.
Side accessor allows to iterate over sides of DOF handler cell.
static Input::Type::Abstract & get_input_type()
MortarMethod
Type of experimental Mortar-like method for non-compatible 1d-2d interaction.
static const Input::Type::Instance & get_input_type_specific()
void output()
Calculate values for output.
static const Input::Type::Instance & get_input_type(FieldSet &eq_data, const std::string &equation_name)
EqFields()
Creation of all fields.
static const Input::Type::Selection & get_bc_type_selection()
Return a Selection corresponding to enum BC_Type.
DarcyFlowMHOutput * output_object
void set_mesh_data_for_bddc(LinSys_BDDC *bddc_ls)
virtual double solved_time() override
void prepare_new_time_step()
void assembly_mh_matrix(MultidimAssembly &assembler)
static const Input::Type::Record & type_field_descriptor()
virtual void read_initial_condition()
void print_matlab_matrix(string matlab_file)
Print darcy flow matrix in matlab format into a file.
virtual void initialize_specific()
std::shared_ptr< Balance > balance_
std::vector< int > get_component_indices_vec(unsigned int component) const
Get vector of all DOF indices of given component (0..side, 1..element, 2..edge)
void allocate_mh_matrix()
void initialize() override
static std::string equation_name()
void solve_nonlinear()
Solve method common to zero_time_step and update solution.
virtual ~DarcyMH() override
std::shared_ptr< EqData > eq_data_
bool use_steady_assembly_
void zero_time_step() override
std::shared_ptr< EqFields > eq_fields_
static const int registrar
Registrar of class to factory.
virtual double solution_precision() const
unsigned int nonlinear_iteration_
virtual void output_data() override
Write computed fields.
virtual void assembly_linear_system()
static const Input::Type::Selection & get_mh_mortar_selection()
Selection for enum MortarMethod.
friend class DarcyFlowMHOutput
void reconstruct_solution_from_schur(MultidimAssembly &assembler)
virtual bool zero_time_term(bool time_global=false)
DarcyMH(Mesh &mesh, const Input::Record in_rec, TimeGovernor *tm=nullptr)
CREATE AND FILL GLOBAL MH MATRIX OF THE WATER MODEL.
virtual void postprocess()
postprocess velocity field (add sources)
static const Input::Type::Record & get_input_type()
void update_solution() override
void create_linear_system(Input::AbstractRecord rec)
void solve_time_step(bool output=true)
Solve the problem without moving to next time and without output.
virtual void setup_time_term()
virtual void assembly_source_term()
Source term is implemented differently in LMH version.
unsigned int begin(int proc) const
get starting local index
unsigned int n_sides() const
Returns number of sides aligned with the edge.
double measure() const
Computes the measure of the element.
SideIter side(const unsigned int loc_index)
unsigned int idx() const
We need this method after replacing Region by RegionIdx, and movinf RegionDB instance into particular...
arma::vec::fixed< spacedim > centre() const
Computes the barycenter.
unsigned int n_neighs_vb() const
Return number of neighbours.
unsigned int n_sides() const
Input::Record input_record_
static Input::Type::Record & record_template()
Template Record with common keys for derived equations.
std::shared_ptr< FieldSet > eq_fieldset_
Compound finite element on dim dimensional simplex.
static const std::string field_descriptor_record_description(const string &record_name)
static constexpr Mask input_copy
static constexpr Mask equation_result
Match result fields. These are never given by input or copy of input.
Dedicated class for storing path to input and output files.
Common base for intersection object.
unsigned int bulk_ele_idx() const
Returns index of bulk element.
static const Input::Type::Record & get_input_type()
void set_from_input(const Input::Record in_rec) override
void set_from_input(const Input::Record in_rec) override
static const Input::Type::Record & get_input_type()
const Vec & get_solution()
virtual const Vec * get_rhs()
void set_negative_definite(bool flag=true)
void set_solution(Vec sol_vec)
void set_matrix_changed()
virtual double get_solution_precision()=0
virtual void start_add_assembly()
virtual void finish_assembly()=0
void rhs_set_value(int row, double val)
virtual const Mat * get_matrix()
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)
virtual PetscErrorCode rhs_zero_entries()
virtual void start_allocation()
virtual SolveInfo solve()=0
void set_positive_definite(bool flag=true)
virtual PetscErrorCode mat_zero_entries()
static Input::Type::Abstract & get_input_type()
const RegionDB & region_db() const
unsigned int n_edges() const
unsigned int n_elements() const
Range< ElementAccessor< 3 > > elements_range() const
Returns range of mesh elements.
BCMesh * bc_mesh() const override
Implement MeshBase::bc_mesh(), getter of boundary mesh.
MixedMeshIntersections & mixed_intersections()
unsigned int n_sides() const
std::vector< std::vector< ILpair > > element_intersections_
static const Input::Type::Record & get_input_type()
The specification of output stream.
RegionSet get_region_set(const std::string &set_name) const
unsigned int bulk_idx() const
Returns index of the region in the bulk set.
Distribution * make_complement_distribution()
get distribution of complement object if complement is defined
void set_from_input(const Input::Record in_rec) override
void set_complement(LinSys_PETSC *ls)
Set complement LinSys object.
Edge edge() const
Returns pointer to the edge connected to the side.
Basic time management functionality for unsteady (and steady) solvers (class Equation).
bool is_end() const
Returns true if the actual time is greater than or equal to the end time.
int set_upper_constraint(double upper, std::string message)
Sets upper constraint for the next time step estimating.
bool is_changed_dt() const
void view(const char *name="") const
double 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.
unsigned int index() const
Class for representation SI units of Fields.
static UnitSI & dimensionless()
Returns dimensionless unit.
void dofs_range(unsigned int n_dofs, unsigned int &min, unsigned int &max, unsigned int component)
Helper method fills range (min and max) of given component.
Output class for darcy_flow_mh model.
Support classes for parallel programing.
Definitions of basic Lagrangean finite elements with polynomial shape functions.
#define FLOW123D_FORCE_LINK_IN_CHILD(x)
#define THROW(whole_exception_expr)
Wrapper for throw. Saves the throwing point.
arma::Col< IntIdx > LocDofVec
int LongIdx
Define type that represents indices of large arrays (elements, nodes, dofs etc.)
Classes with algorithms for computation of intersections of meshes.
Wrappers for linear systems based on MPIAIJ and MATIS format.
Solver based on Multilevel BDDC - using corresponding class of OpenFTL package.
Solver based on the original PETSc solver using MPIAIJ matrix and succesive Schur complement construc...
#define WarningOut()
Macro defining 'warning' record of log.
#define MessageOut()
Macro defining 'message' record of log.
#define DebugOut()
Macro defining 'debug' record of log.
#define MPI_Allreduce(sendbuf, recvbuf, count, datatype, op, comm)
FMT_FUNC int fprintf(std::ostream &os, CStringRef format, ArgList args)
Implementation of range helper class.
Assembly explicit Schur complement for the given linear system. Provides method for resolution of the...
SchurComplement SchurComplement
#define END_TIMER(tag)
Ends a timer with specified tag.
#define START_TIMER(tag)
Starts a timer with specified tag.
void chkerr(unsigned int ierr)
Replacement of new/delete operator in the spirit of xmalloc.
Basic time management class.