ASPECT
interface.h
Go to the documentation of this file.
1 /*
2  Copyright (C) 2011 - 2022 by the authors of the ASPECT code.
3 
4  This file is part of ASPECT.
5 
6  ASPECT is free software; you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation; either version 2, or (at your option)
9  any later version.
10 
11  ASPECT is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with ASPECT; see the file LICENSE. If not see
18  <http://www.gnu.org/licenses/>.
19 */
20 
21 #ifndef _aspect_material_model_interface_h
22 #define _aspect_material_model_interface_h
23 
24 #include <aspect/global.h>
25 #include <aspect/plugins.h>
27 
28 #include <deal.II/base/point.h>
29 
30 // Work around an incorrect instantiation in qprojector.h of deal.II 9.2.0,
31 // which requires including qprojector.h before quadrature.h (and not
32 // after). This file doesn't actually need qprojector.h, so the include can be
33 // removed when we require 9.3.. For more info see
34 // https://github.com/geodynamics/aspect/issues/3728
35 #include <deal.II/base/quadrature.h>
36 
37 #include <deal.II/base/symmetric_tensor.h>
38 #include <deal.II/base/parameter_handler.h>
39 #include <deal.II/dofs/dof_handler.h>
40 #include <deal.II/dofs/dof_accessor.h>
41 #include <deal.II/fe/mapping.h>
42 #include <deal.II/fe/fe_values.h>
43 #include <deal.II/fe/component_mask.h>
44 #include <deal.II/numerics/data_postprocessor.h>
45 #include <deal.II/base/signaling_nan.h>
46 
47 namespace aspect
48 {
49  template <int dim>
50  struct Introspection;
51 
52 
60  namespace MaterialModel
61  {
62  using namespace dealii;
63 
68  namespace NonlinearDependence
69  {
92  {
94 
95  none = 1,
97  pressure = 4,
100 
102  };
103 
104 
109  const Dependence d2)
110  {
111  return Dependence(static_cast<int>(d1) | static_cast<int>(d2));
112  }
113 
115  const Dependence d2)
116  {
117  d1 = (d1 | d2);
118  return d1;
119  }
120 
126  {
132 
138 
144 
150 
156 
162  ModelDependence ();
163  };
164 
173  bool
174  identifies_single_variable(const Dependence dependence);
175  }
176 
181  namespace MaterialProperties
182  {
192  enum Property
193  {
195 
196  none = 1,
198  density = 4,
208 
211  specific_heat |
216  viscosity |
221  };
222 
227  inline Property operator | (const Property d1,
228  const Property d2)
229  {
230  return Property(static_cast<int>(d1) | static_cast<int>(d2));
231  }
232 
234  const Property d2)
235  {
236  return (d1 | d2);
237  }
238  }
239 
240  template <int dim> class AdditionalMaterialInputs;
241 
248  template <int dim>
250  {
261  MaterialModelInputs(const unsigned int n_points,
262  const unsigned int n_comp);
263 
279  MaterialModelInputs(const DataPostprocessorInputs::Vector<dim> &input_data,
280  const Introspection<dim> &introspection,
281  const bool compute_strain_rate = true);
282 
283 
303  MaterialModelInputs(const FEValuesBase<dim,dim> &fe_values,
304  const typename DoFHandler<dim>::active_cell_iterator &cell,
305  const Introspection<dim> &introspection,
306  const LinearAlgebra::BlockVector &solution_vector,
307  const bool compute_strain_rate = true);
308 
323 
327  MaterialModelInputs (MaterialModelInputs &&) noexcept = default;
328 
333  MaterialModelInputs &operator= (const MaterialModelInputs &source) = delete;
334 
338  MaterialModelInputs &operator= (MaterialModelInputs &&) = default;
339 
345  void reinit(const FEValuesBase<dim,dim> &fe_values,
346  const typename DoFHandler<dim>::active_cell_iterator &cell,
347  const Introspection<dim> &introspection,
348  const LinearAlgebra::BlockVector &solution_vector,
349  const bool compute_strain_rate = true);
350 
355  unsigned int n_evaluation_points() const;
356 
363  bool requests_property(const MaterialProperties::Property &property) const;
364 
369  std::vector<Point<dim>> position;
370 
374  std::vector<double> temperature;
375 
379  std::vector<double> pressure;
380 
385  std::vector<Tensor<1,dim>> pressure_gradient;
386 
394  std::vector<Tensor<1,dim>> velocity;
395 
401  std::vector<std::vector<double>> composition;
402 
415  std::vector<SymmetricTensor<2,dim>> strain_rate;
416 
432  typename DoFHandler<dim>::active_cell_iterator current_cell;
433 
442 
448  std::vector<std::unique_ptr<AdditionalMaterialInputs<dim>>> additional_inputs;
449 
456  template <class AdditionalInputType>
457  AdditionalInputType *get_additional_input();
458 
462  template <class AdditionalInputType>
463  const AdditionalInputType *get_additional_input() const;
464  };
465 
466 
467  template <int dim> class AdditionalMaterialOutputs;
468 
469 
476  template <int dim>
478  {
489  MaterialModelOutputs (const unsigned int n_points,
490  const unsigned int n_comp);
491 
506 
510  MaterialModelOutputs (MaterialModelOutputs &&) noexcept = default;
511 
516  MaterialModelOutputs &operator= (const MaterialModelOutputs &source) = delete;
517 
521  MaterialModelOutputs &operator= (MaterialModelOutputs &&) noexcept = default;
522 
527  unsigned int n_evaluation_points() const;
528 
532  std::vector<double> viscosities;
533 
537  std::vector<double> densities;
538 
543  std::vector<double> thermal_expansion_coefficients;
544 
548  std::vector<double> specific_heat;
549 
553  std::vector<double> thermal_conductivities;
554 
559  std::vector<double> compressibilities;
560 
566  std::vector<double> entropy_derivative_pressure;
567 
574  std::vector<double> entropy_derivative_temperature;
575 
612  std::vector<std::vector<double>> reaction_terms;
613 
619  std::vector<std::unique_ptr<AdditionalMaterialOutputs<dim>>> additional_outputs;
620 
628  template <class AdditionalOutputType>
629  AdditionalOutputType *get_additional_output();
630 
634  template <class AdditionalOutputType>
635  const AdditionalOutputType *get_additional_output() const;
636 
641  void move_additional_outputs_from(MaterialModelOutputs<dim> &other);
642  };
643 
644 
645 
660  namespace MaterialAveraging
661  {
701  {
712  };
713 
714 
721  std::string get_averaging_operation_names ();
722 
728  AveragingOperation parse_averaging_operation_name (const std::string &s);
729 
735  template <int dim>
736  void average (const AveragingOperation operation,
737  const typename DoFHandler<dim>::active_cell_iterator &cell,
738  const Quadrature<dim> &quadrature_formula,
739  const Mapping<dim> &mapping,
740  MaterialModelOutputs<dim> &values_out);
741 
747  void average_property (const AveragingOperation operation,
748  const FullMatrix<double> &projection_matrix,
749  const FullMatrix<double> &expansion_matrix,
750  std::vector<double> &values_out);
751  }
752 
753 
772  template <int dim>
774  {
775  public:
776  virtual ~AdditionalMaterialInputs() = default;
777 
782  virtual void
783  fill (const LinearAlgebra::BlockVector &solution,
784  const FEValuesBase<dim> &fe_values,
785  const Introspection<dim> &introspection) = 0;
786  };
787 
788 
815  template <int dim>
817  {
818  public:
823  virtual ~AdditionalMaterialOutputs() = default;
824 
825  virtual void average (const MaterialAveraging::AveragingOperation /*operation*/,
826  const FullMatrix<double> &/*projection_matrix*/,
827  const FullMatrix<double> &/*expansion_matrix*/)
828  {}
829  };
830 
831 
853  template <int dim>
855  {
856  public:
865  NamedAdditionalMaterialOutputs(const std::vector<std::string> &output_names);
866 
877  NamedAdditionalMaterialOutputs(const std::vector<std::string> &output_names,
878  const unsigned int n_points);
879 
883  ~NamedAdditionalMaterialOutputs() override;
884 
889  const std::vector<std::string> &get_names() const;
890 
895  virtual std::vector<double> get_nth_output(const unsigned int idx) const;
896 
898  const FullMatrix<double> &/*projection_matrix*/,
899  const FullMatrix<double> &/*expansion_matrix*/) override
900  {}
901 
902 
907  std::vector<std::vector<double>> output_values;
908 
909  private:
910  const std::vector<std::string> names;
911  };
912 
913 
919  template <int dim>
921  {
922  public:
923  SeismicAdditionalOutputs(const unsigned int n_points);
924 
925  std::vector<double> get_nth_output(const unsigned int idx) const override;
926 
932  std::vector<double> vs;
933 
939  std::vector<double> vp;
940  };
941 
942 
963  template <int dim>
965  {
966  public:
967  ReactionRateOutputs (const unsigned int n_points,
968  const unsigned int n_comp);
969 
970  std::vector<double> get_nth_output(const unsigned int idx) const override;
971 
979  std::vector<std::vector<double>> reaction_rates;
980  };
981 
982 
983 
989  template <int dim>
991  {
992  public:
993  PhaseOutputs(const unsigned int n_points);
994  };
995 
996 
997 
1016  template <int dim>
1018  {
1019  public:
1020  PrescribedFieldOutputs (const unsigned int n_points,
1021  const unsigned int n_comp);
1022 
1023  std::vector<double> get_nth_output(const unsigned int idx) const override;
1024 
1032  std::vector<std::vector<double>> prescribed_field_outputs;
1033  };
1034 
1048  template <int dim>
1050  {
1051  public:
1052  PrescribedTemperatureOutputs (const unsigned int n_points);
1053 
1054  std::vector<double> get_nth_output(const unsigned int idx) const override;
1055 
1063  std::vector<double> prescribed_temperature_outputs;
1064  };
1065 
1072  template <int dim>
1074  {
1075  public:
1076  AdditionalMaterialOutputsStokesRHS(const unsigned int n_points)
1077  : rhs_u(n_points), rhs_p(n_points), rhs_melt_pc(n_points)
1078  {}
1079 
1081  = default;
1082 
1084  const FullMatrix<double> &/*projection_matrix*/,
1085  const FullMatrix<double> &/*expansion_matrix*/) override
1086  {
1087  // TODO: not implemented
1088  }
1089 
1095  std::vector<Tensor<1,dim>> rhs_u;
1096 
1101  std::vector<double> rhs_p;
1102 
1107  std::vector<double> rhs_melt_pc;
1108  };
1109 
1110 
1111 
1127  template <int dim>
1129  {
1130  public:
1134  explicit PrescribedPlasticDilation (const unsigned int n_points);
1135 
1139  virtual std::vector<double> get_nth_output(const unsigned int idx) const;
1140 
1145  std::vector<double> dilation;
1146  };
1147 
1148 
1149 
1156  template <int dim>
1158  {
1159  public:
1160  ElasticOutputs(const unsigned int n_points)
1161  : elastic_force(n_points, numbers::signaling_nan<Tensor<2,dim>>() )
1162  {}
1163 
1164  ~ElasticOutputs() override
1165  = default;
1166 
1168  const FullMatrix<double> &/*projection_matrix*/,
1169  const FullMatrix<double> &/*expansion_matrix*/) override
1170  {
1171  AssertThrow(operation == MaterialAveraging::AveragingOperation::none,ExcNotImplemented());
1172  return;
1173  }
1174 
1180  std::vector<Tensor<2,dim>> elastic_force;
1181  };
1182 
1183 
1184 
1195  template <int dim>
1197  {
1198  public:
1199  EnthalpyOutputs(const unsigned int n_points)
1200  : enthalpies_of_fusion(n_points, numbers::signaling_nan<double>())
1201  {}
1202 
1203  virtual ~EnthalpyOutputs()
1204  = default;
1205 
1206  virtual void average (const MaterialAveraging::AveragingOperation operation,
1207  const FullMatrix<double> &/*projection_matrix*/,
1208  const FullMatrix<double> &/*expansion_matrix*/)
1209  {
1210  AssertThrow(operation == MaterialAveraging::AveragingOperation::none,ExcNotImplemented());
1211  return;
1212  }
1213 
1219  std::vector<double> enthalpies_of_fusion;
1220  };
1221 
1222 
1223 
1239  template <int dim>
1240  class Interface
1241  {
1242  public:
1257 
1262  virtual ~Interface() = default;
1263 
1269  virtual
1270  void
1271  initialize ();
1272 
1277  virtual void update ();
1278 
1291  get_model_dependence () const;
1292 
1302  virtual bool is_compressible () const = 0;
1312  virtual
1313  void evaluate (const MaterialModel::MaterialModelInputs<dim> &in,
1325  static
1326  void
1327  declare_parameters (ParameterHandler &prm);
1328 
1335  virtual
1336  void
1337  parse_parameters (ParameterHandler &prm);
1347  virtual
1348  void
1349  create_additional_named_outputs (MaterialModelOutputs &outputs) const;
1350 
1351 
1359  virtual
1360  void
1361  fill_additional_material_model_inputs(MaterialModel::MaterialModelInputs<dim> &input,
1362  const LinearAlgebra::BlockVector &solution,
1363  const FEValuesBase<dim> &fe_values,
1364  const Introspection<dim> &introspection) const;
1365 
1366  protected:
1383  };
1384 
1400  template <int dim>
1401  void
1402  register_material_model (const std::string &name,
1403  const std::string &description,
1404  void (*declare_parameters_function) (ParameterHandler &),
1405  Interface<dim> *(*factory_function) ());
1406 
1417  template <int dim>
1418  Interface<dim> *
1419  create_material_model (const std::string &model_name);
1420 
1421 
1433  template <int dim>
1434  Interface<dim> *
1435  create_material_model (ParameterHandler &prm);
1436 
1437 
1444  template <int dim>
1445  std::string
1447 
1448 
1454  template <int dim>
1455  void
1456  declare_parameters (ParameterHandler &prm);
1457 
1458 
1459 
1469  template <int dim>
1470  void
1471  write_plugin_graph (std::ostream &output_stream);
1472 
1473 
1474 
1475 // --------------------- template function definitions ----------------------------------
1476 
1477  template <int dim>
1478  template <class AdditionalInputType>
1480  {
1481  for (unsigned int i=0; i<additional_inputs.size(); ++i)
1482  {
1483  AdditionalInputType *result = dynamic_cast<AdditionalInputType *> (additional_inputs[i].get());
1484  if (result)
1485  return result;
1486  }
1487  return nullptr;
1488  }
1489 
1490 
1491  template <int dim>
1492  template <class AdditionalInputType>
1493  const AdditionalInputType *MaterialModelInputs<dim>::get_additional_input() const
1494  {
1495  for (unsigned int i=0; i<additional_inputs.size(); ++i)
1496  {
1497  const AdditionalInputType *result = dynamic_cast<const AdditionalInputType *> (additional_inputs[i].get());
1498  if (result)
1499  return result;
1500  }
1501  return nullptr;
1502  }
1503 
1504 
1505  template <int dim>
1506  template <class AdditionalOutputType>
1508  {
1509  for (unsigned int i=0; i<additional_outputs.size(); ++i)
1510  {
1511  AdditionalOutputType *result = dynamic_cast<AdditionalOutputType *> (additional_outputs[i].get());
1512  if (result)
1513  return result;
1514  }
1515  return nullptr;
1516  }
1517 
1518 
1519  template <int dim>
1520  template <class AdditionalOutputType>
1521  const AdditionalOutputType *MaterialModelOutputs<dim>::get_additional_output() const
1522  {
1523  for (unsigned int i=0; i<additional_outputs.size(); ++i)
1524  {
1525  const AdditionalOutputType *result = dynamic_cast<const AdditionalOutputType *> (additional_outputs[i].get());
1526  if (result)
1527  return result;
1528  }
1529  return nullptr;
1530  }
1531 
1532 
1533  template <int dim>
1535  {
1536  Assert(this->additional_outputs.empty(), ExcMessage("Destination of move needs to be empty!"));
1537  this->additional_outputs = std::move(other.additional_outputs);
1538  }
1539 
1540 
1548 #define ASPECT_REGISTER_MATERIAL_MODEL(classname,name,description) \
1549  template class classname<2>; \
1550  template class classname<3>; \
1551  namespace ASPECT_REGISTER_MATERIAL_MODEL_ ## classname \
1552  { \
1553  aspect::internal::Plugins::RegisterHelper<aspect::MaterialModel::Interface<2>,classname<2>> \
1554  dummy_ ## classname ## _2d (&aspect::MaterialModel::register_material_model<2>, \
1555  name, description); \
1556  aspect::internal::Plugins::RegisterHelper<aspect::MaterialModel::Interface<3>,classname<3>> \
1557  dummy_ ## classname ## _3d (&aspect::MaterialModel::register_material_model<3>, \
1558  name, description); \
1559  }
1560  }
1561 }
1562 
1563 
1564 #endif
void write_plugin_graph(std::ostream &output_stream)
void average_property(const AveragingOperation operation, const FullMatrix< double > &projection_matrix, const FullMatrix< double > &expansion_matrix, std::vector< double > &values_out)
::TrilinosWrappers::MPI::BlockVector BlockVector
Definition: global.h:229
AdditionalInputType * get_additional_input()
Definition: interface.h:1479
std::vector< SymmetricTensor< 2, dim > > strain_rate
Definition: interface.h:415
std::vector< double > thermal_expansion_coefficients
Definition: interface.h:543
virtual void average(const MaterialAveraging::AveragingOperation, const FullMatrix< double > &, const FullMatrix< double > &)
Definition: interface.h:825
MaterialProperties::Property requested_properties
Definition: interface.h:441
std::vector< std::vector< double > > prescribed_field_outputs
Definition: interface.h:1032
std::vector< std::vector< double > > composition
Definition: interface.h:401
std::vector< Point< dim > > position
Definition: interface.h:369
ElasticOutputs(const unsigned int n_points)
Definition: interface.h:1160
AdditionalOutputType * get_additional_output()
Definition: interface.h:1507
std::vector< std::vector< double > > reaction_rates
Definition: interface.h:979
void declare_parameters(ParameterHandler &prm)
bool identifies_single_variable(const Dependence dependence)
void average(const MaterialAveraging::AveragingOperation, const FullMatrix< double > &, const FullMatrix< double > &) override
Definition: interface.h:1083
void average(const AveragingOperation operation, const typename DoFHandler< dim >::active_cell_iterator &cell, const Quadrature< dim > &quadrature_formula, const Mapping< dim > &mapping, MaterialModelOutputs< dim > &values_out)
const std::vector< std::string > names
Definition: interface.h:910
Property operator|=(Property &d1, const Property d2)
Definition: interface.h:233
std::vector< std::vector< double > > reaction_terms
Definition: interface.h:612
EnthalpyOutputs(const unsigned int n_points)
Definition: interface.h:1199
AveragingOperation parse_averaging_operation_name(const std::string &s)
DoFHandler< dim >::active_cell_iterator current_cell
Definition: interface.h:432
Property operator|(const Property d1, const Property d2)
Definition: interface.h:227
AdditionalMaterialOutputsStokesRHS(const unsigned int n_points)
Definition: interface.h:1076
std::vector< double > enthalpies_of_fusion
Definition: interface.h:1219
std::vector< std::unique_ptr< AdditionalMaterialInputs< dim > > > additional_inputs
Definition: interface.h:448
void move_additional_outputs_from(MaterialModelOutputs< dim > &other)
Definition: interface.h:1534
std::vector< Tensor< 1, dim > > pressure_gradient
Definition: interface.h:385
std::vector< std::vector< double > > output_values
Definition: interface.h:907
Interface< dim > * create_material_model(const std::string &model_name)
std::vector< Tensor< 2, dim > > elastic_force
Definition: interface.h:1180
void average(const MaterialAveraging::AveragingOperation operation, const FullMatrix< double > &, const FullMatrix< double > &) override
Definition: interface.h:1167
std::vector< double > entropy_derivative_pressure
Definition: interface.h:566
std::vector< double > thermal_conductivities
Definition: interface.h:553
std::vector< double > entropy_derivative_temperature
Definition: interface.h:574
NonlinearDependence::ModelDependence model_dependence
Definition: interface.h:1382
std::string get_valid_model_names_pattern()
std::vector< Tensor< 1, dim > > velocity
Definition: interface.h:394
virtual void average(const MaterialAveraging::AveragingOperation operation, const FullMatrix< double > &, const FullMatrix< double > &)
Definition: interface.h:1206
void register_material_model(const std::string &name, const std::string &description, void(*declare_parameters_function)(ParameterHandler &), Interface< dim > *(*factory_function)())
void average(const MaterialAveraging::AveragingOperation, const FullMatrix< double > &, const FullMatrix< double > &) override
Definition: interface.h:897
std::vector< std::unique_ptr< AdditionalMaterialOutputs< dim > > > additional_outputs
Definition: interface.h:619