Flow123d  master-f44eb46
darcy_flow_mh_output.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 darcy_flow_mh_output.cc
15  * @ingroup flow
16  * @brief Output class for darcy_flow_mh model.
17  * @author Jan Brezina
18  */
19 
20 #include <vector>
21 #include <iostream>
22 #include <sstream>
23 #include <string>
24 
25 #include <system/global_defs.h>
26 
27 #include "flow/darcy_flow_lmh.hh"
28 #include "flow/assembly_lmh.hh"
32 
33 #include "io/output_time.hh"
34 #include "io/observe.hh"
35 #include "system/system.hh"
36 #include "system/sys_profiler.hh"
37 #include "system/index_types.hh"
38 
39 #include "fields/field_set.hh"
40 #include "fem/dofhandler.hh"
41 #include "fields/field_fe.hh"
42 #include "fields/generic_field.hh"
43 
44 #include "mesh/mesh.h"
45 #include "mesh/partitioning.hh"
46 #include "mesh/accessors.hh"
47 #include "mesh/node_accessor.hh"
48 #include "mesh/range_wrapper.hh"
50 
51 // #include "coupling/balance.hh"
54 
55 namespace it = Input::Type;
56 
57 
58 const it::Instance & DarcyFlowMHOutput::get_input_type(FieldSet& eq_data, const std::string &equation_name) {
60  output_fields += eq_data;
61  return output_fields.make_output_type(equation_name, "");
62 }
63 
64 
66 
67  static it::Record& rec = it::Record("Output_DarcyMHSpecific", "Specific Darcy flow MH output.")
69  .declare_key("compute_errors", it::Bool(), it::Default("false"),
70  "SPECIAL PURPOSE. Computes error norms of the solution, particulary suited for non-compatible coupling models.")
72  "Output file with raw data from MH module.")
73  .close();
74 
77  "Flow_Darcy_MH_specific",
78  "");
79 }
80 
81 
84 {
85 
86 // *this += field_ele_pressure.name("pressure_p0_old").units(UnitSI().m()) // TODO remove: obsolete field
87 // .flags(FieldFlag::equation_result)
88 // .description("Pressure solution - P0 interpolation.");
89 // *this += field_node_pressure.name("pressure_p1").units(UnitSI().m())
90 // .flags(FieldFlag::equation_result)
91 // .description("Pressure solution - P1 interpolation.");
92 // *this += field_ele_piezo_head.name("piezo_head_p0_old").units(UnitSI().m()) // TODO remove: obsolete field
93 // .flags(FieldFlag::equation_result)
94 // .description("Piezo head solution - P0 interpolation.");
95 // *this += field_ele_flux.name("velocity_p0_old").units(UnitSI().m()) // TODO remove: obsolete field
96 // .flags(FieldFlag::equation_result)
97 // .description("Velocity solution - P0 interpolation.");
98  *this += subdomain.name("subdomain")
101  .description("Subdomain ids of the domain decomposition.");
102  *this += region_id.name("region_id")
105  .description("Region ids.");
106 }
107 
108 
110 : EquationOutput()
111 {
112  *this += pressure_diff.name("pressure_diff").units(UnitSI().m())
114  .description("Error norm of the pressure solution. [Experimental]");
115  *this += velocity_diff.name("velocity_diff").units(UnitSI().m().s(-1))
117  .description("Error norm of the velocity solution. [Experimental]");
118  *this += div_diff.name("div_diff").units(UnitSI().s(-1))
120  .description("Error norm of the divergence of the velocity solution. [Experimental]");
121 }
122 
124 : darcy_flow(flow),
125  mesh_(&darcy_flow->mesh()),
126  compute_errors_(false),
128  l2_difference_assembly_(nullptr),
130 {
132  main_mh_in_rec.val<Input::Record>("output_stream"),
134  prepare_output(main_mh_in_rec);
135 
137  raw_eq_data_ = std::make_shared<RawOutputEqData>();
138  raw_eq_data_->flow_data_ = darcy_flow->eq_data_;
139  raw_eq_data_->time_ = &darcy_flow->time();
140  ASSERT_PTR(raw_eq_data_->flow_data_);
141 
142  auto in_rec_specific = main_mh_in_rec.find<Input::Record>("output_specific");
143  if (in_rec_specific) {
144  in_rec_specific->opt_val("compute_errors", compute_errors_);
145 
146  // raw output
147  int rank;
149  if (rank == 0) {
150 
151  // optionally open raw output file
152  FilePath raw_output_file_path;
153  if (in_rec_specific->opt_val("raw_flow_output", raw_output_file_path))
154  {
155  int mpi_size;
156  MPI_Comm_size(MPI_COMM_WORLD, &mpi_size);
157  if(mpi_size > 1)
158  {
159  WarningOut() << "Raw output is not available in parallel computation. MPI size: " << mpi_size << "\n";
160  }
161  else
162  {
163  MessageOut() << "Opening raw flow output: " << raw_output_file_path << "\n";
164  try {
165  raw_output_file_path.open_stream(raw_eq_data_->raw_output_file);
166  } INPUT_CATCH(FilePath::ExcFileOpen, FilePath::EI_Address_String, (*in_rec_specific))
167 
169  }
170  }
171  }
172 
173  auto fields_array = in_rec_specific->val<Input::Array>("fields");
174  if(fields_array.size() > 0){
176  prepare_specific_output(*in_rec_specific);
177  }
178  }
179 }
180 
182 {
183  // we need to add data from the flow equation at this point, not in constructor of OutputFields
185 
186  // read optional user fields
187  // TODO: check of DarcyLMH type is temporary, remove condition at the same time as removal DarcyMH
188  if(DarcyLMH* d = dynamic_cast<DarcyLMH*>(darcy_flow)) {
189  Input::Array user_fields_arr;
190  if (main_mh_in_rec.opt_val("user_fields", user_fields_arr)) {
191  d->init_user_fields(user_fields_arr, this->output_fields);
192  }
193  }
194 
196 
199 
200  Input::Record in_rec_output = main_mh_in_rec.val<Input::Record>("output");
201  //output_stream->add_admissible_field_names(in_rec_output.val<Input::Array>("fields"));
202  //output_stream->mark_output_times(darcy_flow->time());
204 }
205 
207 {
208  diff_eq_data_ = std::make_shared<DiffEqData>();
209  diff_eq_data_->flow_data_ = darcy_flow->eq_data_;
210 
211  { // init DOF handlers represents element DOFs
212  uint p_elem_component = 1;
213  diff_eq_data_->dh_ = std::make_shared<SubDOFHandlerMultiDim>(diff_eq_data_->flow_data_->dh_, p_elem_component);
214  }
215 
216  // mask 2d elements crossing 1d
217  if (diff_eq_data_->flow_data_->mortar_method_ != DarcyLMH::NoMortar) {
218  diff_eq_data_->velocity_mask.resize(mesh_->n_elements(),0);
220  diff_eq_data_->velocity_mask[ isec.bulk_ele_idx() ]++;
221  }
222  }
223 
225 
226  diff_eq_data_->vel_diff_ptr = create_field_fe<3, FieldValue<3>::Scalar>(diff_eq_data_->dh_);
228  diff_eq_data_->pressure_diff_ptr = create_field_fe<3, FieldValue<3>::Scalar>(diff_eq_data_->dh_);
229  output_specific_fields.pressure_diff.set(diff_eq_data_->pressure_diff_ptr, 0);
230  diff_eq_data_->div_diff_ptr = create_field_fe<3, FieldValue<3>::Scalar>(diff_eq_data_->dh_);
232 
235 
236  if (compute_errors_) {
239  }
240 }
241 
242 /*
243  * Input string of specific output python fields.
244  * Array of 3-items for 1D, 2D, 3D
245  */
246 string spec_fields_input = R"YAML(
247  - source_file: analytical_module.py
248  class: AllValues1D
249  used_fields: ["X"]
250  - source_file: analytical_module.py
251  class: AllValues2D
252  used_fields: ["X"]
253  - source_file: analytical_module.py
254  class: AllValues3D
255  used_fields: ["X"]
256 )YAML";
257 
258 
260 {
261  typedef FieldValue<3>::Scalar ScalarSolution;
262  typedef FieldValue<3>::VectorFixed VectorSolution;
263 
265 
266  // Create separate vectors of 1D, 2D, 3D regions
268  for (unsigned int i=1; i<2*mesh_->region_db().bulk_size(); i+=2) {
269  unsigned int dim = mesh_->region_db().get_dim(i);
270  ASSERT_GT(dim, 0).error("Bulk region with dim==0!\n");
271  reg_by_dim[dim-1].push_back( mesh_->region_db().get_label(i) );
272  }
273 
277  in_arr.copy_to(in_recs);
278 
279  // Create instances of FieldPython and set them to Field objects
280  std::string source_file = "analytical_module.py";
281  for (uint i_dim=0; i_dim<3; ++i_dim) {
282  this->set_ref_solution<ScalarSolution>(in_recs[i_dim], flow_eq_fields_->ref_pressure, reg_by_dim[i_dim]);
283  this->set_ref_solution<VectorSolution>(in_recs[i_dim], flow_eq_fields_->ref_velocity, reg_by_dim[i_dim]);
284  this->set_ref_solution<ScalarSolution>(in_recs[i_dim], flow_eq_fields_->ref_divergence, reg_by_dim[i_dim]);
285  }
286 }
287 
288 
289 template <class FieldType>
291  FieldAlgoBaseInitData init_data(output_field.input_name(), output_field.n_comp(), output_field.units(), output_field.limits(), output_field.flags());
292 
293  std::shared_ptr< FieldPython<3, FieldType> > algo = std::make_shared< FieldPython<3, FieldType> >();
294  algo->init_from_input( in_rec, init_data );
295  output_field.set(algo, darcy_flow->time().t(), reg);
296 }
297 
298 #define DARCY_SET_REF_SOLUTION(FTYPE) \
299 template void DarcyFlowMHOutput::set_ref_solution<FTYPE>(const Input::Record &, Field<3, FTYPE> &, std::vector<std::string>)
300 
303 
304 
306 {
307  if (l2_difference_assembly_ != nullptr) delete l2_difference_assembly_;
309 }
310 
311 
312 
313 
314 
315 //=============================================================================
316 // CONVERT SOLUTION, CALCULATE BALANCES, ETC...
317 //=============================================================================
318 
319 
321 {
322  START_TIMER("Darcy fields output");
323 
324  {
325  START_TIMER("post-process output fields");
326 
327  {
328  if (raw_eq_data_->raw_output_file.is_open())
330  }
331  }
332 
333  {
334  START_TIMER("evaluate output fields");
337  }
338 
339  if (compute_errors_)
340  {
341  START_TIMER("compute specific output fields");
342  //compute_l2_difference();
344  }
345 
347  {
348  START_TIMER("evaluate output fields");
351  }
352 
353 
354 }
355 
#define ASSERT_GT(a, b)
Definition of comparative assert macro (Greater Than) only for debug mode.
Definition: asserts.hh:317
#define ASSERT_PTR(ptr)
Definition of assert macro checking non-null pointer (PTR) only for debug mode.
Definition: asserts.hh:341
Standard quantities for output in DarcyFlowMH.
Field< 3, FieldValue< 3 >::Scalar > subdomain
Field< 3, FieldValue< 3 >::Scalar > region_id
Specific quantities for output in DarcyFlowMH - error estimates etc.
Field< 3, FieldValue< 3 >::Scalar > velocity_diff
Field< 3, FieldValue< 3 >::Scalar > div_diff
Field< 3, FieldValue< 3 >::Scalar > pressure_diff
static const Input::Type::Instance & get_input_type_specific()
std::shared_ptr< DiffEqData > diff_eq_data_
DarcyFlowMHOutput(DarcyLMH *flow, Input::Record in_rec)
GenericAssembly< L2DifferenceAssembly > * l2_difference_assembly_
general assembly objects, hold assembly objects of appropriate dimension
void output()
Calculate values for output.
virtual void prepare_output(Input::Record in_rec)
std::shared_ptr< DarcyLMH::EqFields > flow_eq_fields_
void set_ref_solution(const Input::Record &rec, Field< 3, FieldType > &output_field, std::vector< std::string > reg)
OutputSpecificFields output_specific_fields
bool is_output_specific_fields
Output specific field stuff.
std::shared_ptr< RawOutputEqData > raw_eq_data_
static const Input::Type::Instance & get_input_type(FieldSet &eq_data, const std::string &equation_name)
void prepare_specific_output(Input::Record in_rec)
bool compute_errors_
Specific experimental error computing.
std::shared_ptr< OutputTime > output_stream
GenericAssembly< OutputInternalFlowAssembly > * output_internal_assembly_
Mixed-hybrid model of linear Darcy flow, possibly unsteady.
std::shared_ptr< EqData > eq_data_
std::shared_ptr< EqFields > eq_fields_
FieldSet & eq_fieldset()
Definition: equation.hh:206
TimeGovernor & time()
Definition: equation.hh:151
void initialize(std::shared_ptr< OutputTime > stream, Mesh *mesh, Input::Record in_rec, const TimeGovernor &tg)
const Input::Type::Instance & make_output_type_from_record(Input::Type::Record &in_rec, const string &equation_name, const string &aditional_description="")
void output(TimeStep step)
const Input::Type::Instance & make_output_type(const string &equation_name, const string &aditional_description="")
static Input::Type::Record & get_input_type()
const std::string & input_name() const
FieldCommon & description(const string &description)
FieldCommon & flags(FieldFlag::Flags::Mask mask)
FieldCommon & name(const string &name)
std::pair< double, double > limits() const
FieldCommon & units(const UnitSI &units)
Set basic units of the field.
unsigned int n_comp() const
static constexpr Mask equation_result
Match result fields. These are never given by input or copy of input.
Definition: field_flag.hh:55
static constexpr Mask equation_external_output
Match an output field, that can be also copy of other field.
Definition: field_flag.hh:58
Container for various descendants of FieldCommonBase.
Definition: field_set.hh:159
void set_mesh(const Mesh &mesh)
Definition: field_set.hh:297
bool set_time(const TimeStep &time, LimitSide limit_side)
Definition: field_set.cc:190
Class template representing a field with values dependent on: point, element, and region.
Definition: field.hh:92
void set(FieldBasePtr field, double time, std::vector< std::string > region_set_names={"ALL"})
Definition: field.impl.hh:242
Dedicated class for storing path to input and output files.
Definition: file_path.hh:54
void open_stream(Stream &stream) const
Definition: file_path.cc:211
void assemble(std::shared_ptr< DOFHandlerMultiDim > dh) override
General assemble methods.
static auto region_id(Mesh &mesh) -> IndexField
static auto subdomain(Mesh &mesh) -> IndexField
Accessor to input data conforming to declared Array.
Definition: accessors.hh:566
void copy_to(Container &out) const
Reader for (slightly) modified input files.
T get_root_interface() const
Returns the root accessor.
Accessor to the data with type Type::Record.
Definition: accessors.hh:291
bool opt_val(const string &key, Ret &value) const
const Ret val(const string &key) const
Iterator< Ret > find(const string &key) const
Class for declaration of inputs sequences.
Definition: type_base.hh:339
Class for declaration of the input of type Bool.
Definition: type_base.hh:452
Class Input::Type::Default specifies default value of keys of a Input::Type::Record.
Definition: type_record.hh:61
static Default optional()
The factory function to make an empty default value which is optional.
Definition: type_record.hh:124
static FileName output()
The factory function for declaring type FileName for input files.
Definition: type_base.cc:531
Helper class that stores data of generic types.
Definition: type_generic.hh:89
Record type proxy class.
Definition: type_record.hh:182
Record & close() const
Close the Record for further declarations of keys.
Definition: type_record.cc:304
Record & copy_keys(const Record &other)
Copy keys from other record.
Definition: type_record.cc:216
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:503
Class represents intersection of two elements.
const RegionDB & region_db() const
Definition: mesh.h:175
unsigned int n_elements() const
Definition: mesh.h:111
MixedMeshIntersections & mixed_intersections()
Definition: mesh.cc:849
std::vector< IntersectionLocal< 1, 2 > > intersection_storage12_
Stores 1D-2D intersections.
static std::shared_ptr< OutputTime > create_output_stream(const std::string &equation_name, const Input::Record &in_rec, const std::shared_ptr< TimeUnitConversion > &time_unit_conv)
This method delete all object instances of class OutputTime stored in output_streams vector.
Definition: output_time.cc:187
unsigned int get_dim(unsigned int idx) const
Definition: region.cc:219
const std::string & get_label(unsigned int idx) const
Definition: region.cc:203
std::shared_ptr< TimeUnitConversion > get_unit_conversion() const
Getter for time unit conversion object.
double t() const
const TimeStep & step(int index=-1) const
Class for representation SI units of Fields.
Definition: unit_si.hh:40
static UnitSI & dimensionless()
Returns dimensionless unit.
Definition: unit_si.cc:55
Lumped mixed-hybrid model of linear Darcy flow, possibly unsteady.
string spec_fields_input
#define DARCY_SET_REF_SOLUTION(FTYPE)
Output class for darcy_flow_mh model.
Declaration of class which handles the ordering of degrees of freedom (dof) and mappings between loca...
Global macros to enhance readability and debugging, general constants.
#define INPUT_CATCH(ExceptionType, AddressEITag, input_accessor)
Definition: accessors.hh:63
Classes with algorithms for computation of intersections of meshes.
#define WarningOut()
Macro defining 'warning' record of log.
Definition: logger.hh:278
#define MessageOut()
Macro defining 'message' record of log.
Definition: logger.hh:275
unsigned int uint
#define MPI_COMM_WORLD
Definition: mpi.h:123
#define MPI_Comm_size
Definition: mpi.h:235
#define MPI_Comm_rank
Definition: mpi.h:236
Implementation of range helper class.
Helper struct stores data for initizalize descentants of FieldAlgorithmBase.
#define START_TIMER(tag)
Starts a timer with specified tag.