Put fluid density back.

LynxDensityCompressible still need some changes.

include/materials/LynxDensityBase.h

@@ -41,11 +41,11 @@ protected:
 
   bool _has_gravity;
   Real _g;
-  // const std::vector<Real> _fluid_density;
+  const std::vector<Real> _fluid_density;
   const std::vector<Real> _solid_density;
 
   MaterialProperty<RealVectorValue> & _gravity;
-  // MaterialProperty<Real> & _rho_f;
+  MaterialProperty<Real> & _rho_f;
   MaterialProperty<Real> & _rho_s;
   MaterialProperty<Real> & _rho_b;
   MaterialProperty<Real> & _reference_rho_b;

include/materials/LynxDensityThermal.h

@@ -41,7 +41,7 @@ protected:
 
   MaterialProperty<Real> & _drho_dtemp;
   MaterialProperty<Real> & _dinvrho_dtemp;
-  // const std::vector<Real> _beta_fluid;
+  const std::vector<Real> _beta_fluid;
   const std::vector<Real> _beta_solid;
   Real _temp_ref;
   Function * const _temp_ref_fct;

src/materials/LynxDensityBase.C

@@ -30,8 +30,7 @@ validParams<LynxDensityBase>()
   params.addCoupledVar("porosity", "The porosity auxiliary variable.");
   params.addParam<bool>("has_gravity", false, "Model with gravity on?");
   params.addParam<Real>("gravity_acceleration", 9.81, "The magnitude of the gravity acceleration.");
-  // params.addParam<std::vector<Real>>("fluid_density", std::vector<Real>(1, 0.0), "The fluid
-  // density.");
+  params.addParam<std::vector<Real>>("fluid_density", "The fluid density.");
   params.addParam<std::vector<Real>>("solid_density", "The solid density.");
   return params;
 }
@@ -41,11 +40,12 @@ LynxDensityBase::LynxDensityBase(const InputParameters & parameters)
     _porosity(isCoupled("porosity") ? coupledValue("porosity") : _zero),
     _has_gravity(getParam<bool>("has_gravity")),
     _g(_has_gravity ? getParam<Real>("gravity_acceleration") : 0.0),
-    // _fluid_density(getLynxParam<Real>("fluid_density")),
+    _fluid_density(isParamValid("fluid_density") ? getLynxParam<Real>("fluid_density")
+                                                 : std::vector<Real>(_n_composition, 0.0)),
     _solid_density(isParamValid("solid_density") ? getLynxParam<Real>("solid_density")
                                                  : std::vector<Real>(_n_composition, 0.0)),
     _gravity(declareProperty<RealVectorValue>("gravity_vector")),
-    // _rho_f(declareProperty<Real>("fluid_density")),
+    _rho_f(declareProperty<Real>("fluid_density")),
     _rho_s(declareProperty<Real>("solid_density")),
     _rho_b(declareProperty<Real>("bulk_density")),
     _reference_rho_b(declareProperty<Real>("reference_bulk_density"))

src/materials/LynxDensityConstant.C

@@ -40,9 +40,8 @@ void
 LynxDensityConstant::computeQpProperties()
 {
   computeQpGravity();
-  // _rho_f[_qp] = averageProperty(_fluid_density);
+  _rho_f[_qp] = averageProperty(_fluid_density);
   _rho_s[_qp] = averageProperty(_solid_density);
-  // _rho_b[_qp] = _porosity[_qp] * _rho_f[_qp] + (1.0 - _porosity[_qp]) * _rho_s[_qp];
-  _rho_b[_qp] = _rho_s[_qp];
-  _reference_rho_b[_qp] = _rho_s[_qp];
+  _rho_b[_qp] = _porosity[_qp] * _rho_f[_qp] + (1.0 - _porosity[_qp]) * _rho_s[_qp];
+  _reference_rho_b[_qp] = _rho_b[_qp];
 }

src/materials/LynxDensityThermal.C

@@ -29,9 +29,7 @@ validParams<LynxDensityThermal>()
   InputParameters params = validParams<LynxDensityBase>();
   params.addClassDescription(
       "Material calculating densities as a simple linear function of temperature.");
-  params.addRequiredCoupledVar("temperature", "The temperature variable.");
-  // params.addParam<std::vector<Real>>("beta_fluid", std::vector<Real>(1, 0.0), "The fluid thermal
-  // expansion coefficient.");
+  params.addParam<std::vector<Real>>("beta_fluid", "The fluid thermal expansion coefficient.");
   params.addParam<std::vector<Real>>("beta_solid", "The solid thermal expansion coefficient.");
   params.addParam<Real>("reference_temperature", 0.0, "The reference temperature.");
   params.addParam<FunctionName>("reference_temperature_fct",
@@ -44,7 +42,8 @@ LynxDensityThermal::LynxDensityThermal(const InputParameters & parameters)
     _temperature(coupledValue("temperature")),
     _drho_dtemp(declareProperty<Real>("drho_dtemp")),
     _dinvrho_dtemp(declareProperty<Real>("dinvrho_dtemp")),
-    // _beta_fluid(getLynxParam<Real>("beta_fluid")),
+    _beta_fluid(isParamValid("beta_fluid") ? getLynxParam<Real>("beta_fluid")
+                                           : std::vector<Real>(_n_composition, 0.0)),
     _beta_solid(isParamValid("beta_solid") ? getLynxParam<Real>("beta_solid")
                                            : std::vector<Real>(_n_composition, 0.0)),
     _temp_ref(getParam<Real>("reference_temperature")),
@@ -63,20 +62,18 @@ LynxDensityThermal::computeQpProperties()
   if (_temp_ref_fct)
     temp_ref = _temp_ref_fct->value(_t, _q_point[_qp]);
 
-  // _rho_f[_qp] = averageProperty(_fluid_density) * (1.0 - averageProperty(_beta_fluid) *
-  // (_temperature[_qp] - _reference_temperature));
+  _rho_f[_qp] = averageProperty(_fluid_density) *
+                (1.0 - averageProperty(_beta_fluid) * (_temperature[_qp] - temp_ref));
   _rho_s[_qp] = averageProperty(_solid_density) *
                 (1.0 - averageProperty(_beta_solid) * (_temperature[_qp] - temp_ref));
 
-  // _rho_b[_qp] = _porosity[_qp] * _rho_f[_qp] + (1.0 - _porosity[_qp]) * _rho_s[_qp];
-  _rho_b[_qp] = _rho_s[_qp];
+  _rho_b[_qp] = _porosity[_qp] * _rho_f[_qp] + (1.0 - _porosity[_qp]) * _rho_s[_qp];
   _reference_rho_b[_qp] = _rho_b[_qp];
 
-  // Real drho_f_dtemp = -1.0 * averageProperty(_beta_fluid) * averageProperty(_fluid_density);
+  Real drho_f_dtemp = -1.0 * averageProperty(_beta_fluid) * averageProperty(_fluid_density);
   Real drho_s_dtemp = -1.0 * averageProperty(_beta_solid) * averageProperty(_solid_density);
 
-  // _drho_dtemp[_qp] = _porosity[_qp] * drho_f_dtemp + (1.0 - _porosity[_qp]) * drho_s_dtemp;
-  _drho_dtemp[_qp] = drho_s_dtemp;
+  _drho_dtemp[_qp] = _porosity[_qp] * drho_f_dtemp + (1.0 - _porosity[_qp]) * drho_s_dtemp;
 
   _dinvrho_dtemp[_qp] =
       (_rho_b[_qp] > 0.0) ? -1.0 / std::pow(_rho_b[_qp], 2.0) * _drho_dtemp[_qp] : 0.0;