Flow123d  release_2.2.0-914-gf1a3a4f
balance.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 balance.cc
15  * @ingroup transport
16  * @brief Mass balance
17  */
18 
19 #include <iostream>
20 #include <iomanip>
21 
22 #include "system/system.hh"
23 #include "system/sys_profiler.hh"
24 
25 #include <petscmat.h>
26 #include "mesh/mesh.h"
27 #include "io/output_time_set.hh"
28 #include "coupling/balance.hh"
29 #include "fields/unit_si.hh"
30 #include "tools/time_governor.hh"
31 
32 using namespace Input::Type;
33 
35  return Selection("Balance_output_format", "Format of output file for balance.")
36  .add_value(Balance::legacy, "legacy", "Legacy format used by previous program versions.")
37  .add_value(Balance::txt, "txt", "Excel format with tab delimiter.")
38  .add_value(Balance::gnuplot, "gnuplot", "Format compatible with GnuPlot datafile with fixed column width.")
39  .close();
40 }
41 
43  return Record("Balance", "Balance of a conservative quantity, boundary fluxes and sources.")
44  .declare_key("times", OutputTimeSet::get_input_type(), Default("[]"), "" )
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.")
47  //.declare_key("balance_on", Bool(), Default("true"), "Balance is computed if the value is true.")
48  .declare_key("format", Balance::get_format_selection_input_type(), Default("\"txt\""), "Format of output file.")
49  .declare_key("cumulative", Bool(), Default("false"), "Compute cumulative balance over time. "
50  "If true, then balance is calculated at each computational time step, which can slow down the program.")
51  .declare_key("file", FileName::output(), Default::read_time("File name generated from the balanced quantity: <quantity_name>_balance.*"), "File name for output of balance.")
52  .close();
53 }
54 
55 /*
56 std::shared_ptr<Balance> Balance::make_balance(
57  const std::string &file_prefix,
58  const Mesh *mesh,
59  const Input::Record &in_rec,
60  TimeGovernor &tg)
61 {
62  auto ptr = make_shared<Balance>(file_prefix, mesh, in_rec, tg);
63  auto &marks = TimeGovernor::marks();
64  auto balance_output_type = tg.equation_mark_type() | TimeGovernor::marks().type_balance_output();
65  if (marks.begin(balance_output_type) == marks.end(balance_output_type) ) {
66  // no balance output time => force balance off
67  ptr.reset();
68  }
69  return ptr;
70 }*/
71 
72 
73 
74 Balance::Balance(const std::string &file_prefix, const Mesh *mesh)
75  : file_prefix_(file_prefix),
76  mesh_(mesh),
77  last_time_(),
78  initial_(true),
79  allocation_done_(false),
80  balance_on_(true),
81  output_line_counter_(0),
82  output_yaml_header_(false)
83 
84 {
85  MPI_Comm_rank(PETSC_COMM_WORLD, &rank_);
87 }
88 
89 
90 
92 {
93  if (rank_ == 0) {
94  output_.close();
95  output_yaml_.close();
96  }
97  if (! allocation_done_) return;
98 
99  for (unsigned int c=0; c<quantities_.size(); ++c)
100  {
101  chkerr(MatDestroy(&(region_mass_matrix_[c])));
102  chkerr(MatDestroy(&(be_flux_matrix_[c])));
103  chkerr(MatDestroy(&(region_source_matrix_[c])));
104  chkerr(MatDestroy(&(region_source_rhs_[c])));
105  chkerr(VecDestroy(&(be_flux_vec_[c])));
106  chkerr(VecDestroy(&(region_source_vec_[c])));
107  chkerr(VecDestroy(&(region_mass_vec_[c])));
108  }
109  delete[] region_mass_matrix_;
110  delete[] be_flux_matrix_;
111  delete[] be_flux_vec_;
112  delete[] region_source_matrix_;
113  delete[] region_source_rhs_;
114  delete[] region_source_vec_;
115  delete[] region_mass_vec_;
116 
117  chkerr(MatDestroy(&region_be_matrix_));
118  chkerr(VecDestroy(&ones_));
119  chkerr(VecDestroy(&ones_be_));
120 }
121 
122 
124  const Input::Record &in_rec,
125  TimeGovernor &tg)
126 {
128 
129  time_ = &tg;
130 
131  auto &marks = TimeGovernor::marks();
132 
133  output_mark_type_ = tg.equation_mark_type() | marks.type_output(),
134  balance_output_type_ = tg.equation_fixed_mark_type() | marks.type_balance_output();
135 
136  cumulative_ = in_rec.val<bool>("cumulative");
137  output_format_ = in_rec.val<OutputFormat>("format");
138 
139  OutputTimeSet time_set;
140  time_set.read_from_input( in_rec.val<Input::Array>("times"), tg, balance_output_type_);
141  add_output_times_ = in_rec.val<bool>("add_output_times");
142 
143  input_record_ = in_rec;
144 
145 }
146 
147 
148 void Balance::units(const UnitSI &unit)
149 {
151 
152  units_ = unit;
153  units_.undef(false);
154 }
155 
156 unsigned int Balance::add_quantity(const string &name)
157 {
159 
160  Quantity q(quantities_.size(), name);
161  quantities_.push_back(q);
162 
163  return q.index_;
164 }
165 
166 
168 {
169  vector<unsigned int> indices;
170  for (auto name : names)
171  indices.push_back(add_quantity(name));
172  return indices;
173 }
174 
175 
176 void Balance::allocate(unsigned int n_loc_dofs,
177  unsigned int max_dofs_per_boundary)
178 {
180  n_loc_dofs_ = n_loc_dofs;
181  max_dofs_per_boundary_ = max_dofs_per_boundary;
182 }
183 
185 {
186  if (allocation_done_) return;
187 
188  auto &marks = TimeGovernor::marks();
189  if (add_output_times_)
190  marks.add_to_type_all(output_mark_type_, balance_output_type_);
191  // if there are no balance marks turn balance off
192  if (marks.begin(balance_output_type_) == marks.end(balance_output_type_) )
193  {
194  balance_on_ = false;
195  cumulative_ = false;
196  return;
197  }
198 
199  // Max. number of regions to which a single dof can contribute.
200  // TODO: estimate or compute this number directly (from mesh or dof handler).
201  const int n_bulk_regs_per_dof = min(10, (int)mesh_->region_db().bulk_size());
202  const unsigned int n_quant = quantities_.size();
203  const unsigned int n_bdr_reg = mesh_->region_db().boundary_size();
204  const unsigned int n_blk_reg = mesh_->region_db().bulk_size();
205 
206 
207  // construct vector of regions of boundary edges
208  for (unsigned int loc_el = 0; loc_el < mesh_->get_el_ds()->lsize(); loc_el++)
209  {
210  Element *elm = mesh_->element(mesh_->get_el_4_loc()[loc_el]);
211  if (elm->boundary_idx_ != nullptr)
212  {
213  FOR_ELEMENT_SIDES(elm,si)
214  {
215  Boundary *b = elm->side(si)->cond();
216  if (b != nullptr)
217  be_regions_.push_back(b->region().boundary_idx());
218  }
219  }
220  }
221 
222 
223 
224  fluxes_ .resize(n_quant, vector<double>(n_bdr_reg, 0));
225  fluxes_in_ .resize(n_quant, vector<double>(n_bdr_reg, 0));
226  fluxes_out_.resize(n_quant, vector<double>(n_bdr_reg, 0));
227 
228  masses_ .resize(n_quant, vector<double>(n_blk_reg, 0));
229  sources_ .resize(n_quant, vector<double>(n_blk_reg, 0));
230  sources_in_ .resize(n_quant, vector<double>(n_blk_reg, 0));
231  sources_out_.resize(n_quant, vector<double>(n_blk_reg, 0));
232 
233  sum_fluxes_ .resize(n_quant, 0);
234  sum_fluxes_in_ .resize(n_quant, 0);
235  sum_fluxes_out_ .resize(n_quant, 0);
236  sum_masses_ .resize(n_quant, 0);
237  sum_sources_ .resize(n_quant, 0);
238  sum_sources_in_ .resize(n_quant, 0);
239  sum_sources_out_.resize(n_quant, 0);
240 
241  if (cumulative_)
242  {
243  initial_mass_ .resize(n_quant, 0);
244  integrated_fluxes_ .resize(n_quant, 0);
245  integrated_sources_.resize(n_quant, 0);
246  increment_sources_.resize(n_quant, 0);
247  increment_fluxes_.resize(n_quant, 0);
248  }
249 
250 
251 
252  region_mass_matrix_ = new Mat[n_quant];
253  be_flux_matrix_ = new Mat[n_quant];
254  region_source_matrix_ = new Mat[n_quant];
255  region_source_rhs_ = new Mat[n_quant];
256  region_mass_vec_ = new Vec[n_quant];
257  be_flux_vec_ = new Vec[n_quant];
258  region_source_vec_ = new Vec[n_quant];
259 
260  for (unsigned int c=0; c<n_quant; ++c)
261  {
262  chkerr(MatCreateAIJ(PETSC_COMM_WORLD,
263  n_loc_dofs_,
264  (rank_==0)?mesh_->region_db().bulk_size():0,
265  PETSC_DECIDE,
266  PETSC_DECIDE,
267  (rank_==0)?n_bulk_regs_per_dof:0,
268  0,
269  (rank_==0)?0:n_bulk_regs_per_dof,
270  0,
271  &(region_mass_matrix_[c])));
272 
273  chkerr(MatCreateAIJ(PETSC_COMM_WORLD,
274  n_loc_dofs_,
275  (rank_==0)?mesh_->region_db().bulk_size():0,
276  PETSC_DECIDE,
277  PETSC_DECIDE,
278  (rank_==0)?n_bulk_regs_per_dof:0,
279  0,
280  (rank_==0)?0:n_bulk_regs_per_dof,
281  0,
282  &(region_source_matrix_[c])));
283 
284  chkerr(MatCreateAIJ(PETSC_COMM_WORLD,
285  n_loc_dofs_,
286  (rank_==0)?mesh_->region_db().bulk_size():0,
287  PETSC_DECIDE,
288  PETSC_DECIDE,
289  (rank_==0)?n_bulk_regs_per_dof:0,
290  0,
291  (rank_==0)?0:n_bulk_regs_per_dof,
292  0,
293  &(region_source_rhs_[c])));
294 
295  chkerr(MatCreateAIJ(PETSC_COMM_WORLD,
296  be_regions_.size(), // n local rows, number of local boundary edges
297  n_loc_dofs_, // n local cols (local rows of multiplying vector)
298  PETSC_DECIDE, // n global rows
299  PETSC_DECIDE, // n global cols
300  max_dofs_per_boundary_, // allocation, local poriton
301  0,
302  0,
303  0,
304  &(be_flux_matrix_[c])));
305 
306  chkerr(VecCreateMPI(PETSC_COMM_WORLD,
307  (rank_==0)?mesh_->region_db().bulk_size():0,
308  PETSC_DECIDE,
309  &(region_mass_vec_[c])));
310 
311  chkerr(VecCreateMPI(PETSC_COMM_WORLD,
312  be_regions_.size(),
313  PETSC_DECIDE,
314  &(be_flux_vec_[c])));
315 
316  chkerr(VecCreateMPI(PETSC_COMM_WORLD,
317  (rank_==0)?mesh_->region_db().bulk_size():0,
318  PETSC_DECIDE,
319  &(region_source_vec_[c])));
320  }
321 
322  chkerr(MatCreateAIJ(PETSC_COMM_WORLD,
323  be_regions_.size(),
324  (rank_==0)?mesh_->region_db().boundary_size():0,
325  PETSC_DECIDE,
326  PETSC_DECIDE,
327  (rank_==0)?1:0,
328  0,
329  (rank_==0)?0:1,
330  0,
332  ));
333  chkerr(VecGetOwnershipRange(be_flux_vec_[0], &be_offset_, NULL));
334  for (unsigned int loc_el=0; loc_el<be_regions_.size(); ++loc_el)
335  {
336  chkerr(MatSetValue(region_be_matrix_,
337  be_offset_+loc_el,
338  be_regions_[loc_el],
339  1,
340  INSERT_VALUES));
341  }
342  chkerr(MatAssemblyBegin(region_be_matrix_, MAT_FINAL_ASSEMBLY));
343  chkerr(MatAssemblyEnd(region_be_matrix_, MAT_FINAL_ASSEMBLY));
344 
345  double *ones_array;
346  chkerr(VecCreateMPI(PETSC_COMM_WORLD,
347  n_loc_dofs_,
348  PETSC_DECIDE,
349  &ones_));
350  chkerr(VecGetArray(ones_, &ones_array));
351  fill_n(ones_array, n_loc_dofs_, 1);
352  chkerr(VecRestoreArray(ones_, &ones_array));
353 
354  chkerr(VecCreateMPI(PETSC_COMM_WORLD,
355  be_regions_.size(),
356  PETSC_DECIDE,
357  &ones_be_));
358  chkerr(VecGetArray(ones_be_, &ones_array));
359  fill_n(ones_array, be_regions_.size(), 1);
360  chkerr(VecRestoreArray(ones_be_, &ones_array));
361 
362 
363  if (rank_ == 0) {
364  // set default value by output_format_
365  std::string default_file_name;
366  switch (output_format_)
367  {
368  case txt:
369  default_file_name = file_prefix_ + "_balance.txt";
370  break;
371  case gnuplot:
372  default_file_name = file_prefix_ + "_balance.dat";
373  break;
374  case legacy:
375  default_file_name = file_prefix_ + "_balance.txt";
376  break;
377  }
378 
380  try {
382  } INPUT_CATCH(FilePath::ExcFileOpen, FilePath::EI_Address_String, input_record_)
383 
384 
385  // set file name of YAML output
386  string yaml_file_name = file_prefix_ + "_balance.yaml";
388  }
389 
390  allocation_done_ = true;
391 }
392 
393 
394 
396 {
397  if (! balance_on_) return false;
398 
399  auto &marks = TimeGovernor::marks();
400  bool res = marks.current(time_->step(), balance_output_type_) != marks.end(balance_output_type_);
401 
402  //cout << "flag: " << res << " time: " << step.end() << " type:" << hex << balance_output_type_ << endl;
403  return res;
404 }
405 
406 void Balance::start_mass_assembly(unsigned int quantity_idx)
407 {
408  lazy_initialize();
409  if (! balance_on_) return;
410  chkerr(MatZeroEntries(region_mass_matrix_[quantity_idx]));
411  chkerr(VecZeroEntries(region_mass_vec_[quantity_idx]));
412 }
413 
414 
415 void Balance::start_flux_assembly(unsigned int quantity_idx)
416 {
417  lazy_initialize();
418  if (! balance_on_) return;
419  chkerr(MatZeroEntries(be_flux_matrix_[quantity_idx]));
420  chkerr(VecZeroEntries(be_flux_vec_[quantity_idx]));
421 }
422 
423 
424 void Balance::start_source_assembly(unsigned int quantity_idx)
425 {
426  lazy_initialize();
427  if (! balance_on_) return;
428  chkerr(MatZeroEntries(region_source_matrix_[quantity_idx]));
429  chkerr(MatZeroEntries(region_source_rhs_[quantity_idx]));
430  chkerr(VecZeroEntries(region_source_vec_[quantity_idx]));
431 }
432 
433 
434 void Balance::finish_mass_assembly(unsigned int quantity_idx)
435 {
437  if (! balance_on_) return;
438 
439  chkerr(MatAssemblyBegin(region_mass_matrix_[quantity_idx], MAT_FINAL_ASSEMBLY));
440  chkerr(MatAssemblyEnd(region_mass_matrix_[quantity_idx], MAT_FINAL_ASSEMBLY));
441  chkerr(VecAssemblyBegin(region_mass_vec_[quantity_idx]));
442  chkerr(VecAssemblyEnd(region_mass_vec_[quantity_idx]));
443 }
444 
445 void Balance::finish_flux_assembly(unsigned int quantity_idx)
446 {
448  if (! balance_on_) return;
449 
450  chkerr(MatAssemblyBegin(be_flux_matrix_[quantity_idx], MAT_FINAL_ASSEMBLY));
451  chkerr(MatAssemblyEnd(be_flux_matrix_[quantity_idx], MAT_FINAL_ASSEMBLY));
452  chkerr(VecAssemblyBegin(be_flux_vec_[quantity_idx]));
453  chkerr(VecAssemblyEnd(be_flux_vec_[quantity_idx]));
454 }
455 
456 void Balance::finish_source_assembly(unsigned int quantity_idx)
457 {
459  if (! balance_on_) return;
460 
461  chkerr(MatAssemblyBegin(region_source_matrix_[quantity_idx], MAT_FINAL_ASSEMBLY));
462  chkerr(MatAssemblyEnd(region_source_matrix_[quantity_idx], MAT_FINAL_ASSEMBLY));
463  chkerr(MatAssemblyBegin(region_source_rhs_[quantity_idx], MAT_FINAL_ASSEMBLY));
464  chkerr(MatAssemblyEnd(region_source_rhs_[quantity_idx], MAT_FINAL_ASSEMBLY));
465  chkerr(MatMultTranspose(region_source_rhs_[quantity_idx], ones_, region_source_vec_[quantity_idx]));
466 }
467 
468 
469 
470 
471 void Balance::add_mass_matrix_values(unsigned int quantity_idx,
472  unsigned int region_idx,
473  const vector<IdxInt> &dof_indices,
474  const vector<double> &values)
475 {
477  if (! balance_on_) return;
478 
479  PetscInt reg_array[1] = { (int)region_idx };
480 
481  chkerr_assert(MatSetValues(region_mass_matrix_[quantity_idx],
482  dof_indices.size(),
483  &(dof_indices[0]),
484  1,
485  reg_array,
486  &(values[0]),
487  ADD_VALUES));
488 }
489 
490 
491 void Balance::add_flux_matrix_values(unsigned int quantity_idx,
492  unsigned int boundary_idx,
493  const vector<IdxInt> &dof_indices,
494  const vector<double> &values)
495 {
497  if (! balance_on_) return;
498 
499  PetscInt elem_array[1] = { int(be_offset_+boundary_idx) };
500  chkerr_assert(MatSetValues(be_flux_matrix_[quantity_idx],
501  1,
502  elem_array,
503  dof_indices.size(),
504  &(dof_indices[0]),
505  &(values[0]),
506  ADD_VALUES));
507 }
508 
509 
510 void Balance::add_source_matrix_values(unsigned int quantity_idx,
511  unsigned int region_idx,
512  const vector<IdxInt> &dof_indices,
513  const vector<double> &values)
514 {
516  if (! balance_on_) return;
517 
518  PetscInt reg_array[1] = { (int)region_idx };
519 
520  chkerr_assert(MatSetValues(region_source_matrix_[quantity_idx],
521  dof_indices.size(),
522  &(dof_indices[0]),
523  1,
524  reg_array,
525  &(values[0]),
526  ADD_VALUES));
527 }
528 
529 
530 void Balance::add_mass_vec_value(unsigned int quantity_idx,
531  unsigned int region_idx,
532  double value)
533 {
534  chkerr_assert(VecSetValue(region_mass_vec_[quantity_idx],
535  region_idx,
536  value,
537  ADD_VALUES));
538 }
539 
540 
541 void Balance::add_flux_vec_value(unsigned int quantity_idx,
542  unsigned int boundary_idx,
543  double value)
544 {
546  if (! balance_on_) return;
547 
548  chkerr_assert(VecSetValue(be_flux_vec_[quantity_idx],
549  be_offset_+boundary_idx,
550  value,
551  ADD_VALUES));
552 }
553 
554 
555 void Balance::add_source_vec_values(unsigned int quantity_idx,
556  unsigned int region_idx,
557  const vector<IdxInt> &dof_indices,
558  const vector<double> &values)
559 {
561  if (! balance_on_) return;
562 
563  PetscInt reg_array[1] = { (int)region_idx };
564 
565  chkerr_assert(MatSetValues(region_source_rhs_[quantity_idx],
566  dof_indices.size(),
567  &(dof_indices[0]),
568  1,
569  reg_array,
570  &(values[0]),
571  ADD_VALUES));
572 }
573 
574 
575 void Balance::add_cumulative_source(unsigned int quantity_idx, double source)
576 {
578  if (!cumulative_) return;
579 
580  if (rank_ == 0)
581  increment_sources_[quantity_idx] += source;
582 }
583 
584 
585 void Balance::calculate_cumulative(unsigned int quantity_idx,
586  const Vec &solution)
587 {
589  if (!cumulative_) return;
590  if (time_->tlevel() <= 0) return;
591 
592  // sources
593  Vec bulk_vec;
594 
595  chkerr(VecCreateMPIWithArray(PETSC_COMM_WORLD,
596  1,
597  (rank_==0)?mesh_->region_db().bulk_size():0,
598  PETSC_DECIDE,
599  &(sources_[quantity_idx][0]),
600  &bulk_vec));
601 
602  // compute sources on bulk regions: S'.u + s
603  chkerr(VecZeroEntries(bulk_vec));
604  chkerr(MatMultTransposeAdd(region_source_matrix_[quantity_idx],
605  solution, region_source_vec_[quantity_idx], bulk_vec));
606 
607  double sum_sources;
608  chkerr(VecSum(bulk_vec, &sum_sources));
609  VecDestroy(&bulk_vec);
610 
611  if (rank_ == 0)
612  // sum sources in one step
613  increment_sources_[quantity_idx] += sum_sources*time_->dt();
614 
615 
616  // fluxes
617  Vec boundary_vec;
618 
619  chkerr(VecCreateMPIWithArray(PETSC_COMM_WORLD,
620  1,
621  (rank_==0)?mesh_->region_db().boundary_size():0,
622  PETSC_DECIDE,
623  &(fluxes_[quantity_idx][0]),
624  &boundary_vec));
625 
626  // compute fluxes on boundary regions: R'.(F.u + f)
627  chkerr(VecZeroEntries(boundary_vec));
628  Vec temp;
629  chkerr(VecDuplicate(ones_be_, &temp));
630  chkerr(MatMultAdd(be_flux_matrix_[quantity_idx], solution, be_flux_vec_[quantity_idx], temp));
631  // Since internally we keep outgoing fluxes, we change sign
632  // to write to output _incoming_ fluxes.
633  chkerr(VecScale(temp, -1));
634  chkerr(MatMultTranspose(region_be_matrix_, temp, boundary_vec));
635  chkerr(VecDestroy(&temp));
636 
637  double sum_fluxes;
638  chkerr(VecSum(boundary_vec, &sum_fluxes));
639  chkerr(VecDestroy(&boundary_vec));
640 
641  if (rank_ == 0)
642  // sum fluxes in one step
643  increment_fluxes_[quantity_idx] += sum_fluxes*time_->dt();
644 }
645 
646 
647 void Balance::calculate_mass(unsigned int quantity_idx,
648  const Vec &solution,
649  vector<double> &output_array)
650 {
652  if (! balance_on_) return;
653 
654  Vec bulk_vec;
655 
656  chkerr(VecCreateMPIWithArray(PETSC_COMM_WORLD,
657  1,
658  (rank_==0)?mesh_->region_db().bulk_size():0,
659  PETSC_DECIDE,
660  &(output_array[0]),
661  &bulk_vec));
662 
663  // compute mass on regions: M'.u
664  chkerr(VecZeroEntries(bulk_vec));
665  chkerr(MatMultTransposeAdd(region_mass_matrix_[quantity_idx],
666  solution,
667  region_mass_vec_[quantity_idx],
668  bulk_vec));
669  VecDestroy(&bulk_vec);
670 }
671 
672 void Balance::calculate_instant(unsigned int quantity_idx, const Vec& solution)
673 {
674  if ( !is_current() ) return;
675 
676  calculate_mass(quantity_idx, solution, masses_[quantity_idx]);
677 
678  // calculate source
679  Vec bulk_vec;
680 
681  chkerr(VecCreateMPIWithArray(PETSC_COMM_WORLD,
682  1,
683  (rank_==0)?mesh_->region_db().bulk_size():0,
684  PETSC_DECIDE,
685  &(sources_[quantity_idx][0]),
686  &bulk_vec));
687 
688  // compute sources on bulk regions: S'.u + s
689  chkerr(VecZeroEntries(bulk_vec));
690  chkerr(MatMultTransposeAdd(region_source_matrix_[quantity_idx],
691  solution,
692  region_source_vec_[quantity_idx],
693  bulk_vec));
694 
695  // compute positive/negative sources
696  int lsize;
697  Vec mat_r, rhs_r;
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));
703  for (unsigned int r=0; r<mesh_->region_db().bulk_size(); ++r)
704  {
705  MatGetColumnVector(region_source_matrix_[quantity_idx], mat_r, r);
706  MatGetColumnVector(region_source_rhs_[quantity_idx], rhs_r, r);
707 
708  VecGetArrayRead(mat_r, &mat_array);
709  VecGetArrayRead(rhs_r, &rhs_array);
710 
711  sources_in_[quantity_idx][r] = 0;
712  sources_out_[quantity_idx][r] = 0;
713  for (int i=0; i<lsize; ++i)
714  {
715  double f = mat_array[i]*sol_array[i] + rhs_array[i];
716  if (f > 0) sources_in_[quantity_idx][r] += f;
717  else sources_out_[quantity_idx][r] += f;
718  }
719 
720  VecRestoreArrayRead(mat_r, &mat_array);
721  VecRestoreArrayRead(rhs_r, &rhs_array);
722  }
723  chkerr(VecRestoreArrayRead(solution, &sol_array));
724  VecDestroy(&rhs_r);
725  VecDestroy(&mat_r);
726  VecDestroy(&bulk_vec);
727 
728 
729  // calculate flux
730  Vec boundary_vec;
731 
732  chkerr(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1,
733  (rank_==0)?mesh_->region_db().boundary_size():0,
734  PETSC_DECIDE, &(fluxes_[quantity_idx][0]), &boundary_vec));
735 
736  // compute fluxes on boundary regions: R'.(F.u + f)
737  chkerr(VecZeroEntries(boundary_vec));
738  Vec temp;
739  chkerr(VecDuplicate(ones_be_, &temp));
740  chkerr(MatMultAdd(be_flux_matrix_[quantity_idx], solution, be_flux_vec_[quantity_idx], temp));
741  // Since internally we keep outgoing fluxes, we change sign
742  // to write to output _incoming_ fluxes.
743  chkerr(VecScale(temp, -1));
744  chkerr(MatMultTranspose(region_be_matrix_, temp, boundary_vec));
745 
746  // compute positive/negative fluxes
747  fluxes_in_[quantity_idx].assign(mesh_->region_db().boundary_size(), 0);
748  fluxes_out_[quantity_idx].assign(mesh_->region_db().boundary_size(), 0);
749  const double *flux_array;
750 // int lsize;
751  chkerr(VecGetArrayRead(temp, &flux_array));
752  chkerr(VecGetLocalSize(temp, &lsize));
753  for (int e=0; e<lsize; ++e)
754  {
755  if (flux_array[e] < 0)
756  fluxes_out_[quantity_idx][be_regions_[e]] += flux_array[e];
757  else
758  fluxes_in_[quantity_idx][be_regions_[e]] += flux_array[e];
759  }
760  chkerr(VecRestoreArrayRead(temp, &flux_array));
761  VecDestroy(&temp);
762  VecDestroy(&boundary_vec);
763 }
764 
765 
766 
767 
768 
770 {
772  if (! balance_on_) return;
773  if (! is_current() ) return;
774 
775  // gather results from processes and sum them up
776  const unsigned int n_quant = quantities_.size();
777  const unsigned int n_blk_reg = mesh_->region_db().bulk_size();
778  const unsigned int n_bdr_reg = mesh_->region_db().boundary_size();
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++)
782  {
783  for (unsigned int ri=0; ri<n_blk_reg; ri++)
784  {
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];
787  }
788  for (unsigned int ri=0; ri<n_bdr_reg; ri++)
789  {
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];
792  }
793  }
794  MPI_Reduce(&sendbuffer,recvbuffer,buf_size,MPI_DOUBLE,MPI_SUM,0,PETSC_COMM_WORLD);
795  // for other than 0th process update last_time and finish,
796  // on process #0 sum balances over all regions and calculate
797  // cumulative balance over time.
798  if (rank_ == 0)
799  {
800  // update balance vectors
801  for (unsigned int qi=0; qi<n_quant; qi++)
802  {
803  for (unsigned int ri=0; ri<n_blk_reg; ri++)
804  {
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];
807  }
808  for (unsigned int ri=0; ri<n_bdr_reg; ri++)
809  {
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];
812  }
813  }
814  }
815 
816 
817  // The convention for input/output of fluxes is that positive means inward.
818  // Therefore in the following code we switch sign of fluxes.
819  if (rank_ == 0)
820  {
821  sum_fluxes_.assign(n_quant, 0);
822  sum_fluxes_in_.assign(n_quant, 0);
823  sum_fluxes_out_.assign(n_quant, 0);
824  sum_masses_.assign(n_quant, 0);
825  sum_sources_.assign(n_quant, 0);
826  sum_sources_in_.assign(n_quant, 0);
827  sum_sources_out_.assign(n_quant, 0);
828 
829  // sum all boundary fluxes
830  const RegionSet & b_set = mesh_->region_db().get_region_set(".BOUNDARY");
831  for( RegionSet::const_iterator reg = b_set.begin(); reg != b_set.end(); ++reg)
832  {
833  for (unsigned int qi=0; qi<n_quant; qi++)
834  {
835  sum_fluxes_[qi] += fluxes_ [qi][reg->boundary_idx()];
836  sum_fluxes_in_[qi] += fluxes_in_ [qi][reg->boundary_idx()];
837  sum_fluxes_out_[qi] += fluxes_out_[qi][reg->boundary_idx()];
838  }
839  }
840 
841  // sum all volume sources
842  const RegionSet & bulk_set = mesh_->region_db().get_region_set("BULK");
843  for( RegionSet::const_iterator reg = bulk_set.begin(); reg != bulk_set.end(); ++reg)
844  {
845  for (unsigned int qi=0; qi<n_quant; qi++)
846  {
847  sum_masses_[qi] += masses_[qi][reg->bulk_idx()];
848  sum_sources_[qi] += sources_[qi][reg->bulk_idx()];
849  sum_sources_in_[qi] += sources_in_[qi][reg->bulk_idx()];
850  sum_sources_out_[qi] += sources_out_[qi][reg->bulk_idx()];
851  }
852  }
853 
854  // cumulative balance over time
855  if (cumulative_)
856  {
857  // save initial time and mass
858  if (initial_)
859  {
861  for (unsigned int qi=0; qi<n_quant; qi++)
862  initial_mass_[qi] = sum_masses_[qi];
863  initial_ = false;
864  }
865 
866  for (unsigned int qi=0; qi<n_quant; qi++)
867  {
870  }
871  }
872  }
873 
874  last_time_ = time_->t();
875 
876 
877  // perform actual output
878  switch (output_format_)
879  {
880  case txt:
881  output_csv(time_->t(), '\t', "");
882  break;
883  case gnuplot:
884  output_csv(time_->t(), ' ', "#", 30);
885  break;
886  case legacy:
887  output_legacy(time_->t());
888  break;
889  }
890  // output in YAML format
891  output_yaml(time_->t());
892 
893  if (rank_ == 0)
894  {
895  sum_fluxes_.assign(n_quant, 0);
896  sum_sources_.assign(n_quant, 0);
897  increment_fluxes_.assign(n_quant, 0);
898  increment_sources_.assign(n_quant, 0);
899  }
900 }
901 
902 
903 void Balance::output_legacy(double time)
904 {
905  // write output only on process #0
906  if (rank_ != 0) return;
907 
908  const unsigned int n_quant = quantities_.size();
909 
910  // print the head of mass balance file
911  unsigned int c = 6; //column number without label
912  unsigned int w = 14; //column width
913  unsigned int wl = 2*(w-5)+7; //label column width
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";
917 
918  output_ << "# " << setw((w*c+wl-14)/2) << setfill('-') << "--"
919  << " MASS BALANCE "
920  << setw((w*c+wl-14)/2) << setfill('-') << "" << endl
921  << "# Time: " << time << "\n\n\n";
922 
923  // header for table of boundary fluxes
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]"
932  << endl;
933 
934  // draw long line
935  output_ << "# " << setw(w*c+wl) << setfill('-') << "" << setfill(' ') << endl;
936 
937  // print mass fluxes over boundaries
938  const RegionSet & b_set = mesh_->region_db().get_region_set(".BOUNDARY");
939  for( RegionSet::const_iterator reg = b_set.begin(); reg != b_set.end(); ++reg) {
940  for (unsigned int qi=0; qi<n_quant; qi++) {
941  output_ << setw(2) << ""
942  << setw(w) << (int)reg->id()
943  << setw(wl) << reg->label().c_str()
944  << setw(w) << quantities_[qi].name_.c_str()
945  << setw(w) << fluxes_[qi][reg->boundary_idx()]
946  << setw(w) << fluxes_out_[qi][reg->boundary_idx()]
947  << setw(w) << fluxes_in_[qi][reg->boundary_idx()]
948  << endl;
949  }
950  }
951 
952  // draw long line
953  output_ << "# " << setw(w*c+wl) << setfill('-') << "" << setfill(' ') << endl;
954 
955  // total boundary balance
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]"
959  << setw(w) << quantities_[qi].name_.c_str()
960  << setw(w) << sum_fluxes_[qi]
961  << setw(w) << sum_fluxes_out_[qi]
962  << setw(w) << sum_fluxes_in_[qi]
963  << endl;
964  output_ << "\n\n";
965 
966 
967  // header for table of volume sources and masses
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]"
978  << endl;
979 
980  // draw long line
981  output_ << "# " << setw(w*c+wl) << setfill('-') << "" << setfill(' ') << endl;
982 
983  // print balance of volume sources and masses
984  const RegionSet & bulk_set = mesh_->region_db().get_region_set("BULK");
985  for( RegionSet::const_iterator reg = bulk_set.begin(); reg != bulk_set.end(); ++reg)
986  {
987  for (unsigned int qi=0; qi<n_quant; qi++)
988  {
989  output_ << setw(2) << ""
990  << setw(w) << (int)reg->id()
991  << setw(wl) << reg->label().c_str()
992  << setw(w) << quantities_[qi].name_.c_str()
993  << setw(w) << masses_[qi][reg->bulk_idx()]
994  << setw(w) << sources_[qi][reg->bulk_idx()]
995  << endl;
996  }
997  }
998 
999  // draw long line
1000  output_ << "# " << setw(w*c+wl) << setfill('-') << "" << setfill(' ') << endl;
1001 
1002  // total sources balance
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]"
1005  << setw(w) << quantities_[qi].name_.c_str()
1006  << setw(w) << sum_masses_[qi]
1007  << setw(w) << sum_sources_[qi]
1008  << endl;
1009 
1010  if (cumulative_)
1011  {
1012  // Print cumulative sources
1013  output_ << "\n\n# Cumulative mass balance on time interval ["
1014  << setiosflags(ios::left) << time_->init_time() << ","
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]"
1026  << endl;
1027 
1028  for (unsigned int qi=0; qi<n_quant; qi++)
1029  {
1030  double denominator = max(fabs(initial_mass_[qi]+integrated_sources_[qi]-integrated_fluxes_[qi]),fabs(sum_masses_[qi]));
1031  output_ << " " << setiosflags(ios::left)
1032  << setw(w) << quantities_[qi].name_.c_str()
1033  << setw(w) << initial_mass_[qi]
1034  << setw(w) << integrated_sources_[qi]
1035  << setw(w) << integrated_fluxes_[qi]
1036  << setw(w) << initial_mass_[qi]+integrated_sources_[qi]-integrated_fluxes_[qi]
1037  << setw(w) << sum_masses_[qi]
1038  << setw(w) << initial_mass_[qi]+integrated_sources_[qi]-integrated_fluxes_[qi]-sum_masses_[qi]
1039  << setw(w) << fabs(initial_mass_[qi]+integrated_sources_[qi]+integrated_fluxes_[qi]-sum_masses_[qi])/(denominator==0?1:denominator)
1040  << endl;
1041  }
1042  }
1043 
1044  output_ << endl << endl;
1045 }
1046 
1047 
1048 std::string Balance::csv_zero_vals(unsigned int cnt, char delimiter)
1049 {
1050  std::stringstream ss;
1051  for (unsigned int i=0; i<cnt; i++) ss << format_csv_val(0, delimiter);
1052  return ss.str();
1053 }
1054 
1055 
1056 void Balance::output_csv(double time, char delimiter, const std::string& comment_string, unsigned int repeat)
1057 {
1058  // write output only on process #0
1059  if (rank_ != 0) return;
1060 
1061  const unsigned int n_quant = quantities_.size();
1062 
1063  // print data header only on first line
1064  if (repeat==0 && output_line_counter_==0) format_csv_output_header(delimiter, comment_string);
1065 
1066  // print sources and masses over bulk regions
1067  const RegionSet & bulk_set = mesh_->region_db().get_region_set("BULK");
1068  for( RegionSet::const_iterator reg = bulk_set.begin(); reg != bulk_set.end(); ++reg)
1069  {
1070  for (unsigned int qi=0; qi<n_quant; qi++)
1071  {
1072  // print data header (repeat header after every "repeat" lines)
1073  if (repeat && (output_line_counter_%repeat == 0)) format_csv_output_header(delimiter, comment_string);
1074 
1075  output_ << format_csv_val(time, delimiter, true)
1076  << format_csv_val(reg->label(), delimiter)
1077  << format_csv_val(quantities_[qi].name_, delimiter)
1078  << csv_zero_vals(3, delimiter)
1079  << format_csv_val(masses_[qi][reg->bulk_idx()], delimiter)
1080  << format_csv_val(sources_[qi][reg->bulk_idx()], delimiter)
1081  << format_csv_val(sources_in_[qi][reg->bulk_idx()], delimiter)
1082  << format_csv_val(sources_out_[qi][reg->bulk_idx()], delimiter)
1083  << csv_zero_vals(5, delimiter) << endl;
1085  }
1086  }
1087 
1088  // print mass fluxes over boundaries
1089  const RegionSet & b_set = mesh_->region_db().get_region_set(".BOUNDARY");
1090  for( RegionSet::const_iterator reg = b_set.begin(); reg != b_set.end(); ++reg)
1091  {
1092  for (unsigned int qi=0; qi<n_quant; qi++) {
1093  // print data header (repeat header after every "repeat" lines)
1094  if (repeat && (output_line_counter_%repeat == 0)) format_csv_output_header(delimiter, comment_string);
1095 
1096  output_ << format_csv_val(time, delimiter, true)
1097  << format_csv_val(reg->label(), delimiter)
1098  << format_csv_val(quantities_[qi].name_, delimiter)
1099  << format_csv_val(fluxes_[qi][reg->boundary_idx()], delimiter)
1100  << format_csv_val(fluxes_in_[qi][reg->boundary_idx()], delimiter)
1101  << format_csv_val(fluxes_out_[qi][reg->boundary_idx()], delimiter)
1102  << csv_zero_vals(9, delimiter) << endl;
1104  }
1105  }
1106 
1107  if (cumulative_)
1108  {
1109  for (unsigned int qi=0; qi<n_quant; qi++)
1110  {
1111  // print data header (repeat header after every "repeat" lines)
1112  if (repeat && (output_line_counter_%repeat == 0)) format_csv_output_header(delimiter, comment_string);
1113 
1114  double error = sum_masses_[qi] - (initial_mass_[qi] + integrated_sources_[qi] + integrated_fluxes_[qi]);
1115  output_ << format_csv_val(time, delimiter, true)
1116  << format_csv_val("ALL", delimiter)
1117  << format_csv_val(quantities_[qi].name_, delimiter)
1118  << format_csv_val(sum_fluxes_[qi], delimiter)
1119  << format_csv_val(sum_fluxes_in_[qi], delimiter)
1120  << format_csv_val(sum_fluxes_out_[qi], delimiter)
1121  << format_csv_val(sum_masses_[qi], delimiter)
1122  << format_csv_val(sum_sources_[qi], delimiter)
1123  << format_csv_val(sum_sources_in_[qi], delimiter)
1124  << format_csv_val(sum_sources_out_[qi], delimiter)
1125  << format_csv_val(increment_fluxes_[qi], delimiter)
1126  << format_csv_val(increment_sources_[qi], delimiter)
1127  << format_csv_val(integrated_fluxes_[qi], delimiter)
1128  << format_csv_val(integrated_sources_[qi], delimiter)
1129  << format_csv_val(error, delimiter) << endl;
1131  }
1132  }
1133 
1134 }
1135 
1136 
1137 void Balance::format_csv_output_header(char delimiter, const std::string& comment_string)
1138 {
1139  std::stringstream ss;
1140  if (delimiter == ' ') {
1141  ss << setw(output_column_width-comment_string.size()) << "\"time\"";
1142  } else {
1143  ss << "\"time\"";
1144  }
1145 
1146  output_ << comment_string << ss.str()
1147  << format_csv_val("region", delimiter)
1148  << format_csv_val("quantity [" + units_.format_text() + "]", delimiter)
1149  << format_csv_val("flux", delimiter)
1150  << format_csv_val("flux_in", delimiter)
1151  << format_csv_val("flux_out", delimiter)
1152  << format_csv_val("mass", delimiter)
1153  << format_csv_val("source", delimiter)
1154  << format_csv_val("source_in", delimiter)
1155  << format_csv_val("source_out", delimiter)
1156  << format_csv_val("flux_increment", delimiter)
1157  << format_csv_val("source_increment", delimiter)
1158  << format_csv_val("flux_cumulative", delimiter)
1159  << format_csv_val("source_cumulative", delimiter)
1160  << format_csv_val("error", delimiter)
1161  << endl;
1162 }
1163 
1164 std::string Balance::format_csv_val(std::string val, char delimiter, bool initial)
1165 {
1166  std::stringstream ss;
1167  std::replace( val.begin(), val.end(), '\"', '\'');
1168 
1169  if (!initial) ss << delimiter;
1170  if (delimiter == ' ') {
1171  std::stringstream sval;
1172  sval << "\"" << val << "\"";
1173  ss << " " << setw(output_column_width-1) << sval.str();
1174  } else {
1175  ss << "\"" << val << "\"";
1176  }
1177 
1178  return ss.str();
1179 }
1180 
1181 std::string Balance::format_csv_val(double val, char delimiter, bool initial)
1182 {
1183  std::stringstream ss;
1184 
1185  if (!initial) ss << delimiter;
1186  if (delimiter == ' ') {
1187  ss << " " << setw(output_column_width-1) << val;
1188  } else {
1189  ss << val;
1190  }
1191  return ss.str();
1192 }
1193 
1194 void Balance::output_yaml(double time)
1195 {
1196  // write output only on process #0
1197  if (rank_ != 0) return;
1198 
1199  const unsigned int n_quant = quantities_.size();
1200 
1201  // print data header only once
1202  if (!output_yaml_header_) {
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;
1205  output_yaml_ << "data:" << endl;
1206  output_yaml_header_ = true;
1207  }
1208 
1209  output_yaml_ << setfill(' ');
1210 
1211  // print sources and masses over bulk regions
1212  const RegionSet & bulk_set = mesh_->region_db().get_region_set("BULK");
1213  for( RegionSet::const_iterator reg = bulk_set.begin(); reg != bulk_set.end(); ++reg)
1214  {
1215  for (unsigned int qi=0; qi<n_quant; qi++)
1216  {
1217  output_yaml_ << " - time: " << time << endl;
1218  output_yaml_ << setw(4) << "" << "region: " << reg->label() << endl;
1219  output_yaml_ << setw(4) << "" << "quantity: " << quantities_[qi].name_ << endl;
1220  output_yaml_ << setw(4) << "" << "data: " << "[ 0, 0, 0, " << masses_[qi][reg->bulk_idx()] << ", "
1221  << sources_[qi][reg->bulk_idx()] << ", " << sources_in_[qi][reg->bulk_idx()] << ", "
1222  << sources_out_[qi][reg->bulk_idx()] << ", 0, 0, 0, 0, 0 ]" << endl;
1223  }
1224  }
1225 
1226  // print mass fluxes over boundaries
1227  const RegionSet & b_set = mesh_->region_db().get_region_set(".BOUNDARY");
1228  for( RegionSet::const_iterator reg = b_set.begin(); reg != b_set.end(); ++reg)
1229  {
1230  for (unsigned int qi=0; qi<n_quant; qi++) {
1231  output_yaml_ << " - time: " << time << endl;
1232  output_yaml_ << setw(4) << "" << "region: " << reg->label() << endl;
1233  output_yaml_ << setw(4) << "" << "quantity: " << quantities_[qi].name_ << endl;
1234  output_yaml_ << setw(4) << "" << "data: " << "[ " << fluxes_[qi][reg->boundary_idx()] << ", "
1235  << fluxes_in_[qi][reg->boundary_idx()] << ", " << fluxes_out_[qi][reg->boundary_idx()]
1236  << ", 0, 0, 0, 0, 0, 0, 0, 0, 0 ]" << endl;
1237  }
1238  }
1239 
1240  if (cumulative_)
1241  {
1242  for (unsigned int qi=0; qi<n_quant; qi++)
1243  {
1244  double error = sum_masses_[qi] - (initial_mass_[qi] + integrated_sources_[qi] + integrated_fluxes_[qi]);
1245  output_yaml_ << " - time: " << time << endl;
1246  output_yaml_ << setw(4) << "" << "region: ALL" << endl;
1247  output_yaml_ << setw(4) << "" << "quantity: " << quantities_[qi].name_ << endl;
1248  output_yaml_ << setw(4) << "" << "data: " << "[ " << sum_fluxes_[qi] << ", "
1249  << sum_fluxes_in_[qi] << ", " << sum_fluxes_out_[qi] << ", "
1250  << sum_masses_[qi] << ", " << sum_sources_[qi] << ", "
1251  << sum_sources_in_[qi] << ", " << sum_sources_out_[qi] << ", "
1252  << increment_fluxes_[qi] << ", " << increment_sources_[qi] << ", "
1253  << integrated_fluxes_[qi] << ", " << integrated_sources_[qi] << ", "
1254  << error << " ]" << endl;
1255  }
1256  }
1257 }
1258 
1259 
1260 
1261 
1262 
1263 
UnitSI units_
Units of conserved quantities.
Definition: balance.hh:445
void lazy_initialize()
Definition: balance.cc:184
static const Input::Type::Record & get_input_type()
Main balance input record type.
Definition: balance.cc:42
std::vector< double > integrated_sources_
Definition: balance.hh:506
unsigned int add_quantity(const string &name)
Definition: balance.cc:156
std::vector< double > sum_fluxes_out_
Definition: balance.hh:498
void add_mass_matrix_values(unsigned int quantity_idx, unsigned int region_idx, const std::vector< IdxInt > &dof_indices, const std::vector< double > &values)
Definition: balance.cc:471
void calculate_mass(unsigned int quantity_idx, const Vec &solution, vector< double > &output_array)
Definition: balance.cc:647
TimeMark::Type output_mark_type_
TimeMark type for output of particular equation.
Definition: balance.hh:518
Accessor to input data conforming to declared Array.
Definition: accessors.hh:567
void allocate(unsigned int n_loc_dofs, unsigned int max_dofs_per_boundary)
Definition: balance.cc:176
TimeMark::Type balance_output_type_
TimeMark type for balance output of particular equation.
Definition: balance.hh:515
unsigned int * boundary_idx_
Definition: elements.h:94
std::string format_text() const
Definition: unit_si.cc:127
int tlevel() const
void finish_source_assembly(unsigned int quantity_idx)
This method must be called after assembling the matrix and vectors for computing source.
Definition: balance.cc:456
unsigned int bulk_size() const
Definition: region.cc:268
void add_mass_vec_value(unsigned int quantity_idx, unsigned int region_idx, double value)
Definition: balance.cc:530
Class Input::Type::Default specifies default value of keys of a Input::Type::Record.
Definition: type_record.hh:61
Class for declaration of the input of type Bool.
Definition: type_base.hh:458
gnuplot
Definition: balance.hh:134
bool initial_
true before calculating the mass at initial time, otherwise false
Definition: balance.hh:521
Boundary * cond() const
Definition: side_impl.hh:70
IdxInt * get_el_4_loc() const
Definition: mesh.h:195
std::vector< std::vector< double > > sources_in_
Definition: balance.hh:492
void add_cumulative_source(unsigned int quantity_idx, double source)
Definition: balance.cc:575
std::vector< std::vector< double > > fluxes_in_
Definition: balance.hh:488
std::vector< double > sum_fluxes_in_
Definition: balance.hh:497
std::vector< std::vector< double > > fluxes_out_
Definition: balance.hh:489
#define INPUT_CATCH(ExceptionType, AddressEITag, input_accessor)
Definition: accessors.hh:64
const TimeGovernor * time_
Definition: balance.hh:548
void add_source_vec_values(unsigned int quantity_idx, unsigned int region_idx, const std::vector< IdxInt > &dof_values, const std::vector< double > &values)
Definition: balance.cc:555
Vec ones_
auxiliary vectors for summation of matrix columns
Definition: balance.hh:476
#define MPI_SUM
Definition: mpi.h:196
Mat * region_source_rhs_
Matrices for calculation of signed source (n_dofs x n_bulk_regions).
Definition: balance.hh:458
RegionSet get_region_set(const string &set_name) const
Definition: region.cc:328
Definition: mesh.h:99
Input::Record input_record_
Record for current balance.
Definition: balance.hh:546
ofstream output_
Handle for file for output in given OutputFormat of balance and total fluxes over individual regions ...
Definition: balance.hh:431
Vec ones_be_
Definition: balance.hh:476
void read_from_input(Input::Array in_array, const TimeGovernor &tg)
#define FOR_ELEMENT_SIDES(i, j)
Definition: elements.h:157
void chkerr(unsigned int ierr)
Replacement of new/delete operator in the spirit of xmalloc.
Definition: system.hh:147
void output_legacy(double time)
Perform output in old format (for compatibility)
Definition: balance.cc:903
Mat region_be_matrix_
Definition: balance.hh:473
const RegionDB & region_db() const
Definition: mesh.h:170
#define ASSERT(expr)
Allow use shorter versions of macro names if these names is not used with external library...
Definition: asserts.hh:346
Region region()
Definition: boundaries.cc:34
const TimeStep & step(int index=-1) const
bool allocation_done_
true before allocating necessary internal structures (Petsc matrices etc.)
Definition: balance.hh:527
double t() const
void calculate_cumulative(unsigned int quantity_idx, const Vec &solution)
Definition: balance.cc:585
Basic time management functionality for unsteady (and steady) solvers (class Equation).
static TimeMarks & marks()
#define MPI_Reduce(sendbuf, recvbuf, count, datatype, op, root, comm)
Definition: mpi.h:608
void calculate_instant(unsigned int quantity_idx, const Vec &solution)
Definition: balance.cc:672
void add_flux_vec_value(unsigned int quantity_idx, unsigned int boundary_idx, double value)
Definition: balance.cc:541
Record & close() const
Close the Record for further declarations of keys.
Definition: type_record.cc:303
Basic time management class.
std::vector< double > sum_fluxes_
Definition: balance.hh:496
bool balance_on_
If the balance is on. Balance is off in the case of no balance output time marks. ...
Definition: balance.hh:530
static constexpr bool value
Definition: json.hpp:87
std::vector< double > increment_sources_
Definition: balance.hh:509
Vec * be_flux_vec_
Vectors for calculation of flux (n_boundary_edges).
Definition: balance.hh:461
static const Input::Type::Selection & get_format_selection_input_type()
Input selection for file format.
Definition: balance.cc:34
std::vector< double > integrated_fluxes_
Definition: balance.hh:507
unsigned int boundary_idx() const
Returns index of the region in the boundary set.
Definition: region.hh:85
std::vector< Quantity > quantities_
Names of conserved quantities.
Definition: balance.hh:440
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...
Definition: balance.cc:1048
std::vector< std::vector< double > > sources_out_
Definition: balance.hh:493
Mat * be_flux_matrix_
Matrices for calculation of flux (n_boundary_edges x n_dofs).
Definition: balance.hh:452
double last_time_
time of last calculated balance
Definition: balance.hh:512
void open_stream(Stream &stream) const
Definition: file_path.cc:211
double init_time() const
std::vector< std::vector< double > > fluxes_
Definition: balance.hh:487
void finish_mass_assembly(unsigned int quantity_idx)
This method must be called after assembling the matrix for computing mass.
Definition: balance.cc:434
TimeMark::Type equation_fixed_mark_type() const
void format_csv_output_header(char delimiter, const std::string &comment_string)
Print output header.
Definition: balance.cc:1137
int rank_
MPI rank.
Definition: balance.hh:537
Accessor to the data with type Type::Record.
Definition: accessors.hh:292
TimeMark::Type equation_mark_type() const
const Ret val(const string &key) const
void init_from_input(const Input::Record &in_rec, TimeGovernor &tg)
Definition: balance.cc:123
Mat * region_mass_matrix_
Matrices for calculation of mass (n_dofs x n_bulk_regions).
Definition: balance.hh:449
std::vector< std::vector< double > > sources_
Definition: balance.hh:491
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.
std::vector< double > sum_sources_out_
Definition: balance.hh:502
unsigned int max_dofs_per_boundary_
Definition: balance.hh:421
bool cumulative_
if true then cumulative balance is computed
Definition: balance.hh:524
void add_source_matrix_values(unsigned int quantity_idx, unsigned int region_idx, const std::vector< IdxInt > &dof_indices, const std::vector< double > &values)
Definition: balance.cc:510
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:490
Mat * region_source_matrix_
Matrices for calculation of source (n_dofs x n_bulk_regions).
Definition: balance.hh:455
std::vector< double > increment_fluxes_
Definition: balance.hh:508
void chkerr_assert(unsigned int ierr)
Definition: system.hh:158
SideIter side(const unsigned int loc_index)
bool is_current()
Returns true if the current time step is marked for the balance output.
Definition: balance.cc:395
Balance(const std::string &file_prefix, const Mesh *mesh)
Definition: balance.cc:74
void output_csv(double time, char delimiter, const std::string &comment_string, unsigned int repeat=0)
Perform output in csv format.
Definition: balance.cc:1056
Distribution * get_el_ds() const
Definition: mesh.h:189
~Balance()
Definition: balance.cc:91
#define MPI_DOUBLE
Definition: mpi.h:156
#define MPI_Comm_rank
Definition: mpi.h:236
std::vector< unsigned int > add_quantities(const std::vector< string > &names)
Definition: balance.cc:167
Dedicated class for storing path to input and output files.
Definition: file_path.hh:54
const Mesh * mesh_
Definition: balance.hh:442
static Default read_time(const std::string &description)
The factory function to make an default value that will be specified at the time when a key will be r...
Definition: type_record.hh:97
void output_yaml(double time)
Perform output in yaml format.
Definition: balance.cc:1194
void output()
Perform output to file for given time instant.
Definition: balance.cc:769
std::vector< double > sum_sources_
Definition: balance.hh:500
void start_flux_assembly(unsigned int quantity_idx)
Definition: balance.cc:415
std::vector< double > initial_mass_
Definition: balance.hh:503
static const unsigned int output_column_width
Size of column in output (used if delimiter is space)
Definition: balance.hh:387
#define ASSERT_PTR(ptr)
Definition of assert macro checking non-null pointer (PTR)
Definition: asserts.hh:335
const Selection & close() const
Close the Selection, no more values can be added.
unsigned int boundary_size() const
Definition: region.cc:261
double dt() const
void finish_flux_assembly(unsigned int quantity_idx)
This method must be called after assembling the matrix and vector for computing flux.
Definition: balance.cc:445
std::vector< unsigned int > be_regions_
Number of boundary region for each local boundary edge.
Definition: balance.hh:479
ofstream output_yaml_
Definition: balance.hh:434
OutputFormat output_format_
Format of output file.
Definition: balance.hh:437
#define ASSERT_DBG(expr)
Definition: asserts.hh:349
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)
Definition: balance.cc:491
void units(const UnitSI &unit)
Setter for units of conserved quantities.
Definition: balance.cc:148
int be_offset_
Offset for local part of vector of boundary edges.
Definition: balance.hh:482
std::vector< std::vector< double > > masses_
Definition: balance.hh:490
Vec * region_source_vec_
Vectors for calculation of source (n_bulk_regions).
Definition: balance.hh:467
void start_mass_assembly(unsigned int quantity_idx)
Definition: balance.cc:406
Record type proxy class.
Definition: type_record.hh:182
std::vector< double > sum_masses_
Definition: balance.hh:499
Vec * region_mass_vec_
Vectors for calculation of mass (n_bulk_regions).
Definition: balance.hh:464
unsigned int output_line_counter_
hold count of line printed into output_
Definition: balance.hh:540
bool add_output_times_
Add output time marks to balance output time marks.
Definition: balance.hh:533
Class for representation SI units of Fields.
Definition: unit_si.hh:40
unsigned int n_loc_dofs_
Allocation parameters. Set by the allocate method used in the lazy_initialize.
Definition: balance.hh:420
static FileName output()
The factory function for declaring type FileName for input files.
Definition: type_base.cc:533
OutputFormat
Definition: balance.hh:130
void undef(bool val=true)
Definition: unit_si.cc:197
void start_source_assembly(unsigned int quantity_idx)
Definition: balance.cc:424
Template for classes storing finite set of named values.
FilePath balance_output_file_
File path for output_ stream.
Definition: balance.hh:428
std::string file_prefix_
Save prefix passed in in constructor.
Definition: balance.hh:425
bool output_yaml_header_
marks whether YAML output has printed header
Definition: balance.hh:543
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.
Definition: balance.cc:1164
ElementVector element
Vector of elements of the mesh.
Definition: mesh.h:260
std::vector< double > sum_sources_in_
Definition: balance.hh:501
unsigned int lsize(int proc) const
get local size