Making shear and P wave speed user defined quantities for seismic for…

include/bcs/NonReflectingBC.h

@@ -72,13 +72,13 @@ class NonReflectingBC : public IntegratedBC
   const Real _alpha;
 
   /// Density of the soil
-  const MaterialProperty<Real> & _density;
-
-  /// Shear wave speed of the soil
-  const MaterialProperty<Real> & _shear_wave_speed;
+  const Real & _density;
 
   /// P wave speed of the soil
-  const MaterialProperty<Real> & _P_wave_speed;
+  const Real & _p_wave_speed;
+
+  /// Shear wave speed of the soil
+  const Real & _shear_wave_speed;
 };
 
 #endif // NONREFLECTINGBC_H

include/bcs/SeismicForce.h

@@ -46,14 +46,14 @@ class SeismicForce : public IntegratedBC
   /// alpha Parameter for HHT time integration scheme
   const Real _alpha;
 
-  /// density of the soil
-  const MaterialProperty<Real> & _density;
+  /// density of the underlying bedrock
+  const Real & _density;
 
-  /// shear wave speed of the soil
-  const MaterialProperty<Real> & _shear_wave_speed;
+  /// P wave speed of the underlying bedrock
+  const Real & _p_wave_speed;
 
-  /// P wave speed of the soil
-  const MaterialProperty<Real> & _P_wave_speed;
+  /// shear wave speed of the underlying bedrock
+  const Real & _shear_wave_speed;
 };
 
 #endif // SEISMICFORCE_H

src/actions/SeismicForceAction.C

@@ -48,6 +48,12 @@ validParams<SeismicForceAction>()
                                              "input_component.");
   params.addParam<Real>(
       "alpha", 0.0, "The alpha parameter required for HHT time integration scheme.");
+  params.addRequiredRangeCheckedParam<Real>(
+      "density", "density>0.0", "Density of the underlying bedrock.");
+  params.addRequiredRangeCheckedParam<Real>(
+      "p_wave_speed", "p_wave_speed>0.0", "P-wave speed of the underlying bedrock.");
+  params.addRequiredRangeCheckedParam<Real>(
+      "shear_wave_speed", "shear_wave_speed>0.0", "shear wave speed of the underlying bedrock.");
   return params;
 }
 
@@ -97,6 +103,9 @@ SeismicForceAction::act()
   InputParameters params = _factory.getValidParams("SeismicForce");
   params.set<std::vector<BoundaryName>>("boundary") = _boundary;
   params.set<Real>("alpha") = getParam<Real>("alpha");
+  params.set<Real>("density") = getParam<Real>("density");
+  params.set<Real>("p_wave_speed") = getParam<Real>("p_wave_speed");
+  params.set<Real>("shear_wave_speed") = getParam<Real>("shear_wave_speed");
 
   // Loop through each component and create SeismicForce object
   std::string prefix = "mastodon_seismic_force";

src/bcs/NonReflectingBC.C

@@ -52,6 +52,11 @@ NonReflectingBC::commonParameters()
   params.addRequiredParam<Real>("beta", "The beta parameter for Newmark time integration.");
   params.addRequiredParam<Real>("gamma", "The gamma parameter for Newmark time integration.");
   params.addParam<Real>("alpha", 0.0, "The alpha parameter for HHT time integration.");
+  params.addRequiredRangeCheckedParam<Real>("density", "density>0.0", "Density of the material.");
+  params.addRequiredRangeCheckedParam<Real>(
+      "p_wave_speed", "p_wave_speed>0.0", "P-wave speed of the material.");
+  params.addRequiredRangeCheckedParam<Real>(
+      "shear_wave_speed", "shear_wave_speed>0.0", "shear wave speed of the material.");
   return params;
 }
 
@@ -67,9 +72,9 @@ NonReflectingBC::NonReflectingBC(const InputParameters & parameters)
     _beta(getParam<Real>("beta")),
     _gamma(getParam<Real>("gamma")),
     _alpha(getParam<Real>("alpha")),
-    _density(getMaterialProperty<Real>("density")),
-    _shear_wave_speed(getMaterialProperty<Real>("shear_wave_speed")),
-    _P_wave_speed(getMaterialProperty<Real>("P_wave_speed"))
+    _density(getParam<Real>("density")),
+    _p_wave_speed(getParam<Real>("p_wave_speed")),
+    _shear_wave_speed(getParam<Real>("shear_wave_speed"))
 {
 
   // Error checking on variable vectors
@@ -121,18 +126,18 @@ NonReflectingBC::computeQpResidual()
   }
   // residual is test[i][_qp] *( density* V_p * normal component of velocity +
   // density * V_s* tangential component of velocity)
-  return _test[_i][_qp] * _density[_qp] *
-         (_P_wave_speed[_qp] * normal_vel * _normals[_qp](_component) +
-          _shear_wave_speed[_qp] * (vel[_component] - normal_vel * _normals[_qp](_component)));
+  return _test[_i][_qp] * _density *
+         (_p_wave_speed * normal_vel * _normals[_qp](_component) +
+          _shear_wave_speed * (vel[_component] - normal_vel * _normals[_qp](_component)));
 }
 
 Real
 NonReflectingBC::computeQpJacobian()
 {
-  return _test[_i][_qp] * _density[_qp] *
-         ((_P_wave_speed[_qp] - _shear_wave_speed[_qp]) * _normals[_qp](_component) *
+  return _test[_i][_qp] * _density *
+         ((_p_wave_speed - _shear_wave_speed) * _normals[_qp](_component) *
               _normals[_qp](_component) +
-          _shear_wave_speed[_qp]) *
+          _shear_wave_speed) *
          (1. + _alpha) * _gamma / _beta / _dt * _phi[_j][_qp];
 }
 
@@ -150,7 +155,7 @@ NonReflectingBC::computeQpOffDiagJacobian(unsigned int jvar)
     }
 
   if (active)
-    return _test[_i][_qp] * _density[_qp] * (_P_wave_speed[_qp] - _shear_wave_speed[_qp]) *
+    return _test[_i][_qp] * _density * (_p_wave_speed - _shear_wave_speed) *
            _normals[_qp](_component) * _normals[_qp](coupled_component) * (1. + _alpha) * _gamma /
            _beta / _dt * _phi[_j][_qp];
 

src/bcs/SeismicForce.C

@@ -33,6 +33,12 @@ validParams<SeismicForce>()
   params.addParam<FunctionName>("velocity_function",
                                 "The function that describes the input ground velocity.");
   params.addCoupledVar("velocity", "The variable that describes the input velocity.");
+  params.addRequiredRangeCheckedParam<Real>(
+      "density", "density>0.0", "Density of the underlying bedrock.");
+  params.addRequiredRangeCheckedParam<Real>(
+      "p_wave_speed", "p_wave_speed>0.0", "P-wave speed of the underlying bedrock.");
+  params.addRequiredRangeCheckedParam<Real>(
+      "shear_wave_speed", "shear_wave_speed>0.0", "shear wave speed of the underlying bedrock.");
   params.addParam<Real>(
       "alpha", 0.0, "The alpha parameter required for HHT time integration scheme.");
   params.set<bool>("use_displaced_mesh") = true;
@@ -45,9 +51,9 @@ SeismicForce::SeismicForce(const InputParameters & parameters)
     _vel_component(getParam<unsigned int>("vel_component")),
     _factor(getParam<Real>("factor")),
     _alpha(getParam<Real>("alpha")),
-    _density(getMaterialProperty<Real>("density")),
-    _shear_wave_speed(getMaterialProperty<Real>("shear_wave_speed")),
-    _P_wave_speed(getMaterialProperty<Real>("P_wave_speed"))
+    _density(getParam<Real>("density")),
+    _p_wave_speed(getParam<Real>("p_wave_speed")),
+    _shear_wave_speed(getParam<Real>("shear_wave_speed"))
 {
   if (_component >= _mesh.dimension())
     mooseError("Invalid value for 'component' (",
@@ -97,7 +103,7 @@ SeismicForce::computeQpResidual()
   for (unsigned int i = 0; i < _mesh.dimension(); i++)
     tangential_vel[i] = vel[i] - normal_vel * _normals[_qp](i);
 
-  return _factor * _test[_i][_qp] * _density[_qp] *
-         (_P_wave_speed[_qp] * normal_vel * _normals[_qp](_component) +
-          _shear_wave_speed[_qp] * tangential_vel[_component]);
+  return _factor * _test[_i][_qp] * _density *
+         (_p_wave_speed * normal_vel * _normals[_qp](_component) +
+          _shear_wave_speed * tangential_vel[_component]);
 }

test/tests/bcs/nonreflecting_bc/non_reflecting_bc_error.i

@@ -13,17 +13,20 @@
 []
 
 [BCs]
- [./test]
-   type = NonReflectingBC
-   displacements = 'disp_x disp_y'
-   accelerations = 'disp_x disp_y'
-   velocities = 'disp_x disp_y'
-   variable = 'disp_x'
-   component = 0
-   boundary = 'left'
-   beta = 1
-   gamma = 1
- [../]
+  [./test]
+    type = NonReflectingBC
+    displacements = 'disp_x disp_y'
+    accelerations = 'disp_x disp_y'
+    velocities = 'disp_x disp_y'
+    variable = 'disp_x'
+    component = 0
+    boundary = 'left'
+    beta = 1
+    gamma = 1
+    shear_wave_speed = 1.0
+    p_wave_speed = 3.5
+    density = 1.0
+  [../]
 []
 
 [Problem]

test/tests/bcs/nonreflecting_bc/non_reflecting_bc_test.i

@@ -144,23 +144,26 @@
 
 
 [BCs]
- [./front]
-   type = PresetDisplacement
-   variable = disp_x
-   function = x_vel
-   boundary = front
-   acceleration = accel_x
-   velocity = vel_x
-   beta = 0.25
- [../]
- [./NonReflectingBC]
+  [./front]
+    type = PresetDisplacement
+    variable = disp_x
+    function = x_vel
+    boundary = front
+    acceleration = accel_x
+    velocity = vel_x
+    beta = 0.25
+  [../]
+  [./NonReflectingBC]
     [./back]
       displacements = 'disp_x disp_y disp_z'
       velocities = 'vel_x vel_y vel_z'
       accelerations = 'accel_x accel_y accel_z'
       beta = 0.25
       gamma = 0.5
       boundary = 'back'
+      shear_wave_speed = 10.0
+      p_wave_speed = 16.329931618554521
+      density = 1.0
     [../]
   [../]
   [./Periodic]

test/tests/bcs/seismic_force/seismic_force_action_error.i

@@ -25,6 +25,9 @@
       displacements = 'disp_x disp_y'
       input_components = '0 1'
       boundary = 'front'
+      shear_wave_speed = 1.0
+      p_wave_speed = 3.5
+      density = 1.0
     [../]
   [../]
 []

test/tests/bcs/seismic_force/seismic_force_error.i

@@ -24,6 +24,9 @@
     type = SeismicForce
     variable = 'disp_x'
     boundary = 'left'
+    shear_wave_speed = 1.0
+    p_wave_speed = 3.5
+    density = 1.0
   [../]
 []
 

test/tests/bcs/seismic_force/seismic_force_test_x_y.i

@@ -151,6 +151,9 @@
       input_components = '0 1'
       boundary = 'front'
       velocity_functions = 'x_vel y_vel' # input velocity functions
+      shear_wave_speed = 1.0
+      p_wave_speed = 3.5
+      density = 1.0
     [../]
   [../]
   [./back_x]

test/tests/bcs/seismic_force/seismic_force_test_z.i

@@ -191,21 +191,27 @@
 
 [BCs]
   [./SeismicForce]
-   [./top_stress]
-       displacements = 'disp_x disp_y disp_z'
-       input_components = '2'
-       boundary = 'front'
-       velocity_functions = 'z_vel'
-   [../]
+    [./top_stress]
+      displacements = 'disp_x disp_y disp_z'
+      input_components = '2'
+      boundary = 'front'
+      velocity_functions = 'z_vel'
+      shear_wave_speed = 3718.204815052883833
+      p_wave_speed = 7381.527628151251754
+      density = 4.9419e-3
+    [../]
   [../]
   [./NonReflectingBC]
     [./back]
-        displacements = 'disp_x disp_y disp_z'
-        velocities = 'vel_x vel_y vel_z'
-        accelerations = 'accel_x accel_y accel_z'
-        beta = 0.25
-        gamma = 0.5
-        boundary = 'back'
+      displacements = 'disp_x disp_y disp_z'
+      velocities = 'vel_x vel_y vel_z'
+      accelerations = 'accel_x accel_y accel_z'
+      beta = 0.25
+      gamma = 0.5
+      boundary = 'back'
+      shear_wave_speed = 3718.204815052883833
+      p_wave_speed = 7381.527628151251754
+      density = 4.9419e-3
     [../]
   [../]
   [./Periodic]