35 return Selection(
"Balance_output_format",
"Format of output file for balance.")
43 return Record(
"Balance",
"Balance of a conservative quantity, boundary fluxes and sources.")
45 .
declare_key(
"add_output_times",
Bool(),
Default(
"true"),
"Add all output times of the balanced equation to the balance output times set. " 46 "Note that this is not the time set of the output stream.")
50 "If true, then balance is calculated at each computational time step, which can slow down the program.")
75 : file_prefix_(file_prefix),
79 allocation_done_(false),
81 output_line_counter_(0),
82 output_yaml_header_(false)
170 for (
auto name : names)
177 unsigned int max_dofs_per_boundary)
260 for (
unsigned int c=0; c<n_quant; ++c)
262 chkerr(MatCreateAIJ(PETSC_COMM_WORLD,
267 (
rank_==0)?n_bulk_regs_per_dof:0,
269 (rank_==0)?0:n_bulk_regs_per_dof,
273 chkerr(MatCreateAIJ(PETSC_COMM_WORLD,
278 (rank_==0)?n_bulk_regs_per_dof:0,
280 (rank_==0)?0:n_bulk_regs_per_dof,
284 chkerr(MatCreateAIJ(PETSC_COMM_WORLD,
289 (rank_==0)?n_bulk_regs_per_dof:0,
291 (rank_==0)?0:n_bulk_regs_per_dof,
295 chkerr(MatCreateAIJ(PETSC_COMM_WORLD,
306 chkerr(VecCreateMPI(PETSC_COMM_WORLD,
311 chkerr(VecCreateMPI(PETSC_COMM_WORLD,
316 chkerr(VecCreateMPI(PETSC_COMM_WORLD,
322 chkerr(MatCreateAIJ(PETSC_COMM_WORLD,
334 for (
unsigned int loc_el=0; loc_el<
be_regions_.size(); ++loc_el)
346 chkerr(VecCreateMPI(PETSC_COMM_WORLD,
354 chkerr(VecCreateMPI(PETSC_COMM_WORLD,
365 std::string default_file_name;
472 unsigned int region_idx,
479 PetscInt reg_array[1] = { (int)region_idx };
492 unsigned int boundary_idx,
499 PetscInt elem_array[1] = { int(
be_offset_+boundary_idx) };
511 unsigned int region_idx,
518 PetscInt reg_array[1] = { (int)region_idx };
531 unsigned int region_idx,
542 unsigned int boundary_idx,
556 unsigned int region_idx,
563 PetscInt reg_array[1] = { (int)region_idx };
595 chkerr(VecCreateMPIWithArray(PETSC_COMM_WORLD,
603 chkerr(VecZeroEntries(bulk_vec));
608 chkerr(VecSum(bulk_vec, &sum_sources));
609 VecDestroy(&bulk_vec);
619 chkerr(VecCreateMPIWithArray(PETSC_COMM_WORLD,
627 chkerr(VecZeroEntries(boundary_vec));
633 chkerr(VecScale(temp, -1));
635 chkerr(VecDestroy(&temp));
638 chkerr(VecSum(boundary_vec, &sum_fluxes));
639 chkerr(VecDestroy(&boundary_vec));
656 chkerr(VecCreateMPIWithArray(PETSC_COMM_WORLD,
664 chkerr(VecZeroEntries(bulk_vec));
669 VecDestroy(&bulk_vec);
681 chkerr(VecCreateMPIWithArray(PETSC_COMM_WORLD,
689 chkerr(VecZeroEntries(bulk_vec));
698 const double *sol_array, *mat_array, *rhs_array;
699 chkerr(VecGetLocalSize(solution, &lsize));
700 chkerr(VecDuplicate(solution, &mat_r));
701 chkerr(VecDuplicate(solution, &rhs_r));
702 chkerr(VecGetArrayRead(solution, &sol_array));
708 VecGetArrayRead(mat_r, &mat_array);
709 VecGetArrayRead(rhs_r, &rhs_array);
713 for (
int i=0; i<lsize; ++i)
715 double f = mat_array[i]*sol_array[i] + rhs_array[i];
720 VecRestoreArrayRead(mat_r, &mat_array);
721 VecRestoreArrayRead(rhs_r, &rhs_array);
723 chkerr(VecRestoreArrayRead(solution, &sol_array));
726 VecDestroy(&bulk_vec);
732 chkerr(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1,
734 PETSC_DECIDE, &(
fluxes_[quantity_idx][0]), &boundary_vec));
737 chkerr(VecZeroEntries(boundary_vec));
743 chkerr(VecScale(temp, -1));
749 const double *flux_array;
751 chkerr(VecGetArrayRead(temp, &flux_array));
752 chkerr(VecGetLocalSize(temp, &lsize));
753 for (
int e=0; e<lsize; ++e)
755 if (flux_array[e] < 0)
758 fluxes_in_[quantity_idx][be_regions_[e]] += flux_array[e];
760 chkerr(VecRestoreArrayRead(temp, &flux_array));
762 VecDestroy(&boundary_vec);
779 const int buf_size = n_quant*2*n_blk_reg + n_quant*2*n_bdr_reg;
780 double sendbuffer[buf_size], recvbuffer[buf_size];
781 for (
unsigned int qi=0; qi<n_quant; qi++)
783 for (
unsigned int ri=0; ri<n_blk_reg; ri++)
785 sendbuffer[qi*2*n_blk_reg + + ri] =
sources_in_[qi][ri];
786 sendbuffer[qi*2*n_blk_reg + n_blk_reg + ri] =
sources_out_[qi][ri];
788 for (
unsigned int ri=0; ri<n_bdr_reg; ri++)
790 sendbuffer[n_quant*2*n_blk_reg + qi*2*n_bdr_reg + + ri] =
fluxes_in_[qi][ri];
791 sendbuffer[n_quant*2*n_blk_reg + qi*2*n_bdr_reg + n_bdr_reg + ri] =
fluxes_out_[qi][ri];
801 for (
unsigned int qi=0; qi<n_quant; qi++)
803 for (
unsigned int ri=0; ri<n_blk_reg; ri++)
805 sources_in_[qi][ri] = recvbuffer[qi*2*n_blk_reg + + ri];
806 sources_out_[qi][ri] = recvbuffer[qi*2*n_blk_reg + n_blk_reg + ri];
808 for (
unsigned int ri=0; ri<n_bdr_reg; ri++)
810 fluxes_in_[qi][ri] = recvbuffer[n_quant*2*n_blk_reg + qi*2*n_bdr_reg + + ri];
811 fluxes_out_[qi][ri] = recvbuffer[n_quant*2*n_blk_reg + qi*2*n_bdr_reg + n_bdr_reg + ri];
831 for( RegionSet::const_iterator reg = b_set.begin(); reg != b_set.end(); ++reg)
833 for (
unsigned int qi=0; qi<n_quant; qi++)
843 for( RegionSet::const_iterator reg = bulk_set.begin(); reg != bulk_set.end(); ++reg)
845 for (
unsigned int qi=0; qi<n_quant; qi++)
861 for (
unsigned int qi=0; qi<n_quant; qi++)
866 for (
unsigned int qi=0; qi<n_quant; qi++)
906 if (
rank_ != 0)
return;
913 unsigned int wl = 2*(w-5)+7;
914 string bc_head_format =
"# %-*s%-*s%-*s%-*s%-*s%-*s\n",
915 bc_format =
"%*s%-*d%-*s%-*s%-*g%-*g%-*g\n",
916 bc_total_format =
"# %-*s%-*s%-*g%-*g%-*g\n";
918 output_ <<
"# " << setw((w*c+wl-14)/2) << setfill(
'-') <<
"--" 920 << setw((w*c+wl-14)/2) << setfill(
'-') <<
"" << endl
921 <<
"# Time: " << time <<
"\n\n\n";
924 output_ <<
"# Mass flux through boundary [M/T]:\n# " 925 << setiosflags(ios::left) << setfill(
' ')
926 << setw(w) <<
"[boundary_id]" 927 << setw(wl) <<
"[label]" 928 << setw(w) <<
"[substance]" 929 << setw(w) <<
"[total flux]" 930 << setw(w) <<
"[outward flux]" 931 << setw(w) <<
"[inward flux]" 935 output_ <<
"# " << setw(w*c+wl) << setfill(
'-') <<
"" << setfill(
' ') << endl;
939 for( RegionSet::const_iterator reg = b_set.begin(); reg != b_set.end(); ++reg) {
940 for (
unsigned int qi=0; qi<n_quant; qi++) {
942 << setw(w) << (int)reg->id()
943 << setw(wl) << reg->label().c_str()
945 << setw(w) <<
fluxes_[qi][reg->boundary_idx()]
947 << setw(w) <<
fluxes_in_[qi][reg->boundary_idx()]
953 output_ <<
"# " << setw(w*c+wl) << setfill(
'-') <<
"" << setfill(
' ') << endl;
956 for (
unsigned int qi=0; qi<n_quant; qi++)
957 output_ <<
"# " << setiosflags(ios::left)
958 << setw(w+wl) <<
"Total mass flux of substance [M/T]" 968 string src_head_format =
"# %-*s%-*s%-*s%-*s%-*s\n",
969 src_format =
"%*s%-*d%-*s%-*s%-*g%-*g\n",
970 src_total_format =
"# %-*s%-*s%-*g%-*g\n";
971 output_ <<
"# Mass [M] and sources [M/T] on regions:\n" 972 <<
"# " << setiosflags(ios::left)
973 << setw(w) <<
"[region_id]" 974 << setw(wl) <<
"[label]" 975 << setw(w) <<
"[substance]" 976 << setw(w) <<
"[total_mass]" 977 << setw(w) <<
"[total_source]" 981 output_ <<
"# " << setw(w*c+wl) << setfill(
'-') <<
"" << setfill(
' ') << endl;
985 for( RegionSet::const_iterator reg = bulk_set.begin(); reg != bulk_set.end(); ++reg)
987 for (
unsigned int qi=0; qi<n_quant; qi++)
990 << setw(w) << (int)reg->id()
991 << setw(wl) << reg->label().c_str()
993 << setw(w) <<
masses_[qi][reg->bulk_idx()]
994 << setw(w) <<
sources_[qi][reg->bulk_idx()]
1000 output_ <<
"# " << setw(w*c+wl) << setfill(
'-') <<
"" << setfill(
' ') << endl;
1003 for (
unsigned int qi=0; qi<n_quant; qi++)
1004 output_ <<
"# " << setiosflags(ios::left) << setw(w+wl) <<
"Total mass [M] and sources [M/T]" 1013 output_ <<
"\n\n# Cumulative mass balance on time interval [" 1015 << setiosflags(ios::left) << time <<
"]\n" 1016 <<
"# Initial mass [M] + sources integrated over time [M] - flux integrated over time [M] = current mass [M]\n" 1017 <<
"# " << setiosflags(ios::left)
1018 << setw(w) <<
"[substance]" 1019 << setw(w) <<
"[A=init. mass]" 1020 << setw(w) <<
"[B=source]" 1021 << setw(w) <<
"[C=flux]" 1022 << setw(w) <<
"[A+B-C]" 1023 << setw(w) <<
"[D=curr. mass]" 1024 << setw(w) <<
"[A+B-C-D=err.]" 1025 << setw(w) <<
"[rel. error]" 1028 for (
unsigned int qi=0; qi<n_quant; qi++)
1031 output_ <<
" " << setiosflags(ios::left)
1050 std::stringstream ss;
1051 for (
unsigned int i=0; i<cnt; i++) ss <<
format_csv_val(0, delimiter);
1059 if (
rank_ != 0)
return;
1068 for( RegionSet::const_iterator reg = bulk_set.begin(); reg != bulk_set.end(); ++reg)
1070 for (
unsigned int qi=0; qi<n_quant; qi++)
1090 for( RegionSet::const_iterator reg = b_set.begin(); reg != b_set.end(); ++reg)
1092 for (
unsigned int qi=0; qi<n_quant; qi++) {
1109 for (
unsigned int qi=0; qi<n_quant; qi++)
1139 std::stringstream ss;
1140 if (delimiter ==
' ') {
1146 output_ << comment_string << ss.str()
1166 std::stringstream ss;
1167 std::replace( val.begin(), val.end(),
'\"',
'\'');
1169 if (!initial) ss << delimiter;
1170 if (delimiter ==
' ') {
1171 std::stringstream sval;
1172 sval <<
"\"" << val <<
"\"";
1175 ss <<
"\"" << val <<
"\"";
1183 std::stringstream ss;
1185 if (!initial) ss << delimiter;
1186 if (delimiter ==
' ') {
1197 if (
rank_ != 0)
return;
1203 output_yaml_ <<
"column_names: [ flux, flux_in, flux_out, mass, source, source_in, source_out, flux_increment, " 1204 <<
"source_increment, flux_cumulative, source_cumulative, error ]" << endl;
1213 for( RegionSet::const_iterator reg = bulk_set.begin(); reg != bulk_set.end(); ++reg)
1215 for (
unsigned int qi=0; qi<n_quant; qi++)
1218 output_yaml_ << setw(4) <<
"" <<
"region: " << reg->label() << endl;
1220 output_yaml_ << setw(4) <<
"" <<
"data: " <<
"[ 0, 0, 0, " <<
masses_[qi][reg->bulk_idx()] <<
", " 1222 <<
sources_out_[qi][reg->bulk_idx()] <<
", 0, 0, 0, 0, 0 ]" << endl;
1228 for( RegionSet::const_iterator reg = b_set.begin(); reg != b_set.end(); ++reg)
1230 for (
unsigned int qi=0; qi<n_quant; qi++) {
1232 output_yaml_ << setw(4) <<
"" <<
"region: " << reg->label() << endl;
1234 output_yaml_ << setw(4) <<
"" <<
"data: " <<
"[ " <<
fluxes_[qi][reg->boundary_idx()] <<
", " 1236 <<
", 0, 0, 0, 0, 0, 0, 0, 0, 0 ]" << endl;
1242 for (
unsigned int qi=0; qi<n_quant; qi++)
1246 output_yaml_ << setw(4) <<
"" <<
"region: ALL" << endl;
1254 << error <<
" ]" << endl;
UnitSI units_
Units of conserved quantities.
static const Input::Type::Record & get_input_type()
Main balance input record type.
std::vector< double > integrated_sources_
unsigned int add_quantity(const string &name)
std::vector< double > sum_fluxes_out_
void add_mass_matrix_values(unsigned int quantity_idx, unsigned int region_idx, const std::vector< IdxInt > &dof_indices, const std::vector< double > &values)
void calculate_mass(unsigned int quantity_idx, const Vec &solution, vector< double > &output_array)
TimeMark::Type output_mark_type_
TimeMark type for output of particular equation.
void allocate(unsigned int n_loc_dofs, unsigned int max_dofs_per_boundary)
TimeMark::Type balance_output_type_
TimeMark type for balance output of particular equation.
unsigned int * boundary_idx_
std::string format_text() const
void finish_source_assembly(unsigned int quantity_idx)
This method must be called after assembling the matrix and vectors for computing source.
unsigned int bulk_size() const
void add_mass_vec_value(unsigned int quantity_idx, unsigned int region_idx, double value)
bool initial_
true before calculating the mass at initial time, otherwise false
IdxInt * get_el_4_loc() const
std::vector< std::vector< double > > sources_in_
void add_cumulative_source(unsigned int quantity_idx, double source)
std::vector< std::vector< double > > fluxes_in_
std::vector< double > sum_fluxes_in_
std::vector< std::vector< double > > fluxes_out_
const TimeGovernor * time_
void add_source_vec_values(unsigned int quantity_idx, unsigned int region_idx, const std::vector< IdxInt > &dof_values, const std::vector< double > &values)
Vec ones_
auxiliary vectors for summation of matrix columns
Mat * region_source_rhs_
Matrices for calculation of signed source (n_dofs x n_bulk_regions).
RegionSet get_region_set(const string &set_name) const
Input::Record input_record_
Record for current balance.
ofstream output_
Handle for file for output in given OutputFormat of balance and total fluxes over individual regions ...
void read_from_input(Input::Array in_array, const TimeGovernor &tg)
#define FOR_ELEMENT_SIDES(i, j)
void chkerr(unsigned int ierr)
Replacement of new/delete operator in the spirit of xmalloc.
void output_legacy(double time)
Perform output in old format (for compatibility)
const RegionDB & region_db() const
#define ASSERT(expr)
Allow use shorter versions of macro names if these names is not used with external library...
const TimeStep & step(int index=-1) const
bool allocation_done_
true before allocating necessary internal structures (Petsc matrices etc.)
void calculate_cumulative(unsigned int quantity_idx, const Vec &solution)
Basic time management functionality for unsteady (and steady) solvers (class Equation).
static TimeMarks & marks()
#define MPI_Reduce(sendbuf, recvbuf, count, datatype, op, root, comm)
void calculate_instant(unsigned int quantity_idx, const Vec &solution)
void add_flux_vec_value(unsigned int quantity_idx, unsigned int boundary_idx, double value)
Basic time management class.
std::vector< double > sum_fluxes_
bool balance_on_
If the balance is on. Balance is off in the case of no balance output time marks. ...
static constexpr bool value
std::vector< double > increment_sources_
Vec * be_flux_vec_
Vectors for calculation of flux (n_boundary_edges).
static const Input::Type::Selection & get_format_selection_input_type()
Input selection for file format.
std::vector< double > integrated_fluxes_
unsigned int boundary_idx() const
Returns index of the region in the boundary set.
std::vector< Quantity > quantities_
Names of conserved quantities.
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...
std::vector< std::vector< double > > sources_out_
Mat * be_flux_matrix_
Matrices for calculation of flux (n_boundary_edges x n_dofs).
double last_time_
time of last calculated balance
void open_stream(Stream &stream) const
std::vector< std::vector< double > > fluxes_
void finish_mass_assembly(unsigned int quantity_idx)
This method must be called after assembling the matrix for computing mass.
TimeMark::Type equation_fixed_mark_type() const
void format_csv_output_header(char delimiter, const std::string &comment_string)
Print output header.
TimeMark::Type equation_mark_type() const
void init_from_input(const Input::Record &in_rec, TimeGovernor &tg)
Mat * region_mass_matrix_
Matrices for calculation of mass (n_dofs x n_bulk_regions).
std::vector< std::vector< double > > sources_
std::vector< double > sum_sources_out_
unsigned int max_dofs_per_boundary_
bool cumulative_
if true then cumulative balance is computed
void add_source_matrix_values(unsigned int quantity_idx, unsigned int region_idx, const std::vector< IdxInt > &dof_indices, const std::vector< double > &values)
Mat * region_source_matrix_
Matrices for calculation of source (n_dofs x n_bulk_regions).
std::vector< double > increment_fluxes_
void chkerr_assert(unsigned int ierr)
SideIter side(const unsigned int loc_index)
bool is_current()
Returns true if the current time step is marked for the balance output.
Balance(const std::string &file_prefix, const Mesh *mesh)
void output_csv(double time, char delimiter, const std::string &comment_string, unsigned int repeat=0)
Perform output in csv format.
Distribution * get_el_ds() const
std::vector< unsigned int > add_quantities(const std::vector< string > &names)
Dedicated class for storing path to input and output files.
void output_yaml(double time)
Perform output in yaml format.
void output()
Perform output to file for given time instant.
std::vector< double > sum_sources_
void start_flux_assembly(unsigned int quantity_idx)
std::vector< double > initial_mass_
static const unsigned int output_column_width
Size of column in output (used if delimiter is space)
#define ASSERT_PTR(ptr)
Definition of assert macro checking non-null pointer (PTR)
unsigned int boundary_size() const
void finish_flux_assembly(unsigned int quantity_idx)
This method must be called after assembling the matrix and vector for computing flux.
std::vector< unsigned int > be_regions_
Number of boundary region for each local boundary edge.
OutputFormat output_format_
Format of output file.
static const Input::Type::Array get_input_type()
void add_flux_matrix_values(unsigned int quantity_idx, unsigned int boundary_idx, const std::vector< IdxInt > &dof_indices, const std::vector< double > &values)
void units(const UnitSI &unit)
Setter for units of conserved quantities.
int be_offset_
Offset for local part of vector of boundary edges.
std::vector< std::vector< double > > masses_
Vec * region_source_vec_
Vectors for calculation of source (n_bulk_regions).
void start_mass_assembly(unsigned int quantity_idx)
std::vector< double > sum_masses_
Vec * region_mass_vec_
Vectors for calculation of mass (n_bulk_regions).
unsigned int output_line_counter_
hold count of line printed into output_
bool add_output_times_
Add output time marks to balance output time marks.
Class for representation SI units of Fields.
unsigned int n_loc_dofs_
Allocation parameters. Set by the allocate method used in the lazy_initialize.
void undef(bool val=true)
void start_source_assembly(unsigned int quantity_idx)
FilePath balance_output_file_
File path for output_ stream.
std::string file_prefix_
Save prefix passed in in constructor.
bool output_yaml_header_
marks whether YAML output has printed header
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, align text to column.
ElementVector element
Vector of elements of the mesh.
std::vector< double > sum_sources_in_
unsigned int lsize(int proc) const
get local size