OGS
ProcessLib::ComponentTransport Namespace Reference

Classes

struct  IntegrationPointData
 
class  ComponentTransportLocalAssemblerInterface
 
class  LocalAssemblerData
 
class  ComponentTransportProcess
 
struct  ComponentTransportProcessData
 
struct  InterpolationPoint
 
struct  Field
 
struct  LookupTable
 

Functions

void checkMPLProperties (MeshLib::Mesh const &mesh, MaterialPropertyLib::MaterialSpatialDistributionMap const &media_map)
 
std::unique_ptr< ProcesscreateComponentTransportProcess (std::string name, MeshLib::Mesh &mesh, std::unique_ptr< ProcessLib::AbstractJacobianAssembler > &&jacobian_assembler, std::vector< ProcessVariable > const &variables, std::vector< std::unique_ptr< ParameterLib::ParameterBase >> const &parameters, unsigned const integration_order, BaseLib::ConfigTree const &config, std::vector< std::unique_ptr< MeshLib::Mesh >> const &meshes, std::map< int, std::shared_ptr< MaterialPropertyLib::Medium >> const &media, std::unique_ptr< ChemistryLib::ChemicalSolverInterface > &&chemical_solver_interface)
 
std::unique_ptr< LookupTablecreateLookupTable (std::optional< std::string > const tabular_file, std::vector< std::vector< std::reference_wrapper< ProcessVariable >>> const &process_variables)
 
static std::vector< std::size_t > intersection (std::vector< std::size_t > const &vec1, std::vector< std::size_t > const &vec2)
 

Function Documentation

◆ checkMPLProperties()

void ProcessLib::ComponentTransport::checkMPLProperties ( MeshLib::Mesh const &  mesh,
MaterialPropertyLib::MaterialSpatialDistributionMap const &  media_map 
)

Definition at line 28 of file CreateComponentTransportProcess.cpp.

31 {
32  std::array const required_properties_medium = {
37 
38  std::array const required_properties_liquid_phase = {
41 
42  std::array const required_properties_components = {
46 
47  for (auto const& element : mesh.getElements())
48  {
49  auto const element_id = element->getID();
50 
51  auto const& medium = *media_map.getMedium(element_id);
52  checkRequiredProperties(medium, required_properties_medium);
53 
54  // check if liquid phase definition and the corresponding properties
55  // exist
56  auto const& liquid_phase = medium.phase("AqueousLiquid");
57  checkRequiredProperties(liquid_phase, required_properties_liquid_phase);
58 
59  // check if components and the corresponding properties exist
60  auto const number_of_components = liquid_phase.numberOfComponents();
61  for (std::size_t component_id = 0; component_id < number_of_components;
62  ++component_id)
63  {
64  if (!liquid_phase.hasComponent(component_id))
65  {
66  OGS_FATAL(
67  "The component {:d} in the AqueousLiquid phase isn't "
68  "specified.",
69  component_id);
70  }
71  auto const& component = liquid_phase.component(component_id);
72  checkRequiredProperties(component, required_properties_components);
73  }
74  }
75 }
#define OGS_FATAL(...)
Definition: Error.h:26
@ longitudinal_dispersivity
used to compute the hydrodynamic dispersion tensor.
Definition: PropertyType.h:58
@ transversal_dispersivity
used to compute the hydrodynamic dispersion tensor.
Definition: PropertyType.h:100
@ retardation_factor
specify retardation factor used in component transport process.
Definition: PropertyType.h:81
void checkRequiredProperties(Component const &c, Container const &required_properties)
Definition: Component.h:96

References MaterialPropertyLib::checkRequiredProperties(), MaterialPropertyLib::decay_rate, MaterialPropertyLib::density, MeshLib::Mesh::getElements(), MaterialPropertyLib::MaterialSpatialDistributionMap::getMedium(), MaterialPropertyLib::longitudinal_dispersivity, OGS_FATAL, MaterialPropertyLib::permeability, MaterialPropertyLib::pore_diffusion, MaterialPropertyLib::porosity, MaterialPropertyLib::retardation_factor, MaterialPropertyLib::transversal_dispersivity, and MaterialPropertyLib::viscosity.

Referenced by createComponentTransportProcess().

◆ createComponentTransportProcess()

std::unique_ptr< Process > ProcessLib::ComponentTransport::createComponentTransportProcess ( std::string  name,
MeshLib::Mesh mesh,
std::unique_ptr< ProcessLib::AbstractJacobianAssembler > &&  jacobian_assembler,
std::vector< ProcessVariable > const &  variables,
std::vector< std::unique_ptr< ParameterLib::ParameterBase >> const &  parameters,
unsigned const  integration_order,
BaseLib::ConfigTree const &  config,
std::vector< std::unique_ptr< MeshLib::Mesh >> const &  meshes,
std::map< int, std::shared_ptr< MaterialPropertyLib::Medium >> const &  media,
std::unique_ptr< ChemistryLib::ChemicalSolverInterface > &&  chemical_solver_interface 
)
Input File Parameter:
prj__processes__process__type
Input File Parameter:
prj__processes__process__ComponentTransport__coupling_scheme
Input File Parameter:
prj__processes__process__ComponentTransport__process_variables
Input File Parameter:
prj__processes__process__ComponentTransport__process_variables__pressure
Input File Parameter:
prj__processes__process__ComponentTransport__process_variables__concentration
Input File Parameter:
prj__processes__process__ComponentTransport__specific_body_force
Input File Parameter:
prj__processes__process__ComponentTransport__non_advective_form
Input File Parameter:
prj__processes__process__ComponentTransport__chemically_induced_porosity_change
Input File Parameter:
prj__processes__process__ComponentTransport__temperature_field
Input File Parameter:
prj__processes__process__ComponentTransport__tabular_file
Input File Parameter:
prj__processes__process__calculatesurfaceflux

Definition at line 77 of file CreateComponentTransportProcess.cpp.

89 {
91  config.checkConfigParameter("type", "ComponentTransport");
92 
93  DBUG("Create ComponentTransportProcess.");
94 
95  auto const coupling_scheme =
97  config.getConfigParameter<std::string>("coupling_scheme",
98  "monolithic_scheme");
99  const bool use_monolithic_scheme = (coupling_scheme != "staggered");
100 
101  // Process variable.
102 
104  auto const pv_config = config.getConfigSubtree("process_variables");
105 
106  std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>
107  process_variables;
108 
109  // Collect all process variables in a vector before allocation
110  // pressure first, concentration then
111  auto const collected_process_variables = findProcessVariables(
112  variables, pv_config,
113  {
114  "pressure",
116  "concentration"});
117 
118  // Check number of components for each process variable
119  auto it = std::find_if(
120  collected_process_variables.cbegin(),
121  collected_process_variables.cend(),
122  [](std::reference_wrapper<ProcessLib::ProcessVariable> const& pv)
123  { return pv.get().getNumberOfGlobalComponents() != 1; });
124 
125  if (it != collected_process_variables.end())
126  {
127  OGS_FATAL(
128  "Number of components for process variable '{:s}' should be 1 "
129  "rather "
130  "than {:d}.",
131  it->get().getName(),
132  it->get().getNumberOfGlobalComponents());
133  }
134  int const hydraulic_process_id = 0;
135  int const first_transport_process_id = use_monolithic_scheme ? 0 : 1;
136 
137  // Allocate the collected process variables into a two-dimensional vector,
138  // depending on what scheme is adopted
139  if (use_monolithic_scheme) // monolithic scheme.
140  {
141  process_variables.push_back(std::move(collected_process_variables));
142  }
143  else // staggered scheme.
144  {
145  std::vector<std::reference_wrapper<ProcessLib::ProcessVariable>>
146  per_process_variable;
147 
148  if (!chemical_solver_interface)
149  {
150  for (auto& pv : collected_process_variables)
151  {
152  per_process_variable.emplace_back(pv);
153  process_variables.push_back(std::move(per_process_variable));
154  }
155  }
156  else
157  {
158  auto sort_by_component =
159  [&per_process_variable,
160  collected_process_variables](auto const& c_name)
161  {
162  auto pv = std::find_if(collected_process_variables.begin(),
163  collected_process_variables.end(),
164  [&c_name](auto const& v) -> bool
165  { return v.get().getName() == c_name; });
166 
167  if (pv == collected_process_variables.end())
168  {
169  OGS_FATAL(
170  "Component {:s} given in "
171  "<chemical_system>/<solution>/"
172  "<components> is not found in specified "
173  "coupled processes (see "
174  "<process>/<process_variables>/"
175  "<concentration>).",
176  c_name);
177  }
178 
179  per_process_variable.emplace_back(*pv);
180  return std::move(per_process_variable);
181  };
182 
183  auto const components =
184  chemical_solver_interface->getComponentList();
185  // pressure
186  per_process_variable.emplace_back(collected_process_variables[0]);
187  process_variables.push_back(std::move(per_process_variable));
188  // concentration
189  assert(components.size() + 1 == collected_process_variables.size());
190  std::transform(components.begin(), components.end(),
191  std::back_inserter(process_variables),
192  sort_by_component);
193  }
194  }
195 
196  // Specific body force parameter.
197  Eigen::VectorXd specific_body_force;
198  std::vector<double> const b =
200  config.getConfigParameter<std::vector<double>>("specific_body_force");
201  assert(!b.empty() && b.size() < 4);
202  if (b.size() < mesh.getDimension())
203  {
204  OGS_FATAL(
205  "specific body force (gravity vector) has {:d} components, mesh "
206  "dimension is {:d}",
207  b.size(), mesh.getDimension());
208  }
209  bool const has_gravity = MathLib::toVector(b).norm() > 0;
210  if (has_gravity)
211  {
212  specific_body_force.resize(b.size());
213  std::copy_n(b.data(), b.size(), specific_body_force.data());
214  }
215 
216  bool const non_advective_form =
218  config.getConfigParameter<bool>("non_advective_form", false);
219 
220  bool chemically_induced_porosity_change =
222  config.getConfigParameter<bool>("chemically_induced_porosity_change",
223  false);
224 
225  auto const temperature = ParameterLib::findOptionalTagParameter<double>(
227  config, "temperature_field", parameters, 1, &mesh);
228 
229  auto media_map =
231 
232  auto lookup_table = ComponentTransport::createLookupTable(
234  config.getConfigParameterOptional<std::string>("tabular_file"),
235  process_variables);
236 
237  DBUG("Check the media properties of ComponentTransport process ...");
238  checkMPLProperties(mesh, *media_map);
239  DBUG("Media properties verified.");
240 
241  ComponentTransportProcessData process_data{
242  std::move(media_map),
243  specific_body_force,
244  has_gravity,
245  non_advective_form,
246  temperature,
247  chemically_induced_porosity_change,
248  chemical_solver_interface.get(),
249  std::move(lookup_table),
250  hydraulic_process_id,
251  first_transport_process_id};
252 
253  SecondaryVariableCollection secondary_variables;
254 
255  ProcessLib::createSecondaryVariables(config, secondary_variables);
256 
257  std::unique_ptr<ProcessLib::SurfaceFluxData> surfaceflux;
258  auto surfaceflux_config =
260  config.getConfigSubtreeOptional("calculatesurfaceflux");
261  if (surfaceflux_config)
262  {
264  *surfaceflux_config, meshes);
265  }
266 
267  return std::make_unique<ComponentTransportProcess>(
268  std::move(name), mesh, std::move(jacobian_assembler), parameters,
269  integration_order, std::move(process_variables),
270  std::move(process_data), std::move(secondary_variables),
271  use_monolithic_scheme, std::move(surfaceflux),
272  std::move(chemical_solver_interface));
273 }
void DBUG(char const *fmt, Args const &... args)
Definition: Logging.h:27
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
Definition: Mesh.h:71
std::unique_ptr< MaterialSpatialDistributionMap > createMaterialSpatialDistributionMap(std::map< int, std::shared_ptr< Medium >> const &media, MeshLib::Mesh const &mesh)
Eigen::Map< const Vector > toVector(std::vector< double > const &data, Eigen::VectorXd::Index size)
Creates an Eigen mapped vector from the given data vector.
std::unique_ptr< LookupTable > createLookupTable(std::optional< std::string > const tabular_file, std::vector< std::vector< std::reference_wrapper< ProcessVariable >>> const &process_variables)
void checkMPLProperties(MeshLib::Mesh const &mesh, MaterialPropertyLib::MaterialSpatialDistributionMap const &media_map)
std::vector< std::reference_wrapper< ProcessVariable > > findProcessVariables(std::vector< ProcessVariable > const &variables, BaseLib::ConfigTree const &pv_config, std::initializer_list< std::string > tags)
void createSecondaryVariables(BaseLib::ConfigTree const &config, SecondaryVariableCollection &secondary_variables)
static std::unique_ptr< ProcessLib::SurfaceFluxData > createSurfaceFluxData(BaseLib::ConfigTree const &calculatesurfaceflux_config, std::vector< std::unique_ptr< MeshLib::Mesh >> const &meshes)

References BaseLib::ConfigTree::checkConfigParameter(), checkMPLProperties(), createLookupTable(), MaterialPropertyLib::createMaterialSpatialDistributionMap(), ProcessLib::createSecondaryVariables(), ProcessLib::SurfaceFluxData::createSurfaceFluxData(), DBUG(), ProcessLib::findProcessVariables(), BaseLib::ConfigTree::getConfigParameter(), BaseLib::ConfigTree::getConfigParameterOptional(), BaseLib::ConfigTree::getConfigSubtree(), BaseLib::ConfigTree::getConfigSubtreeOptional(), MeshLib::Mesh::getDimension(), MaterialPropertyLib::name, OGS_FATAL, and MathLib::toVector().

Referenced by ProjectData::parseProcesses().

◆ createLookupTable()

std::unique_ptr< LookupTable > ProcessLib::ComponentTransport::createLookupTable ( std::optional< std::string > const  tabular_file,
std::vector< std::vector< std::reference_wrapper< ProcessVariable >>> const &  process_variables 
)

Definition at line 55 of file CreateLookupTable.cpp.

59 {
60  if (!tabular_file)
61  {
62  return nullptr;
63  }
64 
65  auto const path_to_tabular_file =
67 
68  if (!BaseLib::IsFileExisting(path_to_tabular_file))
69  {
70  OGS_FATAL(
71  "Not found the tabular file with the specified file path: {:s}",
72  path_to_tabular_file);
73  }
74 
75  INFO("Found the tabular file: {:s}", path_to_tabular_file);
76 
77  std::ifstream in(path_to_tabular_file);
78  if (!in)
79  {
80  OGS_FATAL("Couldn't open the tabular file: {:s}.",
81  path_to_tabular_file);
82  }
83 
84  // read field names
85  std::string line;
86  std::getline(in, line);
87  std::vector<std::string> field_names;
88  boost::split(field_names, line, boost::is_any_of("\t "));
89 
90  // categorize entry fields
91  std::vector<std::string> input_field_names;
92  std::copy_if(field_names.begin(), field_names.end(),
93  std::back_inserter(input_field_names),
94  [](auto const& field_name) -> bool
95  { return field_name.find("_new") == std::string::npos; });
96 
97  // read table entries
98  std::map<std::string, std::vector<double>> tabular_data;
99  while (std::getline(in, line))
100  {
101  std::vector<std::string> field_data;
102  boost::split(field_data, line, boost::is_any_of("\t "));
103 
104  assert(field_data.size() == field_names.size());
105  for (std::size_t field_id = 0; field_id < field_data.size(); ++field_id)
106  {
107  tabular_data[field_names[field_id]].push_back(
108  std::stod(field_data[field_id]));
109  }
110  }
111  in.close();
112 
113  std::vector<Field> input_fields;
114  input_fields.reserve(input_field_names.size());
115  for (auto const& field_name : input_field_names)
116  {
117  auto pv = std::find_if(process_variables.begin(),
118  process_variables.end(),
119  [&field_name](auto const& v) -> bool
120  {
121  return v[0].get().getName() == field_name ||
122  v[0].get().getName() + "_prev" ==
123  field_name;
124  });
125 
126  if (pv == process_variables.end())
127  {
128  OGS_FATAL(
129  "Not found field name {:s} in the group of process variables "
130  "defined in the project file.",
131  field_name);
132  }
133 
134  auto const process_id =
135  static_cast<int>(std::distance(process_variables.cbegin(), pv));
136 
137  auto seed_points = tabular_data[field_name];
138  BaseLib::makeVectorUnique(seed_points);
139 
140  std::vector<std::vector<std::size_t>> point_id_groups;
141  for (auto const seed_point : seed_points)
142  {
143  auto const point_id_group =
144  getIndexVector(tabular_data[field_name], seed_point);
145  point_id_groups.push_back(point_id_group);
146  }
147 
148  input_fields.emplace_back(point_id_groups, seed_points, field_name,
149  process_id);
150  }
151 
152  return std::make_unique<LookupTable>(std::move(input_fields),
153  std::move(tabular_data));
154 }
void INFO(char const *fmt, Args const &... args)
Definition: Logging.h:32
std::string const & getProjectDirectory()
Returns the directory where the prj file resides.
Definition: FileTools.cpp:217
bool IsFileExisting(const std::string &strFilename)
Returns true if given file exists.
Definition: FileTools.cpp:43
std::string joinPaths(std::string const &pathA, std::string const &pathB)
Definition: FileTools.cpp:212
void makeVectorUnique(std::vector< T > &v)
Definition: Algorithm.h:209
std::vector< std::size_t > getIndexVector(std::vector< double > const &data, double const value)

References anonymous_namespace{CreateLookupTable.cpp}::getIndexVector(), BaseLib::getProjectDirectory(), INFO(), BaseLib::IsFileExisting(), BaseLib::joinPaths(), BaseLib::makeVectorUnique(), and OGS_FATAL.

Referenced by createComponentTransportProcess().

◆ intersection()

static std::vector<std::size_t> ProcessLib::ComponentTransport::intersection ( std::vector< std::size_t > const &  vec1,
std::vector< std::size_t > const &  vec2 
)
static

Definition at line 21 of file LookupTable.cpp.

23 {
24  std::unordered_set<std::size_t> set(vec1.begin(), vec1.end());
25  std::vector<std::size_t> vec;
26  for (auto const a : vec2)
27  {
28  if (set.contains(a))
29  {
30  vec.push_back(a);
31  set.erase(a);
32  }
33  }
34  return vec;
35 }
void set(PETScVector &x, double const a)
Definition: LinAlg.cpp:32

References MathLib::LinAlg::set().

Referenced by MeshGeoToolsLib::GeoMapper::getMeshElevation(), and ProcessLib::ComponentTransport::LookupTable::getTableEntryID().