Skip to content

Commit

Permalink
follow
Browse files Browse the repository at this point in the history
  • Loading branch information
@lhy11009
lhy11009 committed Nov 7, 2024
1 parent d07d104 commit 35030e0
Show file tree
Hide file tree
Showing 8 changed files with 361 additions and 10 deletions.
14 changes: 7 additions & 7 deletions source/material_model/utilities.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1163,7 +1163,7 @@ namespace aspect
void
PhaseFunctionDiscrete<dim>::initialize()
{
AssertThrow (this->introspection().get_number_of_fields_of_type(CompositionalFieldDescription::chemical_composition) == material_file_names.size(),
AssertThrow (this->introspection().get_number_of_fields_of_type(CompositionalFieldDescription::chemical_composition)+1 == material_file_names.size(),
ExcMessage(" The 'discrete phase function' requires that the number of compositional fields "
"matches the number of lookup tables."));

Expand Down Expand Up @@ -1200,15 +1200,15 @@ namespace aspect
double function_value;

// lookup the most abundant phase
unsigned n_phases_composition = 0;
unsigned n_comp = 0;
for (auto i : this->introspection().chemical_composition_field_indices())
unsigned int start_phase_transition_index = 0;
unsigned int n_comp = 0;
for (unsigned i = 0 ; i < this->introspection().n_chemical_composition_fields() + 1 ; i++)
{
n_phases_composition += n_phases_per_chemical_composition[i];
if (in.phase_index < n_phases_composition){
if (in.phase_index < start_phase_transition_index + n_phase_transitions_per_chemical_composition[i]){
n_comp = i;
break;
}
start_phase_transition_index += n_phase_transitions_per_chemical_composition[i];
}

Assert (in.temperature >= minimum_temperature[n_comp] && in.temperature < maximum_temperature[n_comp], ExcInternalError());
Expand All @@ -1227,7 +1227,7 @@ namespace aspect
const unsigned maph = unsigned(material_lookup[n_comp]->get_data(temperature_pressure, 7));

// determine the value of phase function
if (in.phase_index < maph)
if (in.phase_index - start_phase_transition_index < maph)
function_value = 1.0;
else
function_value = 0.0;
Expand Down
8 changes: 5 additions & 3 deletions tests/visco_plastic_phases_discrete.prm
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,8 @@
# This prm file is used to generate a phase diagram.
# and test the implementation of discrete phase functions
# in the visco-plastic material model.
# In the outputs, the material statistics are checked for
# correst viscosity values.

set Dimension = 2
set CFL number = 1.0
Expand Down Expand Up @@ -95,12 +97,12 @@ subsection Material model
set Densities = background: 3300.0|3394.4|3442.1|3453.2|3617.6|3691.5|3774.7|3929.1, spharz: 3235.0|3372.3|3441.7|3441.7|3680.8|3717.8|3759.4|3836.6
set Heat capacities = 1250.0
set Thermal expansivities = 0.0
# The options to look up the dominant phase
set Use dominant phase for viscosity = true
set Viscous flow law = diffusion
set Viscosity averaging scheme = harmonic
# for discrete phase function
set Data directory = $ASPECT_SOURCE_DIR/data/material-model/entropy-table/pyrtable/
set Material file names = material_table_temperature_pressure_small.txt
set Viscous flow law = diffusion
set Viscosity averaging scheme = harmonic
set Phase transition phases = background:1
set Prefactors for diffusion creep = background:5e-21|5e-19, spharz:5e-21
set Grain size exponents for diffusion creep = 0.0
Expand Down
147 changes: 147 additions & 0 deletions tests/visco_plastic_phases_discrete_comp.prm
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,147 @@
# This prm file is used to generate a phase diagram.
# and test the implementation of discrete phase functions
# in the visco-plastic material model. Specifically, this test
# checks the implementation for multiple compositions.
# In the outputs, the material statistics are checked for
# correst viscosity values.

set Dimension = 2
set CFL number = 1.0
set End time = 0
set Adiabatic surface temperature = 1673.0
set Use years in output instead of seconds = true
set Nonlinear solver scheme = no Advection, no Stokes

# Model geometry
subsection Geometry model
set Model name = box

subsection Box
set X repetitions = 50
set Y repetitions = 50
set X extent = 800e3
set Y extent = 800e3
end
end

# Mesh refinement specifications
subsection Mesh refinement
set Initial adaptive refinement = 0
set Initial global refinement = 0
set Time steps between mesh refinement = 0
end

subsection Gravity model
set Model name = vertical

subsection Vertical
set Magnitude = 10.0
end
end

# Temperature boundary and initial conditions
# The initial temperature increases linearly from
# the left to the right boundary.
subsection Boundary temperature model
set Fixed temperature boundary indicators = left, right
set List of model names = box

subsection Box
set Left temperature = 273
set Right temperature = 2273
end
end

subsection Initial temperature model
set Model name = function

subsection Function
set Coordinate system = cartesian
set Variable names = x, y
set Function constants = XMAX=800e3, Tl=273.0, Tr=2273
set Function expression = Tl * (x - XMAX)/(-XMAX) + Tr * x / XMAX
end
end

# Comment the following 3 sections for steinberg model and pyrolitic lookup table
# Fields of composition
subsection Compositional fields
set Number of fields = 1
set Names of fields = spharz
set Compositional field methods = field
end

# Initial composition model
subsection Initial composition model
set List of model names = function

subsection Function
set Coordinate system = cartesian
set Function expression = 1.0
end
end

# Boundary composition model
subsection Boundary composition model
set List of model names = box
set Fixed composition boundary indicators = top, bottom, left, right

subsection Box
set Bottom composition = 1.0
set Left composition = 1.0
set Right composition = 1.0
set Top composition = 1.0
end
end

# Value for material model
# A set of parameters are assigned to the conventional
# phase function interface. A single phase with index 1
# is added to the background composition with a smaller
# viscosity by using the discrete phase transition.
subsection Material model
set Model name = visco plastic

subsection Visco Plastic
set Reference temperature = 273
set Maximum viscosity = 1e25
set Minimum viscosity = 1e13
set Phase transition depths = background:410e3|520e3|560e3|670e3|670e3|670e3|670e3, spharz: 410e3|520e3|560e3|670e3|670e3|670e3|670e3
set Phase transition widths = background:5e3|5e3|5e3|5e3|5e3|5e3|5e3, spharz: 5e3|5e3|5e3|5e3|5e3|5e3|5e3
set Phase transition temperatures = background:1662.0|1662.0|1662.0|1662.0|1662.0|1662.0|1662.0, spharz: 1662.0|1662.0|1662.0|1662.0|1662.0|1662.0|1662.0
set Phase transition Clapeyron slopes = background:4e6|4.1e6|4e6|-2e6|4e6|-3.1e6|1.3e6, spharz: 4e6|4.1e6|4e6|-2e6|4e6|-3.1e6|1.3e6
set Densities = background: 3300.0|3394.4|3442.1|3453.2|3617.6|3691.5|3774.7|3929.1, spharz: 3235.0|3372.3|3441.7|3441.7|3680.8|3717.8|3759.4|3836.6
set Heat capacities = 1250.0
set Thermal expansivities = 0.0
# The options to look up the dominant phase
set Use dominant phase for viscosity = true
set Material file names = material_table_temperature_pressure_small.txt, material_table_temperature_pressure_small.txt
set Phase transition phases = background:1, spharz:1
set Viscous flow law = diffusion
set Viscosity averaging scheme = harmonic
set Data directory = $ASPECT_SOURCE_DIR/data/material-model/entropy-table/pyrtable/
set Prefactors for diffusion creep = background:5e-21|5e-19, spharz:5e-25|5e-15
set Grain size exponents for diffusion creep = 0.0
set Activation energies for diffusion creep = 0.0
set Activation volumes for diffusion creep = 0.0
end
end

subsection Boundary velocity model
set Tangential velocity boundary indicators = top, bottom, left, right
end

subsection Postprocess
set List of postprocessors = visualization, material statistics

subsection Visualization
set Output format = gnuplot
set List of output variables = material properties
set Time between graphical output = 0e6

subsection Material properties
set List of material properties = density, thermal expansivity, specific heat, viscosity
end
end
end

14 changes: 14 additions & 0 deletions tests/visco_plastic_phases_discrete_comp/screen-output
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,14 @@

Number of active cells: 2,500 (on 1 levels)
Number of degrees of freedom: 43,405 (20,402+2,601+10,201+10,201)

*** Timestep 0: t=0 years, dt=0 years

Postprocessing:
Writing graphical output: output-visco_plastic_phases_discrete_comp/solution/solution-00000
Average density / Average viscosity / Total mass: 3428 kg/m^3, 7.6e+23 Pa s, 2.194e+15 kg

Termination requested by criterion: end time



13 changes: 13 additions & 0 deletions tests/visco_plastic_phases_discrete_comp/statistics
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,13 @@
# 1: Time step number
# 2: Time (years)
# 3: Time step size (years)
# 4: Number of mesh cells
# 5: Number of Stokes degrees of freedom
# 6: Number of temperature degrees of freedom
# 7: Number of degrees of freedom for all compositions
# 8: Number of nonlinear iterations
# 9: Visualization file name
# 10: Average density (kg/m^3)
# 11: Average viscosity (Pa s)
# 12: Total mass (kg)
0 0.000000000000e+00 0.000000000000e+00 2500 23003 10201 10201 0 output-visco_plastic_phases_discrete_comp/solution/solution-00000 3.42821554e+03 7.60000000e+23 2.19405794e+15
148 changes: 148 additions & 0 deletions tests/visco_plastic_phases_discrete_mixed.prm
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,148 @@
# This prm file is used to generate a phase diagram.
# and test the implementation of discrete phase functions
# in the visco-plastic material model. Specifically, this test
# checks the implementation for multiple compositions with
# the two compositions mixed by half and half.
# In the outputs, the material statistics are checked for
# correst viscosity values.

set Dimension = 2
set CFL number = 1.0
set End time = 0
set Adiabatic surface temperature = 1673.0
set Use years in output instead of seconds = true
set Nonlinear solver scheme = no Advection, no Stokes

# Model geometry
subsection Geometry model
set Model name = box

subsection Box
set X repetitions = 50
set Y repetitions = 50
set X extent = 800e3
set Y extent = 800e3
end
end

# Mesh refinement specifications
subsection Mesh refinement
set Initial adaptive refinement = 0
set Initial global refinement = 0
set Time steps between mesh refinement = 0
end

subsection Gravity model
set Model name = vertical

subsection Vertical
set Magnitude = 10.0
end
end

# Temperature boundary and initial conditions
# The initial temperature increases linearly from
# the left to the right boundary.
subsection Boundary temperature model
set Fixed temperature boundary indicators = left, right
set List of model names = box

subsection Box
set Left temperature = 273
set Right temperature = 2273
end
end

subsection Initial temperature model
set Model name = function

subsection Function
set Coordinate system = cartesian
set Variable names = x, y
set Function constants = XMAX=800e3, Tl=273.0, Tr=2273
set Function expression = Tl * (x - XMAX)/(-XMAX) + Tr * x / XMAX
end
end

# Comment the following 3 sections for steinberg model and pyrolitic lookup table
# Fields of composition
subsection Compositional fields
set Number of fields = 1
set Names of fields = spharz
set Compositional field methods = field
end

# Initial composition model
subsection Initial composition model
set List of model names = function

subsection Function
set Coordinate system = cartesian
set Function expression = 0.5
end
end

# Boundary composition model
subsection Boundary composition model
set List of model names = box
set Fixed composition boundary indicators = top, bottom, left, right

subsection Box
set Bottom composition = 0.5
set Left composition = 0.5
set Right composition = 0.5
set Top composition = 0.5
end
end

# Value for material model
# A set of parameters are assigned to the conventional
# phase function interface. A single phase with index 1
# is added to the background composition with a smaller
# viscosity by using the discrete phase transition.
subsection Material model
set Model name = visco plastic

subsection Visco Plastic
set Reference temperature = 273
set Maximum viscosity = 1e25
set Minimum viscosity = 1e13
set Phase transition depths = background:410e3|520e3|560e3|670e3|670e3|670e3|670e3, spharz: 410e3|520e3|560e3|670e3|670e3|670e3|670e3
set Phase transition widths = background:5e3|5e3|5e3|5e3|5e3|5e3|5e3, spharz: 5e3|5e3|5e3|5e3|5e3|5e3|5e3
set Phase transition temperatures = background:1662.0|1662.0|1662.0|1662.0|1662.0|1662.0|1662.0, spharz: 1662.0|1662.0|1662.0|1662.0|1662.0|1662.0|1662.0
set Phase transition Clapeyron slopes = background:4e6|4.1e6|4e6|-2e6|4e6|-3.1e6|1.3e6, spharz: 4e6|4.1e6|4e6|-2e6|4e6|-3.1e6|1.3e6
set Densities = background: 3300.0|3394.4|3442.1|3453.2|3617.6|3691.5|3774.7|3929.1, spharz: 3235.0|3372.3|3441.7|3441.7|3680.8|3717.8|3759.4|3836.6
set Heat capacities = 1250.0
set Thermal expansivities = 0.0
# The options to look up the dominant phase
set Use dominant phase for viscosity = true
set Material file names = material_table_temperature_pressure_small.txt, material_table_temperature_pressure_small.txt
set Phase transition phases = background:1, spharz:1
set Viscous flow law = diffusion
set Viscosity averaging scheme = harmonic
set Data directory = $ASPECT_SOURCE_DIR/data/material-model/entropy-table/pyrtable/
set Prefactors for diffusion creep = background:5e-21|5e-19, spharz:5e-25|5e-15
set Grain size exponents for diffusion creep = 0.0
set Activation energies for diffusion creep = 0.0
set Activation volumes for diffusion creep = 0.0
end
end

subsection Boundary velocity model
set Tangential velocity boundary indicators = top, bottom, left, right
end

subsection Postprocess
set List of postprocessors = visualization, material statistics

subsection Visualization
set Output format = gnuplot
set List of output variables = material properties
set Time between graphical output = 0e6

subsection Material properties
set List of material properties = density, thermal expansivity, specific heat, viscosity
end
end
end

Loading

0 comments on commit 35030e0

Please sign in to comment.