Flow123d
JS_before_hm-1804-gf2ad740aa
|
Go to the documentation of this file.
21 #include <unordered_map>
43 return Selection(
"Balance_output_format",
"Format of output file for balance.")
51 return Record(
"Balance",
"Balance of a conservative quantity, boundary fluxes and sources.")
53 .
declare_key(
"add_output_times",
Bool(),
Default(
"true"),
"Add all output times of the balanced equation to the balance output times set. "
54 "Note that this is not the time set of the output stream.")
58 "If true, then balance is calculated at each computational time step, which can slow down the program.")
87 : file_prefix_(file_prefix),
91 allocation_done_(false),
93 output_line_counter_(0),
94 output_yaml_header_(false)
176 for (
auto name : names)
183 unsigned int max_dofs_per_boundary)
192 unsigned int max_dofs_per_boundary)
232 unsigned int be_id = 0;
233 for (
unsigned int loc_el = 0; loc_el <
mesh_->
get_el_ds()->lsize(); loc_el++)
238 for(
unsigned int si=0; si<elm->
n_sides(); si++)
287 for (
unsigned int c=0; c<n_quant; ++c)
289 chkerr(MatCreateAIJ(PETSC_COMM_WORLD,
294 (
rank_==0)?n_bulk_regs_per_dof:0,
296 (
rank_==0)?0:n_bulk_regs_per_dof,
300 chkerr(MatCreateSeqAIJ(PETSC_COMM_SELF,
307 chkerr(MatCreateSeqAIJ(PETSC_COMM_SELF,
314 chkerr(MatCreateAIJ(PETSC_COMM_WORLD,
325 chkerr(VecCreateMPI(PETSC_COMM_WORLD,
330 chkerr(VecCreateMPI(PETSC_COMM_WORLD,
341 std::string default_file_name;
455 uint m = mat_values.size();
457 for (
uint i=0; i<m; i++)
488 uint m = mat_values.size();
490 for (
uint i=0; i<m; i++)
511 unsigned int region_idx,
519 PetscInt reg_array[1] = { (int)region_idx };
522 loc_dof_indices.size(),
523 loc_dof_indices.memptr(),
526 &(mult_mat_values[0]),
530 loc_dof_indices.size(),
531 loc_dof_indices.memptr(),
534 &(add_mat_values[0]),
557 double temp_source = 0;
558 int lsize, n_cols_mat, n_cols_rhs;
560 const double *vals_mat, *vals_rhs, *sol_array;
564 chkerr(VecGetArrayRead(solution, &sol_array));
570 for (
int i=0; i<lsize; ++i){
576 for (
int j=0; j<n_cols_mat; ++j)
577 temp_source += vals_mat[j]*sol_array[i] + vals_rhs[j];
582 chkerr(VecRestoreArrayRead(solution, &sol_array));
589 chkerr(VecCreateMPI(PETSC_COMM_WORLD,
597 chkerr(VecSum(temp, &sum_fluxes));
598 chkerr(VecDestroy(&temp));
617 chkerr(VecCreateMPIWithArray(PETSC_COMM_WORLD,
625 chkerr(VecZeroEntries(bulk_vec));
630 chkerr(VecDestroy(&bulk_vec));
646 int lsize, n_cols_mat, n_cols_rhs;
648 const double *vals_mat, *vals_rhs, *sol_array;
652 chkerr(VecGetArrayRead(solution, &sol_array));
657 for (
int i=0; i<lsize; ++i)
665 for (
int j=0; j<n_cols_mat; ++j)
669 double f = vals_mat[j]*sol_array[i] + vals_rhs[j];
674 chkerr(VecRestoreArrayRead(solution, &sol_array));
678 chkerr(VecCreateMPI(PETSC_COMM_WORLD,
686 chkerr(VecScale(temp, -1));
691 const double *flux_array;
693 chkerr(VecGetArrayRead(temp, &flux_array));
694 chkerr(VecGetLocalSize(temp, &lsize));
695 for (
int e=0; e<lsize; ++e)
697 if (flux_array[e] < 0)
702 chkerr(VecRestoreArrayRead(temp, &flux_array));
703 chkerr(VecDestroy(&temp));
720 const int buf_size = n_quant*2*n_blk_reg + n_quant*2*n_bdr_reg + n_quant;
721 double sendbuffer[buf_size], recvbuffer[buf_size];
722 for (
unsigned int qi=0; qi<n_quant; qi++)
724 for (
unsigned int ri=0; ri<n_blk_reg; ri++)
726 sendbuffer[qi*2*n_blk_reg + + ri] =
sources_in_[qi][ri];
727 sendbuffer[qi*2*n_blk_reg + n_blk_reg + ri] =
sources_out_[qi][ri];
729 for (
unsigned int ri=0; ri<n_bdr_reg; ri++)
731 sendbuffer[n_quant*2*n_blk_reg + qi*2*n_bdr_reg + + ri] =
fluxes_in_[qi][ri];
732 sendbuffer[n_quant*2*n_blk_reg + qi*2*n_bdr_reg + n_bdr_reg + ri] =
fluxes_out_[qi][ri];
747 for (
unsigned int qi=0; qi<n_quant; qi++)
749 for (
unsigned int ri=0; ri<n_blk_reg; ri++)
751 sources_in_[qi][ri] = recvbuffer[qi*2*n_blk_reg + + ri];
752 sources_out_[qi][ri] = recvbuffer[qi*2*n_blk_reg + n_blk_reg + ri];
754 for (
unsigned int ri=0; ri<n_bdr_reg; ri++)
756 fluxes_in_[qi][ri] = recvbuffer[n_quant*2*n_blk_reg + qi*2*n_bdr_reg + + ri];
757 fluxes_out_[qi][ri] = recvbuffer[n_quant*2*n_blk_reg + qi*2*n_bdr_reg + n_bdr_reg + ri];
781 for( RegionSet::const_iterator reg = b_set.begin(); reg != b_set.end(); ++reg)
783 for (
unsigned int qi=0; qi<n_quant; qi++)
793 for( RegionSet::const_iterator reg = bulk_set.begin(); reg != bulk_set.end(); ++reg)
795 for (
unsigned int qi=0; qi<n_quant; qi++)
811 for (
unsigned int qi=0; qi<n_quant; qi++)
816 for (
unsigned int qi=0; qi<n_quant; qi++)
856 if (
rank_ != 0)
return;
863 unsigned int wl = 2*(w-5)+7;
864 string bc_head_format =
"# %-*s%-*s%-*s%-*s%-*s%-*s\n",
865 bc_format =
"%*s%-*d%-*s%-*s%-*g%-*g%-*g\n",
866 bc_total_format =
"# %-*s%-*s%-*g%-*g%-*g\n";
868 output_ <<
"# " << setw((w*c+wl-14)/2) << setfill(
'-') <<
"--"
870 << setw((w*c+wl-14)/2) << setfill(
'-') <<
"" << endl
875 output_ <<
"# Mass flux through boundary [M/T]:\n# "
876 << setiosflags(ios::left) << setfill(
' ')
877 << setw(w) <<
"[boundary_id]"
878 << setw(wl) <<
"[label]"
879 << setw(w) <<
"[substance]"
880 << setw(w) <<
"[total flux]"
881 << setw(w) <<
"[outward flux]"
882 << setw(w) <<
"[inward flux]"
886 output_ <<
"# " << setw(w*c+wl) << setfill(
'-') <<
"" << setfill(
' ') << endl;
890 for( RegionSet::const_iterator reg = b_set.begin(); reg != b_set.end(); ++reg) {
891 for (
unsigned int qi=0; qi<n_quant; qi++) {
893 << setw(w) << (int)reg->id()
894 << setw(wl) << reg->label().c_str()
898 << setw(w) <<
fluxes_in_[qi][reg->boundary_idx()]
904 output_ <<
"# " << setw(w*c+wl) << setfill(
'-') <<
"" << setfill(
' ') << endl;
907 for (
unsigned int qi=0; qi<n_quant; qi++)
908 output_ <<
"# " << setiosflags(ios::left)
909 << setw(w+wl) <<
"Total mass flux of substance [M/T]"
919 string src_head_format =
"# %-*s%-*s%-*s%-*s%-*s\n",
920 src_format =
"%*s%-*d%-*s%-*s%-*g%-*g\n",
921 src_total_format =
"# %-*s%-*s%-*g%-*g\n";
922 output_ <<
"# Mass [M] and sources [M/T] on regions:\n"
923 <<
"# " << setiosflags(ios::left)
924 << setw(w) <<
"[region_id]"
925 << setw(wl) <<
"[label]"
926 << setw(w) <<
"[substance]"
927 << setw(w) <<
"[total_mass]"
928 << setw(w) <<
"[total_source]"
932 output_ <<
"# " << setw(w*c+wl) << setfill(
'-') <<
"" << setfill(
' ') << endl;
936 for( RegionSet::const_iterator reg = bulk_set.begin(); reg != bulk_set.end(); ++reg)
938 for (
unsigned int qi=0; qi<n_quant; qi++)
941 << setw(w) << (int)reg->id()
942 << setw(wl) << reg->label().c_str()
944 << setw(w) <<
masses_[qi][reg->bulk_idx()]
951 output_ <<
"# " << setw(w*c+wl) << setfill(
'-') <<
"" << setfill(
' ') << endl;
954 for (
unsigned int qi=0; qi<n_quant; qi++)
955 output_ <<
"# " << setiosflags(ios::left) << setw(w+wl) <<
"Total mass [M] and sources [M/T]"
964 output_ <<
"\n\n# Cumulative mass balance on time interval ["
966 << setiosflags(ios::left) << time <<
"]\n"
967 <<
"# Initial mass [M] + sources integrated over time [M] - flux integrated over time [M] = current mass [M]\n"
968 <<
"# " << setiosflags(ios::left)
969 << setw(w) <<
"[substance]"
970 << setw(w) <<
"[A=init. mass]"
971 << setw(w) <<
"[B=source]"
972 << setw(w) <<
"[C=flux]"
973 << setw(w) <<
"[A+B-C]"
974 << setw(w) <<
"[D=curr. mass]"
975 << setw(w) <<
"[A+B-C-D=err.]"
976 << setw(w) <<
"[rel. error]"
979 for (
unsigned int qi=0; qi<n_quant; qi++)
982 output_ <<
" " << setiosflags(ios::left)
1001 std::stringstream ss;
1002 for (
unsigned int i=0; i<cnt; i++) ss <<
format_csv_val(0, delimiter);
1010 if (
rank_ != 0)
return;
1019 for( RegionSet::const_iterator reg = bulk_set.begin(); reg != bulk_set.end(); ++reg)
1021 for (
unsigned int qi=0; qi<n_quant; qi++)
1041 for( RegionSet::const_iterator reg = b_set.begin(); reg != b_set.end(); ++reg)
1043 for (
unsigned int qi=0; qi<n_quant; qi++) {
1060 for (
unsigned int qi=0; qi<n_quant; qi++)
1090 std::stringstream ss;
1091 if (delimiter ==
' ') {
1098 output_ << comment_string << ss.str()
1118 std::stringstream ss;
1119 std::replace( val.begin(), val.end(),
'\"',
'\'');
1121 if (!initial) ss << delimiter;
1122 if (delimiter ==
' ') {
1123 std::stringstream sval;
1124 sval <<
"\"" << val <<
"\"";
1127 ss <<
"\"" << val <<
"\"";
1135 std::stringstream ss;
1137 if (!initial) ss << delimiter;
1138 if (delimiter ==
' ') {
1156 output_yaml_ <<
"column_names: [ flux, flux_in, flux_out, mass, source, source_in, source_out, flux_increment, "
1157 <<
"source_increment, flux_cumulative, source_cumulative, error ]" << endl;
1166 for( RegionSet::const_iterator reg = bulk_set.begin(); reg != bulk_set.end(); ++reg)
1168 for (
unsigned int qi=0; qi<n_quant; qi++)
1171 output_yaml_ << setw(4) <<
"" <<
"region: " << reg->label() << endl;
1173 output_yaml_ << setw(4) <<
"" <<
"data: " <<
"[ 0, 0, 0, " <<
masses_[qi][reg->bulk_idx()] <<
", "
1176 <<
", 0, 0, 0, 0, 0 ]" << endl;
1182 for( RegionSet::const_iterator reg = b_set.begin(); reg != b_set.end(); ++reg)
1184 for (
unsigned int qi=0; qi<n_quant; qi++) {
1186 output_yaml_ << setw(4) <<
"" <<
"region: " << reg->label() << endl;
1191 <<
", 0, 0, 0, 0, 0, 0, 0, 0, 0 ]" << endl;
1197 for (
unsigned int qi=0; qi<n_quant; qi++)
1201 output_yaml_ << setw(4) <<
"" <<
"region: ALL" << endl;
1209 << error <<
" ]" << endl;
arma::Col< IntIdx > LocDofVec
Distribution * get_el_ds() const
void calculate_instant(unsigned int quantity_idx, const Vec &solution)
void allocate(const std::shared_ptr< DOFHandlerMultiDim > &dh, unsigned int max_dofs_per_boundary)
static const unsigned int output_column_width
Size of column in output (used if delimiter is space)
std::vector< unsigned int > be_regions_
Maps local boundary edge to its region boundary index.
Side side() const
Return Side of given cell and side_idx.
const DOFHandlerMultiDim * dh() const
Return DOF handler.
unsigned int n_loc_dofs_seq_
std::vector< double > sum_fluxes_out_
Balance(const std::string &file_prefix, const Mesh *mesh)
Basic time management class.
UnitSI units_
Units of conserved quantities.
const std::vector< LongIdx > & get_local_to_global_map() const
Get the map between local dof indices and the global ones.
std::string format_text() const
const TimeGovernor * time_
std::shared_ptr< TimeUnitConversion > get_unit_conversion() const
Getter for time unit conversion object.
#define ASSERT(expr)
Allow use shorter versions of macro names if these names is not used with external library.
std::string csv_zero_vals(unsigned int cnt, char delimiter)
Return part of output represented by zero values. Count of zero values is given by cnt parameter.
const unsigned int index_
Internal index within list of quantities.
unsigned int add_quantity(const string &name)
TimeMark::Type balance_output_type_
TimeMark type for balance output of particular equation.
TimeMark::Type equation_fixed_mark_type() const
Support classes for parallel programing.
std::vector< double > integrated_sources_
bool is_boundary() const
Returns true for side on the boundary.
void format_csv_output_header(char delimiter, const std::string &comment_string)
Print output header.
std::vector< double > sum_sources_
void finish_source_assembly(unsigned int quantity_idx)
This method must be called after assembling the matrix and vectors for computing source.
Dedicated class for storing path to input and output files.
void chkerr(unsigned int ierr)
Replacement of new/delete operator in the spirit of xmalloc.
bool output_yaml_header_
marks whether YAML output has printed header
void start_flux_assembly(unsigned int quantity_idx)
std::vector< double > increment_fluxes_
void init_from_input(const Input::Record &in_rec, TimeGovernor &tg)
static bool do_yaml_output_
Declaration of class which handles the ordering of degrees of freedom (dof) and mappings between loca...
std::vector< Quantity > quantities_
Names of conserved quantities.
void open_stream(Stream &stream) const
void units(const UnitSI &unit)
Setter for units of conserved quantities.
void add_flux_values(unsigned int quantity_idx, const DHCellSide &side, const LocDofVec &loc_dof_indices, const std::vector< double > &mat_values, double vec_value)
unsigned int boundary_size() const
void start_mass_assembly(unsigned int quantity_idx)
Vec * region_mass_vec_
Vectors for calculation of mass (n_bulk_regions).
Mat * be_flux_matrix_
Matrices for calculation of flux (n_boundary_edges x n_dofs).
void finish_mass_assembly(unsigned int quantity_idx)
This method must be called after assembling the matrix for computing mass.
Side accessor allows to iterate over sides of DOF handler cell.
Class used for marking specified times at which some events occur.
std::ofstream output_
Handle for file for output in given OutputFormat of balance and total fluxes over individual regions ...
void output_yaml(double time)
Perform output in yaml format.
const DHCellAccessor & cell() const
Return DHCellAccessor appropriate to the side.
Mat * region_mass_matrix_
Matrices for calculation of mass (n_dofs x n_bulk_regions).
bool balance_on_
If the balance is on. Balance is off in the case of no balance output time marks.
unsigned int bulk_size() const
bool cumulative_
if true then cumulative balance is computed
std::vector< std::vector< double > > masses_
Basic time management functionality for unsteady (and steady) solvers (class Equation).
const TimeStep & step(int index=-1) const
Vec * be_flux_vec_
Vectors for calculation of flux (n_boundary_edges).
const ElementAccessor< 3 > elm() const
Return ElementAccessor to element of loc_ele_idx_.
void chkerr_assert(unsigned int ierr)
std::vector< double > sum_fluxes_in_
Class for representation SI units of Fields.
const RegionDB & region_db() const
std::vector< double > sum_masses_
std::vector< std::vector< double > > fluxes_
RegionSet get_region_set(const std::string &set_name) const
#define MPI_Reduce(sendbuf, recvbuf, count, datatype, op, root, comm)
unsigned int output_line_counter_
hold count of line printed into output_
void output()
Perform output to file for given time instant.
std::vector< double > initial_mass_
bool add_output_times_
Add output time marks to balance output time marks.
void calculate_mass(unsigned int quantity_idx, const Vec &solution, vector< double > &output_array)
void read_from_input(Input::Array in_array, const TimeGovernor &tg)
Mat * region_source_matrix_
Matrices for calculation of source (n_dofs x n_bulk_regions).
bool is_current()
Returns true if the current time step is marked for the balance output.
std::string format_csv_val(std::string val, char delimiter, bool initial=false)
Format string value of csv output. Wrap string into quotes and if delimiter is space,...
int LongIdx
Define type that represents indices of large arrays (elements, nodes, dofs etc.)
void calculate_cumulative(unsigned int quantity_idx, const Vec &solution)
OutputFormat output_format_
Format of output file.
std::vector< std::vector< double > > fluxes_in_
int be_offset_
Offset for local part of vector of boundary edges.
std::vector< std::vector< double > > sources_in_
void add_source_values(unsigned int quantity_idx, unsigned int region_idx, const LocDofVec &loc_dof_indices, const std::vector< double > &mult_mat_values, const std::vector< double > &add_mat_values)
void undef(bool val=true)
std::vector< double > sum_fluxes_
std::vector< double > increment_sources_
std::string file_prefix_
Save prefix passed in in constructor.
unsigned int n_sides() const
double last_time_
time of last calculated balance
static const Input::Type::Array get_input_type()
static TimeMarks & marks()
static void set_yaml_output()
Set global variable to output balance files into YAML format (in addition to the table format).
Cell accessor allow iterate over DOF handler cells.
void output_csv(double time, char delimiter, const std::string &comment_string, unsigned int repeat=0)
Perform output in csv format.
std::vector< std::vector< double > > fluxes_out_
static const Input::Type::Selection & get_format_selection_input_type()
Input selection for file format.
LongIdx get_boundary_edge_uid(SideIter side)
unsigned int n_loc_dofs_par_
Allocation parameters. Set by the allocate method used in the lazy_initialize.
unsigned int boundary_idx() const
Returns index of the region in the boundary set.
Input::Record input_record_
Record for current balance.
TimeMark::Type output_mark_type_
TimeMark type for output of particular equation.
TimeMark::Type equation_mark_type() const
void start_source_assembly(unsigned int quantity_idx)
std::vector< double > sum_sources_in_
bool allocation_done_
true before allocating necessary internal structures (Petsc matrices etc.)
std::vector< double > sum_sources_out_
FilePath balance_output_file_
File path for output_ stream.
std::ofstream output_yaml_
bool initial_
true before calculating the mass at initial time, otherwise false
unsigned int max_dofs_per_boundary_
void add_mass_values(unsigned int quantity_idx, const DHCellAccessor &dh_cell, const LocDofVec &loc_dof_indices, const std::vector< double > &mat_values, double vec_value)
void add_cumulative_source(unsigned int quantity_idx, double source)
#define ASSERT_PTR(ptr)
Definition of assert macro checking non-null pointer (PTR)
SideIter side(const unsigned int loc_index)
void finish_flux_assembly(unsigned int quantity_idx)
This method must be called after assembling the matrix and vector for computing flux.
virtual ElementAccessor< 3 > element_accessor(unsigned int idx) const
Create and return ElementAccessor to element of given idx.
std::vector< double > integrated_fluxes_
unsigned int bulk_idx() const
Returns index of the region in the bulk set.
Mat * region_source_rhs_
Matrices for calculation of signed source (n_dofs x n_bulk_regions).
std::vector< std::vector< double > > sources_out_
unsigned int * boundary_idx_
std::unordered_map< LongIdx, unsigned int > be_id_map_
static const Input::Type::Record & get_input_type()
Main balance input record type.
void output_legacy(double time)
Perform output in old format (for compatibility)
RegionIdx region_idx() const
std::vector< unsigned int > add_quantities(const std::vector< string > &names)
LongIdx * get_el_4_loc() const