Change seismic velocity default to signaling nan

Also remove seismic velocities from depth average if not provided

doc/modules/changes/20180427_gassmoeller_2

@@ -0,0 +1,5 @@
+Changed: The default value for AdditionalSeismicOutputs was changed from -1 to
+a signaling Nan. This ensures every plugin that creates seismic outputs needs
+to use a material model that actually fills those outputs.
+<br>
+(Rene Gassmoeller, 2018/04/27)

source/material_model/interface.cc

@@ -763,8 +763,8 @@ namespace aspect
     SeismicAdditionalOutputs<dim>::SeismicAdditionalOutputs (const unsigned int n_points)
       :
       NamedAdditionalMaterialOutputs<dim>(make_seismic_additional_outputs_names()),
-      vs(n_points, -1.0),
-      vp(n_points, -1.0)
+      vs(n_points, numbers::signaling_nan<double>()),
+      vp(n_points, numbers::signaling_nan<double>())
     {}
 
 

source/postprocess/depth_average.cc

@@ -366,11 +366,40 @@ namespace aspect
             if ( output_all_variables || std::find( output_variables.begin(), output_variables.end(), "sinking velocity") != output_variables.end() )
               variables.push_back("sinking_velocity");
 
-            if ( output_all_variables || std::find( output_variables.begin(), output_variables.end(), "Vs") != output_variables.end() )
-              variables.push_back("Vs");
+            // handle seismic velocities, because they may, or may not be provided by the material model
+            {
+              MaterialModel::MaterialModelOutputs<dim> out(1, this->n_compositional_fields());
+              this->get_material_model().create_additional_named_outputs(out);
+
+              const bool material_model_provides_seismic_output =
+                (out.template get_additional_output<MaterialModel::SeismicAdditionalOutputs<dim> >() != 0);
+
+              if (output_all_variables && material_model_provides_seismic_output)
+                {
+                  variables.push_back("Vs");
+                  variables.push_back("Vp");
+                }
 
-            if ( output_all_variables || std::find( output_variables.begin(), output_variables.end(), "Vp") != output_variables.end() )
-              variables.push_back("Vp");
+              if (std::find( output_variables.begin(), output_variables.end(), "Vs") != output_variables.end())
+                {
+                  AssertThrow(material_model_provides_seismic_output,
+                              ExcMessage("You requested seismic velocities from the 'Depth average' postprocessor, "
+                                         "but the material model does not provide seismic velocities. Either remove 'Vs' and "
+                                         "'Vp' from the 'List of output variables' parameter, or use a material model that "
+                                         "provides these velocities."));
+                  variables.push_back("Vs");
+                }
+
+              if (std::find( output_variables.begin(), output_variables.end(), "Vp") != output_variables.end())
+                {
+                  AssertThrow(material_model_provides_seismic_output,
+                              ExcMessage("You requested seismic velocities from the 'Depth average' postprocessor, "
+                                         "but the material model does not provide seismic velocities. Either remove 'Vs' and "
+                                         "'Vp' from the 'List of output variables' parameter, or use a material model that "
+                                         "provides these velocities."));
+                  variables.push_back("Vp");
+                }
+            }
 
             if ( output_all_variables || std::find( output_variables.begin(), output_variables.end(), "viscosity") != output_variables.end() )
               variables.push_back("viscosity");

source/postprocess/visualization/seismic_anomalies.cc

@@ -405,6 +405,19 @@ namespace aspect
       void
       SeismicVsAnomaly<dim>::parse_parameters (ParameterHandler &prm)
       {
+        // First make sure the material model actually provides seismic velocities
+        MaterialModel::MaterialModelOutputs<dim> out(1, this->n_compositional_fields());
+        this->get_material_model().create_additional_named_outputs(out);
+
+        const bool material_model_provides_seismic_output =
+          (out.template get_additional_output<MaterialModel::SeismicAdditionalOutputs<dim> >() != 0);
+
+        AssertThrow(material_model_provides_seismic_output,
+                    ExcMessage("You requested the 'Seismic' postprocessor, "
+                               "but the material model does not provide seismic velocities. Either remove 'Vs' and "
+                               "'Vp' from the 'List of output variables' parameter, or use a material model that "
+                               "provides these velocities."));
+
         prm.enter_subsection("Postprocess");
         {
           prm.enter_subsection("Visualization");
@@ -468,13 +481,25 @@ namespace aspect
       void
       SeismicVpAnomaly<dim>::parse_parameters (ParameterHandler &prm)
       {
+        // First make sure the material model actually provides seismic velocities
+        MaterialModel::MaterialModelOutputs<dim> out(1, this->n_compositional_fields());
+        this->get_material_model().create_additional_named_outputs(out);
+
+        const bool material_model_provides_seismic_output =
+          (out.template get_additional_output<MaterialModel::SeismicAdditionalOutputs<dim> >() != 0);
+
+        AssertThrow(material_model_provides_seismic_output,
+                    ExcMessage("You requested the 'Seismic' postprocessor, "
+                               "but the material model does not provide seismic velocities. Either remove 'Vs' and "
+                               "'Vp' from the 'List of output variables' parameter, or use a material model that "
+                               "provides these velocities."));
+
         prm.enter_subsection("Postprocess");
         {
           prm.enter_subsection("Visualization");
           {
             prm.enter_subsection("Vp anomaly");
             {
-
               // Average velocity scheme
               if (prm.get ("Average velocity scheme") == "reference profile")
                 average_velocity_scheme = reference_profile;
@@ -522,7 +547,8 @@ namespace aspect
                                                   "capture step changes in velocities, but small enough to "
                                                   "maintain a reasonable number of evaluation points per slice. "
                                                   "Bear in mind that lateral averaging subsamples the "
-                                                  "finite element mesh.")
+                                                  "finite element mesh. Note that this plugin requires a "
+                                                  "material model that provides seismic velocities.")
 
       ASPECT_REGISTER_VISUALIZATION_POSTPROCESSOR(SeismicVpAnomaly,
                                                   "Vp anomaly",
@@ -544,7 +570,8 @@ namespace aspect
                                                   "capture step changes in velocities, but small enough to "
                                                   "maintain a reasonable number of evaluation points per slice. "
                                                   "Bear in mind that lateral averaging subsamples the "
-                                                  "finite element mesh.")
+                                                  "finite element mesh. Note that this plugin requires a "
+                                                  "material model that provides seismic velocities.")
     }
   }
 }