OGS
ProcessLib::ComponentTransport Namespace Reference

Classes

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

Functions

void checkMPLProperties (MeshLib::Mesh const &mesh, MaterialPropertyLib::MaterialSpatialDistributionMap const &media_map)
 
std::unique_ptr< ProcesscreateComponentTransportProcess (std::string const &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 void intersection (std::vector< std::size_t > &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 31 of file CreateComponentTransportProcess.cpp.

34{
35 std::array const required_properties_medium = {
40
41 std::array const required_properties_liquid_phase = {
44
45 std::array const required_properties_components = {
49
50 for (auto const& element_id : mesh.getElements() | MeshLib::views::ids)
51 {
52 auto const& medium = *media_map.getMedium(element_id);
53 checkRequiredProperties(medium, required_properties_medium);
54
55 // check if liquid phase definition and the corresponding properties
56 // exist
57 auto const& liquid_phase = medium.phase("AqueousLiquid");
58 checkRequiredProperties(liquid_phase, required_properties_liquid_phase);
59
60 // check if components and the corresponding properties exist
61 auto const number_of_components = liquid_phase.numberOfComponents();
62 for (std::size_t component_id = 0; component_id < number_of_components;
63 ++component_id)
64 {
65 if (!liquid_phase.hasComponent(component_id))
66 {
68 "The component {:d} in the AqueousLiquid phase isn't "
69 "specified.",
70 component_id);
71 }
72 auto const& component = liquid_phase.component(component_id);
73 checkRequiredProperties(component, required_properties_components);
74 }
75 }
76}
#define OGS_FATAL(...)
Definition Error.h:26
void checkRequiredProperties(Component const &c, std::span< PropertyType const > const required_properties)
Definition Component.cpp:60
@ longitudinal_dispersivity
used to compute the hydrodynamic dispersion tensor.
@ transversal_dispersivity
used to compute the hydrodynamic dispersion tensor.
@ retardation_factor
specify retardation factor used in component transport process.
constexpr ranges::views::view_closure ids
For an element of a range view return its id.
Definition Mesh.h:225

References MaterialPropertyLib::decay_rate, MaterialPropertyLib::density, MeshLib::Mesh::getElements(), MaterialPropertyLib::MaterialSpatialDistributionMap::getMedium(), MeshLib::views::ids, 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 const & 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

Process Variables

Input File Parameter
prj__processes__process__ComponentTransport__process_variables

Primary process variables as they appear in the global component vector.

Input File Parameter
prj__processes__process__ComponentTransport__process_variables__pressure
Input File Parameter
prj__processes__process__ComponentTransport__process_variables__concentration

Temperature is optional.

Input File Parameter
prj__processes__process__ComponentTransport__process_variables__temperature

Process Parameters

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__ComponentTransport__aperture_size
Input File Parameter
prj__processes__process__ComponentTransport__aperture_size__parameter
Input File Parameter
prj__processes__process__ComponentTransport__linear
Input File Parameter
prj__processes__process__ComponentTransport__linear_solver_compute_only_upon_timestep_change
Input File Parameter
prj__processes__process__calculatesurfaceflux

Definition at line 78 of file CreateComponentTransportProcess.cpp.

90{
92 config.checkConfigParameter("type", "ComponentTransport");
93
94 DBUG("Create ComponentTransportProcess.");
95
96 auto const coupling_scheme =
98 config.getConfigParameter<std::string>("coupling_scheme",
99 "monolithic_scheme");
100 const bool use_monolithic_scheme = (coupling_scheme != "staggered");
101
103
105 auto const pv_config = config.getConfigSubtree("process_variables");
106
107 std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>
108 process_variables;
109
111 std::vector collected_process_variables = findProcessVariables(
112 variables, pv_config,
113 {
114 "pressure",
116 "concentration"});
117
119 std::vector const temperature_variable = findProcessVariables(
120 variables, pv_config,
122 "temperature", true /*temperature is optional*/);
123 bool const isothermal = temperature_variable.empty();
124 if (!isothermal)
125 {
126 assert(temperature_variable.size() == 1);
127 collected_process_variables.insert(
128 ++collected_process_variables.begin(), temperature_variable[0]);
129 }
130
131 // Check number of components for each process variable
132 auto it = std::find_if(
133 collected_process_variables.cbegin(),
134 collected_process_variables.cend(),
135 [](std::reference_wrapper<ProcessLib::ProcessVariable> const& pv)
136 { return pv.get().getNumberOfGlobalComponents() != 1; });
137
138 if (it != collected_process_variables.end())
139 {
140 OGS_FATAL(
141 "Number of components for process variable '{:s}' should be 1 "
142 "rather "
143 "than {:d}.",
144 it->get().getName(),
145 it->get().getNumberOfGlobalComponents());
146 }
147
148 // Allocate the collected process variables into a two-dimensional vector,
149 // depending on what scheme is adopted
150 if (use_monolithic_scheme) // monolithic scheme.
151 {
152 if (!isothermal)
153 {
154 OGS_FATAL(
155 "Currently, non-isothermal component transport process can "
156 "only be simulated in staggered scheme.");
157 }
158
159 process_variables.push_back(std::move(collected_process_variables));
160 }
161 else // staggered scheme.
162 {
163 std::vector<std::reference_wrapper<ProcessLib::ProcessVariable>>
164 per_process_variable;
165
166 if (!chemical_solver_interface)
167 {
168 for (auto& pv : collected_process_variables)
169 {
170 per_process_variable.emplace_back(pv);
171 process_variables.push_back(std::move(per_process_variable));
172 }
173 }
174 else
175 {
176 auto sort_by_component =
177 [&per_process_variable,
178 collected_process_variables](auto const& c_name)
179 {
180 auto pv = std::find_if(collected_process_variables.begin(),
181 collected_process_variables.end(),
182 [&c_name](auto const& v) -> bool
183 { return v.get().getName() == c_name; });
184
185 if (pv == collected_process_variables.end())
186 {
187 OGS_FATAL(
188 "Component {:s} given in "
189 "<chemical_system>/<solution>/"
190 "<components> is not found in specified "
191 "coupled processes (see "
192 "<process>/<process_variables>/"
193 "<concentration>).",
194 c_name);
195 }
196
197 per_process_variable.emplace_back(*pv);
198 return std::move(per_process_variable);
199 };
200
201 auto const components =
202 chemical_solver_interface->getComponentList();
203 // pressure
204 per_process_variable.emplace_back(collected_process_variables[0]);
205 process_variables.push_back(std::move(per_process_variable));
206 // temperature
207 if (!isothermal)
208 {
209 per_process_variable.emplace_back(
210 collected_process_variables[1]);
211 process_variables.push_back(std::move(per_process_variable));
212 }
213 // concentration
214 assert(components.size() + (isothermal ? 1 : 2) ==
215 collected_process_variables.size());
216 std::transform(components.begin(), components.end(),
217 std::back_inserter(process_variables),
218 sort_by_component);
219 }
220 }
221
223 std::vector<double> const b =
225 config.getConfigParameter<std::vector<double>>("specific_body_force");
226 assert(!b.empty() && b.size() < 4);
227 // Specific body force parameter.
228 Eigen::VectorXd specific_body_force(b.size());
229 int const mesh_space_dimension =
231 if (static_cast<int>(b.size()) != mesh_space_dimension)
232 {
233 OGS_FATAL(
234 "specific body force (gravity vector) has {:d} components, mesh "
235 "dimension is {:d}",
236 b.size(), mesh_space_dimension);
237 }
238 bool const has_gravity = MathLib::toVector(b).norm() > 0;
239 if (has_gravity)
240 {
241 std::copy_n(b.data(), b.size(), specific_body_force.data());
242 }
243
244 bool const non_advective_form =
246 config.getConfigParameter<bool>("non_advective_form", false);
247
248 bool chemically_induced_porosity_change =
250 config.getConfigParameter<bool>("chemically_induced_porosity_change",
251 false);
252
253 auto const temperature_field = ParameterLib::findOptionalTagParameter<
254 double>(
256 config, "temperature_field", parameters, 1, &mesh);
257 if (!isothermal && temperature_field != nullptr)
258 {
259 OGS_FATAL("Temperature field is set for non-isothermal setup.")
260 }
261
262 auto media_map =
264
265 auto lookup_table = ComponentTransport::createLookupTable(
267 config.getConfigParameterOptional<std::string>("tabular_file"),
268 process_variables);
269
270 DBUG("Check the media properties of ComponentTransport process ...");
271 checkMPLProperties(mesh, media_map);
272 DBUG("Media properties verified.");
273
274 auto stabilizer = NumLib::createNumericalStabilization(mesh, config);
275
276 auto const* aperture_size_parameter = &ParameterLib::findParameter<double>(
278 auto const aperture_config =
280 config.getConfigSubtreeOptional("aperture_size");
281 if (aperture_config)
282 {
283 aperture_size_parameter = &ParameterLib::findParameter<double>(
285 *aperture_config, "parameter", parameters, 1);
286 }
287
288 auto const is_linear =
290 config.getConfigParameter("linear", false);
291
292 auto const ls_compute_only_upon_timestep_change =
294 config.getConfigParameter(
295 "linear_solver_compute_only_upon_timestep_change", false);
296
297 auto const rotation_matrices = MeshLib::getElementRotationMatrices(
298 mesh_space_dimension, mesh.getDimension(), mesh.getElements());
299 std::vector<Eigen::VectorXd> projected_specific_body_force_vectors;
300 projected_specific_body_force_vectors.reserve(rotation_matrices.size());
301
302 std::transform(rotation_matrices.begin(), rotation_matrices.end(),
303 std::back_inserter(projected_specific_body_force_vectors),
304 [&specific_body_force](const auto& R)
305 { return R * R.transpose() * specific_body_force; });
306
307 ComponentTransportProcessData process_data{
308 std::move(media_map),
309 has_gravity,
310 non_advective_form,
311 temperature_field,
312 chemically_induced_porosity_change,
313 chemical_solver_interface.get(),
314 std::move(lookup_table),
315 std::move(stabilizer),
316 projected_specific_body_force_vectors,
317 mesh_space_dimension,
318 *aperture_size_parameter,
319 isothermal,
320 };
321
322 SecondaryVariableCollection secondary_variables;
323
324 ProcessLib::createSecondaryVariables(config, secondary_variables);
325
326 std::unique_ptr<ProcessLib::SurfaceFluxData> surfaceflux;
327 auto surfaceflux_config =
329 config.getConfigSubtreeOptional("calculatesurfaceflux");
330 if (surfaceflux_config)
331 {
333 *surfaceflux_config, meshes);
334 }
335
336 return std::make_unique<ComponentTransportProcess>(
337 std::move(name), mesh, std::move(jacobian_assembler), parameters,
338 integration_order, std::move(process_variables),
339 std::move(process_data), std::move(secondary_variables),
340 use_monolithic_scheme, std::move(surfaceflux),
341 std::move(chemical_solver_interface), is_linear,
342 ls_compute_only_upon_timestep_change);
343}
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:30
std::vector< Node * > const & getNodes() const
Get the nodes-vector for the mesh.
Definition Mesh.h:106
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition Mesh.h:109
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
Definition Mesh.h:88
static PROCESSLIB_EXPORT const std::string constant_one_parameter_name
Definition Process.h:46
MaterialSpatialDistributionMap createMaterialSpatialDistributionMap(std::map< int, std::shared_ptr< Medium > > const &media, MeshLib::Mesh const &mesh)
constexpr int size(int const displacement_dim)
Vectorized tensor size for given displacement dimension.
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::vector< Eigen::MatrixXd > getElementRotationMatrices(int const space_dimension, int const mesh_dimension, std::vector< Element * > const &elements)
Element rotation matrix computation.
int getSpaceDimension(std::vector< Node * > const &nodes)
Computes dimension of the embedding space containing the set of given points.
NumericalStabilization createNumericalStabilization(MeshLib::Mesh const &mesh, BaseLib::ConfigTree const &config)
Parameter< ParameterDataType > * findOptionalTagParameter(BaseLib::ConfigTree const &process_config, std::string const &tag, std::vector< std::unique_ptr< ParameterBase > > const &parameters, int const num_components, MeshLib::Mesh const *const mesh=nullptr)
Definition Utils.h:162
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(), ProcessLib::Process::constant_one_parameter_name, createLookupTable(), MaterialPropertyLib::createMaterialSpatialDistributionMap(), NumLib::createNumericalStabilization(), ProcessLib::createSecondaryVariables(), ProcessLib::SurfaceFluxData::createSurfaceFluxData(), DBUG(), ParameterLib::findOptionalTagParameter(), ProcessLib::findProcessVariables(), BaseLib::ConfigTree::getConfigParameter(), BaseLib::ConfigTree::getConfigParameterOptional(), BaseLib::ConfigTree::getConfigSubtree(), BaseLib::ConfigTree::getConfigSubtreeOptional(), MeshLib::Mesh::getDimension(), MeshLib::getElementRotationMatrices(), MeshLib::Mesh::getElements(), MeshLib::Mesh::getNodes(), MeshLib::getSpaceDimension(), 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 {
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(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:35
std::string const & getProjectDirectory()
Returns the directory where the prj file resides.
bool IsFileExisting(const std::string &strFilename)
Returns true if given file exists.
Definition FileTools.cpp:47
std::string joinPaths(std::string const &pathA, std::string const &pathB)
void makeVectorUnique(std::vector< T > &v)
Definition Algorithm.h:175
std::vector< std::size_t > getIndexVector(std::vector< double > const &data, double const value)

References BaseLib::getProjectDirectory(), INFO(), BaseLib::IsFileExisting(), BaseLib::joinPaths(), BaseLib::makeVectorUnique(), and OGS_FATAL.

Referenced by createComponentTransportProcess().

◆ intersection()

static void ProcessLib::ComponentTransport::intersection ( std::vector< std::size_t > & 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 vec1.clear();
26 for (auto const a : vec2)
27 {
28 if (set.contains(a))
29 {
30 vec1.push_back(a);
31 set.erase(a);
32 }
33 }
34}

Referenced by ProcessLib::ComponentTransport::LookupTable::getTableEntryID().