53 return Selection(
"DG_variant",
"Type of penalty term.")
54 .
add_value(non_symmetric,
"non-symmetric",
"non-symmetric weighted interior penalty DG method")
55 .
add_value(incomplete,
"incomplete",
"incomplete weighted interior penalty DG method")
56 .
add_value(symmetric,
"symmetric",
"symmetric weighted interior penalty DG method")
79 std::string equation_name = std::string(Model::ModelEqData::name()) +
"_DG";
80 return Model::get_input_type(
"DG",
"Discontinuous Galerkin (DG) solver")
83 "Solver for the linear system.")
86 .make_field_descriptor_type(equation_name)),
88 "Input fields of the equation.")
90 "Variant of the interior penalty discontinuous Galerkin method.")
92 "Polynomial order for the finite element in DG method (order 0 is suitable if there is no diffusion/dispersion).")
94 "If true, use DG projection of the initial condition field."
95 "Otherwise, evaluate initial condition field directly (well suited for reading native data).")
97 EqFields().output_fields.make_output_type(equation_name,
""),
98 IT::Default(
"{ \"fields\": [ " + Model::ModelEqData::default_output_field() +
"] }"),
99 "Specification of output fields and output times.")
103 template<
class Model>
105 Input::register_class< TransportDG<Model>,
Mesh &,
const Input::Record>(std::string(Model::ModelEqData::name()) +
"_DG") +
110 template<
class Model>
114 .
name(
"fracture_sigma")
116 "Coefficient of diffusive transfer through fractures (for each substance).")
124 "Penalty parameter influencing the discontinuity of the solution (for each substance). "
125 "Its default value 1 is sufficient in most cases. Higher value diminishes the inter-element jumps.")
138 .
description(
"Subdomain ids of the domain decomposition.");
144 this->add_coords_field();
145 this->set_default_fieldset();
150 template<
typename Model>
152 :
Model(init_mesh, in_rec),
166 Model::init_balance(in_rec);
177 eq_data_->dg_order = in_rec.
val<
unsigned int>(
"dg_order");
179 Model::init_from_input(in_rec);
182 shared_ptr<DiscreteSpace> ds = make_shared<EqualOrderDiscreteSpace>(Model::mesh_, fe);
183 eq_data_->dh_ = make_shared<DOFHandlerMultiDim>(*Model::mesh_);
190 template<
class Model>
193 eq_fields_->set_components(eq_data_->substances_.names());
194 eq_fields_->set_input_list( input_rec.val<
Input::Array>(
"input_fields"), *(Model::time_) );
195 eq_data_->set_time_governor(Model::time_);
196 eq_data_->balance_ = this->balance();
197 eq_fields_->initialize();
200 eq_data_->max_edg_sides = max(Model::mesh_->max_edge_sides(1), max(Model::mesh_->max_edge_sides(2), Model::mesh_->max_edge_sides(3)));
201 ret_sources.resize(eq_data_->n_substances());
202 ret_sources_prev.resize(eq_data_->n_substances());
205 if (input_rec.opt_val(
"user_fields", user_fields_arr)) {
206 this->init_user_fields(user_fields_arr, eq_fields_->output_fields);
209 eq_data_->output_vec.resize(eq_data_->n_substances());
210 eq_fields_->output_field.set_components(eq_data_->substances_.names());
211 eq_fields_->output_field.set_mesh(*Model::mesh_);
212 eq_fields_->output_fields.set_mesh(*Model::mesh_);
215 eq_fields_->output_field.setup_components();
216 for (
unsigned int sbi=0; sbi<eq_data_->n_substances(); sbi++)
219 auto output_field_ptr = create_field_fe< 3, FieldValue<3>::Scalar >(eq_data_->dh_);
220 eq_fields_->output_field[sbi].set(output_field_ptr, 0);
221 eq_data_->output_vec[sbi] = output_field_ptr->vec();
225 eq_fields_->output_fields.initialize(Model::output_stream_, Model::mesh_, input_rec.val<
Input::Record>(
"output"), this->time());
228 std::string petsc_default_opts;
229 if (eq_data_->dh_->distr()->np() == 1)
230 petsc_default_opts =
"-ksp_type bcgs -pc_type ilu -pc_factor_levels 2 -ksp_diagonal_scale_fix -pc_factor_fill 6.0";
232 petsc_default_opts =
"-ksp_type bcgs -ksp_diagonal_scale_fix -pc_type asm -pc_asm_overlap 4 -sub_pc_type ilu -sub_pc_factor_levels 3 -sub_pc_factor_fill 6.0";
235 eq_data_->ls =
new LinSys*[eq_data_->n_substances()];
236 eq_data_->ls_dt =
new LinSys*[eq_data_->n_substances()];
237 eq_data_->conc_fe.resize(eq_data_->n_substances());
240 shared_ptr<DiscreteSpace> ds = make_shared<EqualOrderDiscreteSpace>(Model::mesh_, fe);
241 eq_data_->dh_p0 = make_shared<DOFHandlerMultiDim>(*Model::mesh_);
242 eq_data_->dh_p0->distribute_dofs(ds);
244 stiffness_matrix.resize(eq_data_->n_substances(),
nullptr);
245 mass_matrix.resize(eq_data_->n_substances(),
nullptr);
246 rhs.resize(eq_data_->n_substances(),
nullptr);
247 mass_vec.resize(eq_data_->n_substances(),
nullptr);
248 eq_data_->ret_vec.resize(eq_data_->n_substances(),
nullptr);
250 for (
unsigned int sbi = 0; sbi < eq_data_->n_substances(); sbi++) {
251 eq_data_->ls[sbi] =
new LinSys_PETSC(eq_data_->dh_->distr().get(), petsc_default_opts);
253 eq_data_->ls[sbi]->set_solution(eq_data_->output_vec[sbi].petsc_vec());
255 eq_data_->ls_dt[sbi] =
new LinSys_PETSC(eq_data_->dh_->distr().get(), petsc_default_opts);
258 eq_data_->conc_fe[sbi] = create_field_fe< 3, FieldValue<3>::Scalar >(eq_data_->dh_p0);
260 VecDuplicate(eq_data_->ls[sbi]->get_solution(), &eq_data_->ret_vec[sbi]);
264 init_projection = input_rec.val<
bool>(
"init_projection");
278 Model::balance_->allocate(eq_data_->dh_->distr()->lsize(), mass_assembly_->eval_points()->max_size());
280 eq_fields_->init_condition.setup_components();
281 for (
unsigned int sbi=0; sbi<eq_data_->n_substances(); sbi++)
288 template<
class Model>
293 if (rhs.size() > 0) {
296 for (
unsigned int i=0; i<eq_data_->n_substances(); i++)
298 if (eq_data_->ls !=
nullptr) {
299 delete eq_data_->ls[i];
300 delete eq_data_->ls_dt[i];
303 if (stiffness_matrix.size() > 0) {
304 if (stiffness_matrix[i])
305 chkerr(MatDestroy(&stiffness_matrix[i]));
307 chkerr(MatDestroy(&mass_matrix[i]));
309 chkerr(VecDestroy(&rhs[i]));
311 chkerr(VecDestroy(&mass_vec[i]));
312 if (eq_data_->ret_vec[i])
313 chkerr(VecDestroy(&eq_data_->ret_vec[i]));
316 if (eq_data_->ls !=
nullptr) {
317 delete[] eq_data_->ls;
318 delete[] eq_data_->ls_dt;
319 eq_data_->ls =
nullptr;
327 if (mass_assembly_ !=
nullptr) {
328 delete mass_assembly_;
329 delete stiffness_assembly_;
330 delete sources_assembly_;
331 delete bdr_cond_assembly_;
332 delete init_assembly_;
340 template<
class Model>
344 eq_fields_->mark_input_times( *(Model::time_) );
346 std::stringstream ss;
353 set_initial_condition();
354 for (
unsigned int sbi = 0; sbi < eq_data_->n_substances(); sbi++)
355 ( (
LinSys_PETSC *)eq_data_->ls[sbi] )->set_initial_guess_nonzero();
358 if (!allocation_done) preallocate();
361 for (
unsigned int sbi=0; sbi<eq_data_->n_substances(); ++sbi)
363 Model::balance_->calculate_instant(eq_data_->subst_idx_[sbi], eq_data_->ls[sbi]->get_solution());
366 ret_sources_prev[sbi] = 0;
373 template<
class Model>
377 for (
unsigned int i=0; i<eq_data_->n_substances(); i++)
380 eq_data_->ls[i]->start_allocation();
381 stiffness_matrix[i] = NULL;
385 eq_data_->ls_dt[i]->start_allocation();
386 mass_matrix[i] = NULL;
387 VecZeroEntries(eq_data_->ret_vec[i]);
389 stiffness_assembly_->assemble(eq_data_->dh_);
390 mass_assembly_->assemble(eq_data_->dh_);
391 sources_assembly_->assemble(eq_data_->dh_);
392 bdr_cond_assembly_->assemble(eq_data_->dh_);
393 for (
unsigned int i=0; i<eq_data_->n_substances(); i++)
395 VecAssemblyBegin(eq_data_->ret_vec[i]);
396 VecAssemblyEnd(eq_data_->ret_vec[i]);
399 allocation_done =
true;
404 template<
class Model>
409 Model::time_->next_time();
410 Model::time_->view(
"TDG");
419 for (
unsigned int i=0; i<eq_data_->n_substances(); i++)
421 eq_data_->ls_dt[i]->start_add_assembly();
422 eq_data_->ls_dt[i]->mat_zero_entries();
423 VecZeroEntries(eq_data_->ret_vec[i]);
425 mass_assembly_->assemble(eq_data_->dh_);
426 for (
unsigned int i=0; i<eq_data_->n_substances(); i++)
428 eq_data_->ls_dt[i]->finish_assembly();
429 VecAssemblyBegin(eq_data_->ret_vec[i]);
430 VecAssemblyEnd(eq_data_->ret_vec[i]);
432 if (mass_matrix[i] == NULL)
434 VecDuplicate(eq_data_->ls[i]->get_solution(), &mass_vec[i]);
435 MatMult(*(eq_data_->ls_dt[i]->get_matrix()), eq_data_->ls[i]->get_solution(), mass_vec[i]);
436 MatConvert(*( eq_data_->ls_dt[i]->get_matrix() ), MATSAME, MAT_INITIAL_MATRIX, &mass_matrix[i]);
439 MatCopy(*( eq_data_->ls_dt[i]->get_matrix() ), mass_matrix[i], DIFFERENT_NONZERO_PATTERN);
444 if (stiffness_matrix[0] == NULL
446 || eq_fields_->flow_flux.changed())
450 for (
unsigned int i=0; i<eq_data_->n_substances(); i++)
452 eq_data_->ls[i]->start_add_assembly();
453 eq_data_->ls[i]->mat_zero_entries();
455 stiffness_assembly_->assemble(eq_data_->dh_);
456 for (
unsigned int i=0; i<eq_data_->n_substances(); i++)
458 eq_data_->ls[i]->finish_assembly();
460 if (stiffness_matrix[i] == NULL)
461 MatConvert(*( eq_data_->ls[i]->get_matrix() ), MATSAME, MAT_INITIAL_MATRIX, &stiffness_matrix[i]);
463 MatCopy(*( eq_data_->ls[i]->get_matrix() ), stiffness_matrix[i], DIFFERENT_NONZERO_PATTERN);
470 || eq_fields_->flow_flux.changed())
472 for (
unsigned int i=0; i<eq_data_->n_substances(); i++)
474 eq_data_->ls[i]->start_add_assembly();
475 eq_data_->ls[i]->rhs_zero_entries();
477 sources_assembly_->assemble(eq_data_->dh_);
478 bdr_cond_assembly_->assemble(eq_data_->dh_);
479 for (
unsigned int i=0; i<eq_data_->n_substances(); i++)
481 eq_data_->ls[i]->finish_assembly();
483 if (rhs[i] ==
nullptr) VecDuplicate(*( eq_data_->ls[i]->get_rhs() ), &rhs[i]);
484 VecCopy(*( eq_data_->ls[i]->get_rhs() ), rhs[i]);
506 for (
unsigned int i=0; i<eq_data_->n_substances(); i++)
508 MatConvert(stiffness_matrix[i], MATSAME, MAT_INITIAL_MATRIX, &m);
509 MatAXPY(m, 1./Model::time_->dt(), mass_matrix[i], SUBSET_NONZERO_PATTERN);
510 eq_data_->ls[i]->set_matrix(m, DIFFERENT_NONZERO_PATTERN);
512 VecDuplicate(rhs[i], &w);
513 VecWAXPY(w, 1./Model::time_->dt(), mass_vec[i], rhs[i]);
514 eq_data_->ls[i]->set_rhs(w);
519 eq_data_->ls[i]->solve();
522 MatMult(*(eq_data_->ls_dt[i]->get_matrix()), eq_data_->ls[i]->get_solution(), mass_vec[i]);
533 calculate_cumulative_balance();
539 template<
class Model>
543 for (
auto cell : eq_data_->dh_->own_range() )
545 LocDofVec loc_dof_indices = cell.get_loc_dof_indices();
546 unsigned int n_dofs=loc_dof_indices.n_rows;
548 DHCellAccessor dh_p0_cell = eq_data_->dh_p0->cell_accessor_from_element(cell.elm_idx());
551 for (
unsigned int sbi=0; sbi<eq_data_->n_substances(); ++sbi)
553 eq_data_->conc_fe[sbi]->vec().set(dof_p0, 0);
555 for (
unsigned int j=0; j<n_dofs; ++j)
556 eq_data_->conc_fe[sbi]->vec().add( dof_p0, eq_data_->ls[sbi]->get_solution_array()[loc_dof_indices[j]] );
558 eq_data_->conc_fe[sbi]->vec().set( dof_p0, max(eq_data_->conc_fe[sbi]->vec().get(dof_p0)/n_dofs, 0.) );
566 template<
class Model>
575 eq_fields_->output_fields.set_time( this->time().step(),
LimitSide::left);
577 eq_fields_->output_fields.output(this->time().step());
579 Model::output_data();
582 for (
unsigned int sbi=0; sbi<eq_data_->n_substances(); ++sbi)
583 Model::balance_->calculate_instant(eq_data_->subst_idx_[sbi], eq_data_->ls[sbi]->get_solution());
584 Model::balance_->output();
591 template<
class Model>
594 if (Model::balance_->cumulative())
596 for (
unsigned int sbi=0; sbi<eq_data_->n_substances(); ++sbi)
598 Model::balance_->calculate_cumulative(eq_data_->subst_idx_[sbi], eq_data_->ls[sbi]->get_solution());
601 VecDot(eq_data_->ret_vec[sbi], eq_data_->ls[sbi]->get_solution(), &ret_sources[sbi]);
603 Model::balance_->add_cumulative_source(eq_data_->subst_idx_[sbi], (ret_sources[sbi]-ret_sources_prev[sbi])/Model::time_->dt());
604 ret_sources_prev[sbi] = ret_sources[sbi];
612 template<
class Model>
617 for (
unsigned int sbi=0; sbi<eq_data_->n_substances(); sbi++)
618 eq_data_->ls[sbi]->start_allocation();
620 init_assembly_->assemble(eq_data_->dh_);
622 for (
unsigned int sbi=0; sbi<eq_data_->n_substances(); sbi++)
623 eq_data_->ls[sbi]->start_add_assembly();
625 init_assembly_->assemble(eq_data_->dh_);
627 for (
unsigned int sbi=0; sbi<eq_data_->n_substances(); sbi++)
629 eq_data_->ls[sbi]->finish_assembly();
630 eq_data_->ls[sbi]->solve();
634 init_assembly_->assemble(eq_data_->dh_);
638 template<
class Model>
641 if (solution_changed)
643 for (
auto cell : eq_data_->dh_->own_range() )
645 LocDofVec loc_dof_indices = cell.get_loc_dof_indices();
646 unsigned int n_dofs=loc_dof_indices.n_rows;
648 DHCellAccessor dh_p0_cell = eq_data_->dh_p0->cell_accessor_from_element(cell.elm_idx());
651 for (
unsigned int sbi=0; sbi<eq_data_->n_substances(); ++sbi)
653 double old_average = 0;
654 for (
unsigned int j=0; j<n_dofs; ++j)
655 old_average += eq_data_->ls[sbi]->get_solution_array()[loc_dof_indices[j]];
656 old_average /= n_dofs;
658 for (
unsigned int j=0; j<n_dofs; ++j)
659 eq_data_->ls[sbi]->get_solution_array()[loc_dof_indices[j]]
660 += eq_data_->conc_fe[sbi]->vec().get(dof_p0) - old_average;
665 for (
unsigned int sbi=0; sbi<eq_data_->n_substances(); ++sbi)
666 MatMult(*(eq_data_->ls_dt[sbi]->get_matrix()), eq_data_->ls[sbi]->get_solution(), mass_vec[sbi]);
Discontinuous Galerkin method for equation of transport with dispersion.
Cell accessor allow iterate over DOF handler cells.
LocDofVec get_loc_dof_indices() const
Returns the local indices of dofs associated to the cell on the local process.
static Input::Type::Record & user_fields_template(std::string equation_name)
Template Record with common key user_fields for derived equations.
static bool print_message_table(ostream &stream, std::string equation_name)
FieldCommon & description(const string &description)
FieldCommon & flags_add(FieldFlag::Flags::Mask mask)
FieldCommon & flags(FieldFlag::Flags::Mask mask)
FieldCommon & name(const string &name)
FieldCommon & units(const UnitSI &units)
Set basic units of the field.
FieldCommon & input_default(const string &input_default)
static constexpr Mask in_time_term
A field is part of time term of the equation.
static constexpr Mask equation_external_output
Match an output field, that can be also copy of other field.
static constexpr Mask in_rhs
A field is part of the right hand side of the equation.
static constexpr Mask in_main_matrix
A field is part of main "stiffness matrix" of the equation.
Generic class of assemblation.
static auto region_id(Mesh &mesh) -> IndexField
static auto subdomain(Mesh &mesh) -> IndexField
static const Input::Type::Record & get_input_type()
Field< 3, FieldValue< 3 >::Scalar > region_id
MultiField< 3, FieldValue< 3 >::Scalar > dg_penalty
Penalty enforcing inter-element continuity of solution (for each substance).
EquationOutput output_fields
Field< 3, FieldValue< 3 >::Scalar > subdomain
MultiField< 3, FieldValue< 3 >::Scalar > fracture_sigma
Transition parameter for diffusive transfer on fractures (for each substance).
Transport with dispersion implemented using discontinuous Galerkin method.
void update_solution() override
Computes the solution in one time instant.
static const Input::Type::Record & get_input_type()
Declare input record type for the equation TransportDG.
bool allocation_done
Indicates whether matrices have been preallocated.
Input::Record input_rec
Array for storing the output solution data.
void zero_time_step() override
Initialize solution in the zero time.
std::shared_ptr< EqFields > eq_fields_
Fields for model parameters.
void initialize() override
GenericAssembly< MassAssemblyDim > * mass_assembly_
general assembly objects, hold assembly objects of appropriate dimension
static const Input::Type::Selection & get_dg_variant_selection_input_type()
Input type for the DG variant selection.
TransportDG(Mesh &init_mesh, const Input::Record in_rec)
Constructor.
~TransportDG() override
Destructor.
void compute_p0_interpolation()
Compute P0 interpolation of the solution (used in reaction term).
void output_data()
Postprocesses the solution and writes to output file.
void calculate_cumulative_balance()
void set_initial_condition()
Calculates the dispersivity (diffusivity) tensor from the velocity field.
std::shared_ptr< EqData > eq_data_
Data for model parameters.
void update_after_reactions(bool solution_changed)
static UnitSI & dimensionless()
Returns dimensionless unit.
Discontinuous Galerkin method for equation of transport with dispersion.
Definitions of basic Lagrangean finite elements with polynomial shape functions.
Definitions of Raviart-Thomas finite elements.
Class FEValues calculates finite element data on the actual cells such as shape function values,...
#define FLOW123D_FORCE_LINK_IN_CHILD(x)
Discontinuous Galerkin method for equation of transport with dispersion.
arma::Col< IntIdx > LocDofVec
static constexpr bool value
Solver based on the original PETSc solver using MPIAIJ matrix and succesive Schur complement construc...
#define WarningOut()
Macro defining 'warning' record of log.
Definitions of particular quadrature rules on simplices.
#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.
Discontinuous Galerkin method for equation of transport with dispersion.