Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

allow ascii data reference profile material model to output seismic v… #1570

Merged
merged 1 commit into from May 12, 2017
Merged

allow ascii data reference profile material model to output seismic v… #1570

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
305 changes: 305 additions & 0 deletions data/adiabatic-conditions/ascii-data/example_seismic_velocities.txt
View file
Open in desktop

Large diffs are not rendered by default.

17 changes: 17 additions & 0 deletions include/aspect/material_model/ascii_reference_profile.h
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
@@ -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
44 changes: 44 additions & 0 deletions tests/burnman_seismic_test/screen-output
View file
Open in desktop
@@ -0,0 +1,44 @@
-----------------------------------------------------------------------------
-----------------------------------------------------------------------------


Setting up GPlates boundary velocity plugin.

Input point 1 spherical coordinates: 1.571 4.870
Input point 1 normalized cartesian coordinates: 0.157 -0.988 0.000
Input point 1 rotated model coordinates: 0.157 -0.988 0.000
Input point 2 spherical coordinates: 1.571 5.240
Input point 2 normalized cartesian coordinates: 0.503 -0.864 0.000
Input point 2 rotated model coordinates: 0.503 -0.864 0.000

Model will be rotated by 0.00 degrees around axis 0.00 0.00 1.00
The ParaView rotation angles are: 0.00 0.00 0.00
The inverse ParaView rotation angles are: 0.00 0.00 0.00

Loading GPlates data boundary file ASPECT_DIR/data/velocity-boundary-conditions/gplates/current_day.gpml.


Loading new velocity file did not succeed.
Assuming constant boundary conditions for rest of model run.

Number of active cells: 3,072 (on 5 levels)
Number of degrees of freedom: 41,280 (25,344+3,264+12,672)

*** Timestep 0: t=0 years
Solving temperature system... 0 iterations.
Rebuilding Stokes preconditioner...
Solving Stokes system... 156+0 iterations.

Postprocessing:
Writing graphical output: output-burnman_seismic_test/solution/solution-00000
RMS, max velocity: 0.319 m/year, 1.61 m/year
Temperature min/avg/max: 273 K, 2082 K, 3500 K
Heat fluxes through boundary parts: -1.315e+06 W, 3.033e+06 W

Termination requested by criterion: end time


+---------------------------------------------+------------+------------+
+---------------------------------+-----------+------------+------------+
+---------------------------------+-----------+------------+------------+