Skip to content

Commit

Permalink
allow ascii data reference profile material model to output seismic v…
Browse files Browse the repository at this point in the history 
…elocities
  • Loading branch information
@jdannberg
jdannberg committed May 11, 2017
1 parent e3819f1 commit a57f69a
Show file tree
Hide file tree
Showing 9 changed files with 22,180 additions and 2 deletions.
305 changes: 305 additions & 0 deletions data/adiabatic-conditions/ascii-data/example_seismic_velocities.txt
View file

Large diffs are not rendered by default.

17 changes: 17 additions & 0 deletions include/aspect/material_model/ascii_reference_profile.h
View file
Expand Up @@ -107,6 +107,14 @@ namespace aspect
* @}
*/

/**
* Add the named outputs for seismic velocities.
*/
virtual
void
create_additional_named_outputs (MaterialModel::MaterialModelOutputs<dim> &out) const;


private:
/**
* Use truncated anelastic approximation?
Expand Down Expand Up @@ -154,6 +162,15 @@ namespace aspect
unsigned int thermal_expansivity_index;
unsigned int specific_heat_index;
unsigned int compressibility_index;

/**
* The column indices of the seismic velocities and their temperature
* derivatives columns in the data file.
*/
unsigned int seismic_vp_index;
unsigned int seismic_vs_index;
unsigned int seismic_dvp_dT_index;
unsigned int seismic_dvs_dT_index;
};

}
Expand Down
10 changes: 9 additions & 1 deletion include/aspect/utilities.h
View file
Expand Up @@ -793,12 +793,20 @@ namespace aspect

/**
* Returns the column index of a column with the given name
* @p column_name. Returns numbers::invalid_unsigned_int if no such
* @p column_name. Throws an exception if no such
* column exists or no names were provided in the file.
*/
unsigned int
get_column_index_from_name(const std::string &column_name) const;

/**
* Returns the column index of a column with the given name
* @p column_name. Returns an invalid unsigned int if no such
* column exists or no names were provided in the file.
*/
unsigned int
maybe_get_column_index_from_name(const std::string &column_name) const;

/**
* Returns a string that contains the name of the column with index
* @p column_index. Returns an empty string if no such
Expand Down
44 changes: 43 additions & 1 deletion source/material_model/ascii_reference_profile.cc
View file
Expand Up @@ -33,7 +33,11 @@ namespace aspect
density_index(numbers::invalid_unsigned_int),
thermal_expansivity_index(numbers::invalid_unsigned_int),
specific_heat_index(numbers::invalid_unsigned_int),
compressibility_index(numbers::invalid_unsigned_int)
compressibility_index(numbers::invalid_unsigned_int),
seismic_vp_index(numbers::invalid_unsigned_int),
seismic_vs_index(numbers::invalid_unsigned_int),
seismic_dvp_dT_index(numbers::invalid_unsigned_int),
seismic_dvs_dT_index(numbers::invalid_unsigned_int)
{}

template <int dim>
Expand All @@ -46,6 +50,14 @@ namespace aspect
thermal_expansivity_index = profile.get_column_index_from_name("thermal_expansivity");
specific_heat_index = profile.get_column_index_from_name("specific_heat");
compressibility_index = profile.get_column_index_from_name("compressibility");

// these are only optional entries in the data file, read them in if they exist,
// but keep the invalid unsigned int entry if the columns do not exist
// the strings have to be all lower case for the lookup class to find them
seismic_vp_index = profile.maybe_get_column_index_from_name("seismic_vp");
seismic_vs_index = profile.maybe_get_column_index_from_name("seismic_vs");
seismic_dvp_dT_index = profile.maybe_get_column_index_from_name("seismic_dvp_dt");
seismic_dvs_dT_index = profile.maybe_get_column_index_from_name("seismic_dvs_dt");
}

template <int dim>
Expand Down Expand Up @@ -91,6 +103,21 @@ namespace aspect

for (unsigned int c=0; c<in.composition[i].size(); ++c)
out.reaction_terms[i][c] = 0.0;

// fill seismic velocities outputs if they exist
if (SeismicAdditionalOutputs<dim> *seismic_out = out.template get_additional_output<SeismicAdditionalOutputs<dim> >())
{
if (seismic_vp_index != numbers::invalid_unsigned_int)
seismic_out->vp[i] = profile.get_data_component(profile_position,seismic_vp_index);
if (seismic_vs_index != numbers::invalid_unsigned_int)
seismic_out->vs[i] = profile.get_data_component(profile_position,seismic_vs_index);
if (seismic_dvp_dT_index != numbers::invalid_unsigned_int)
seismic_out->vp[i] += profile.get_data_component(profile_position,seismic_dvp_dT_index)
* temperature_deviation;
if (seismic_dvs_dT_index != numbers::invalid_unsigned_int)
seismic_out->vs[i] += profile.get_data_component(profile_position,seismic_dvs_dT_index)
* temperature_deviation;
}
}
}

Expand Down Expand Up @@ -199,6 +226,21 @@ namespace aspect
this->model_dependence.thermal_conductivity = NonlinearDependence::none;

}



template <int dim>
void
AsciiReferenceProfile<dim>::create_additional_named_outputs (MaterialModel::MaterialModelOutputs<dim> &out) const
{
if (out.template get_additional_output<SeismicAdditionalOutputs<dim> >() == NULL)
{
const unsigned int n_points = out.viscosities.size();
out.additional_outputs.push_back(
std_cxx11::shared_ptr<MaterialModel::AdditionalMaterialOutputs<dim> >
(new MaterialModel::SeismicAdditionalOutputs<dim> (n_points)));
}
}
}
}

Expand Down
16 changes: 16 additions & 0 deletions source/utilities.cc
View file
Expand Up @@ -1805,6 +1805,22 @@ namespace aspect
return lookup->get_column_index_from_name(column_name);
}

template <int dim>
unsigned int
AsciiDataProfile<dim>::maybe_get_column_index_from_name(const std::string &column_name) const
{
try
{
// read the entries in if they exist
return lookup->get_column_index_from_name(column_name);
}
catch (...)
{
// return an invalid unsigned int entry if the column does not exist
return numbers::invalid_unsigned_int;
}
}

template <int dim>
std::string
AsciiDataProfile<dim>::get_column_name_from_index(const unsigned int column_index) const
Expand Down
149 changes: 149 additions & 0 deletions tests/burnman_seismic_test.prm
View file
@@ -0,0 +1,149 @@
# A cookbook based on the burman_test cookbook, but modified testing
# the output of seismic velocity data generated by burnman and reading
# them in through an ascii data file.

set Dimension = 2
set Use years in output instead of seconds = true
set End time = 0
set Linear solver tolerance = 1e-6

subsection Material model
set Model name = ascii reference profile

subsection Ascii reference profile
set Thermal viscosity exponent = 10.0
set Viscosity prefactors = 1.0, 0.1, 1.0, 10.0

subsection Ascii data model
set Data directory = $ASPECT_SOURCE_DIR/data/adiabatic-conditions/ascii-data/
set Data file name = example_seismic_velocities.txt
end
end

set Material averaging = harmonic average
end

# The geometry is a spherical shell with the inner and
# outer radius of the Earth's mantle.
subsection Geometry model
set Model name = spherical shell

subsection Spherical shell
set Inner radius = 3481000
set Outer radius = 6336000
end
end

subsection Adiabatic conditions model
set Model name = ascii data

subsection Ascii data model
set Data directory = $ASPECT_SOURCE_DIR/data/adiabatic-conditions/ascii-data/
set Data file name = example_seismic_velocities.txt
end
end


# The gravity model reads an ascii data file that contains a gravity
# profile consistent with the material properties computed using the
# Birch-Murnaghan equation of state. It automatically uses the same file
# specified in the adiabatic conditions model.
subsection Gravity model
set Model name = ascii data
subsection Ascii data model
set Data directory = $ASPECT_SOURCE_DIR/data/adiabatic-conditions/ascii-data/
set Data file name = example_seismic_velocities.txt
end
end

# The model uses the anelastic liquid approximation.
subsection Formulation
set Formulation = anelastic liquid approximation
end

# The present-day plate velocities are imposed on the upper model
# boundary.
subsection Model settings
set Prescribed velocity boundary indicators = outer:gplates
set Tangential velocity boundary indicators = inner

set Fixed temperature boundary indicators = inner, outer
end


subsection Boundary velocity model
subsection GPlates model
set Data directory = $ASPECT_SOURCE_DIR/data/velocity-boundary-conditions/gplates/
set Velocity file name = current_day.gpml
set Data file time step = 1e6
set Point one = 1.5708,4.87
set Point two = 1.5708,5.24
end
end


subsection Boundary temperature model
set Model name = spherical constant
subsection Spherical constant
set Inner temperature = 3500
set Outer temperature = 273
end
end

# The initial temperature is an adiabatic profile, which is
# taken from the adabatic conditions model. In this case, this is
# is the ascii data model that loads a file containing profiles
# corresponding to the Birch-Murnaghan equation of state.
# In addition to the adiabatic profile, there are thermal boundary
# layers at the sirface and the core-mantle boundary.
subsection Initial temperature model
set List of model names = adiabatic, function

subsection Adiabatic
set Age top boundary layer = 5e7
end

subsection Function
set Variable names = x,z
set Function constants = pi=3.1415926
set Function expression = if(sqrt(x*x+(z-4500000)*(z-4500000))<1000000,800,0)
end
end


subsection Mesh refinement
set Refinement fraction = 0.4
set Coarsening fraction = 0.05
set Initial adaptive refinement = 0
set Initial global refinement = 4
set Strategy = temperature
set Time steps between mesh refinement = 5
end


subsection Heating model
set List of model names = adiabatic heating, shear heating

subsection Adiabatic heating
set Use simplified adiabatic heating = true
end
end

subsection Postprocess
set List of postprocessors = visualization, velocity statistics, temperature statistics, heat flux statistics

subsection Visualization
set Output format = gnuplot
set Time between graphical output = 0
set Number of grouped files = 0
set List of output variables = material properties, gravity, adiabat, nonadiabatic temperature, nonadiabatic pressure, named additional outputs
end

subsection Depth average
set Time between graphical output = 1e6
end
end

subsection Checkpointing
set Time between checkpoint = 1800
end

0 comments on commit a57f69a

Please sign in to comment.