Skip to content

Commit

Permalink
fix freezing in melt simple material model and improve documentation
Browse files Browse the repository at this point in the history
  • Loading branch information
@jdannberg
jdannberg committed Mar 20, 2018
1 parent 4b5ef67 commit 04f5c1a
Show file tree
Hide file tree
Showing 9 changed files with 365 additions and 71 deletions.
6 changes: 6 additions & 0 deletions doc/modules/changes/20180316_jdannberg
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,6 @@
Fixed: The freezing of melt in the melt simple material model is
now more consistent: The freezing rate parameter determines the
freezing of all melt, and melt freezes if both the porosity and
the depletion are larger than the equilibrium melt fraction.
<br>
(Juliane Dannberg, 2018/03/16)
94 changes: 74 additions & 20 deletions source/material_model/melt_simple.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -258,7 +258,7 @@ namespace aspect
out.densities[i] = (reference_rho_s + delta_rho)
* temperature_dependence * std::exp(compressibility * (in.pressure[i] - this->get_surface_pressure()));

if (this->include_melt_transport())
if (this->include_melt_transport() && in.strain_rate.size())
{
AssertThrow(this->introspection().compositional_name_exists("peridotite"),
ExcMessage("Material model Melt simple only works if there is a "
Expand All @@ -271,44 +271,68 @@ namespace aspect

// calculate the melting rate as difference between the equilibrium melt fraction
// and the solution of the previous time step
double melting_rate = 0.0;
double porosity_change = 0.0;
if (fractional_melting)
{
// solidus is lowered by previous melting events (fractional melting)
const double solidus_change = (in.composition[i][peridotite_idx] - in.composition[i][porosity_idx]) * depletion_solidus_change;
const double eq_melt_fraction = melt_fraction(in.temperature[i] - solidus_change, this->get_adiabatic_conditions().pressure(in.position[i]));
melting_rate = eq_melt_fraction - old_porosity[i];
porosity_change = eq_melt_fraction - old_porosity[i];
}
else
// batch melting
melting_rate = melt_fraction(in.temperature[i], this->get_adiabatic_conditions().pressure(in.position[i]))
- std::max(maximum_melt_fractions[i], 0.0);
porosity_change = melt_fraction(in.temperature[i], this->get_adiabatic_conditions().pressure(in.position[i]))
- std::max(maximum_melt_fractions[i], 0.0);

// remove melt that gets close to the surface
if (this->get_geometry_model().depth(in.position[i]) < extraction_depth)
melting_rate = -old_porosity[i] * (in.position[i](1) - (this->get_geometry_model().maximal_depth() - extraction_depth))/extraction_depth;
porosity_change = -old_porosity[i] * (in.position[i](1) - (this->get_geometry_model().maximal_depth() - extraction_depth))/extraction_depth;

// freezing of melt below the solidus
{
const double freezing = freezing_rate * this->get_timestep() / year_in_seconds
* 0.5 * (melt_fraction(in.temperature[i], this->get_adiabatic_conditions().pressure(in.position[i])) - old_porosity[i]
- std::abs(melt_fraction(in.temperature[i], this->get_adiabatic_conditions().pressure(in.position[i])) - old_porosity[i]));
melting_rate += freezing;
// If the porosity is larger than the equilibrium melt fraction, melt should freeze again.
// Because we do not track the melt composition, we have to use a workaround here for freezing of melt:
// We reduce the porosity until either it reaches the equilibrium melt fraction, or the depletion
// (peridotite field), which decreases as melt freezes, reaches the same value as the equilibrium
// melt fraction, whatever happens earlier. An exception is when the melt fraction is zero; in this case
// all melt should freeze.
const double eq_melt_fraction = melt_fraction(in.temperature[i], this->get_adiabatic_conditions().pressure(in.position[i]));

// If the porosity change is not negative, there is no freezing, and the change in porosity
// is covered by the melting relation above.

// porosity reaches the equilibrium melt fraction:
const double porosity_change_wrt_melt_fraction = std::min(eq_melt_fraction - old_porosity[i] - porosity_change,0.0);

// depletion reaches the equilibrium melt fraction:
const double porosity_change_wrt_depletion = std::min((eq_melt_fraction - maximum_melt_fractions[i])
* (1.0 - old_porosity[i]) / (1.0 - maximum_melt_fractions[i]),0.0);
double freezing_amount = std::max(porosity_change_wrt_melt_fraction, porosity_change_wrt_depletion);

if (eq_melt_fraction == 0.0)
freezing_amount = - old_porosity[i];

porosity_change += freezing_amount;

// Adapt time scale of freezing with respect to melting.
if (porosity_change < 0 && this->get_parameters().use_operator_splitting)
porosity_change *= freezing_rate * melting_time_scale;
else if (porosity_change < 0)
porosity_change *= std::min(1.0, freezing_rate * this->get_timestep());
}

// do not allow negative porosity
if (old_porosity[i] + melting_rate < 0)
melting_rate = -old_porosity[i];
porosity_change = std::max(porosity_change, -old_porosity[i]);

// because depletion is a volume-based, and not a mass-based property that is advected,
// additional scaling factors on the right hand side apply
for (unsigned int c=0; c<in.composition[i].size(); ++c)
{
if (c == peridotite_idx && this->get_timestep_number() > 0 && (in.strain_rate.size()))
out.reaction_terms[i][c] = melting_rate * (1 - in.composition[i][peridotite_idx])
if (c == peridotite_idx && this->get_timestep_number() > 0)
out.reaction_terms[i][c] = porosity_change * (1 - in.composition[i][peridotite_idx])
/ (1 - in.composition[i][porosity_idx]);
else if (c == porosity_idx && this->get_timestep_number() > 0 && (in.strain_rate.size()))
out.reaction_terms[i][c] = melting_rate
else if (c == porosity_idx && this->get_timestep_number() > 0)
out.reaction_terms[i][c] = porosity_change
* out.densities[i] / this->get_timestep();
else
out.reaction_terms[i][c] = 0.0;
Expand All @@ -321,7 +345,7 @@ namespace aspect
if (c == peridotite_idx && this->get_timestep_number() > 0)
reaction_rate_out->reaction_rates[i][c] = out.reaction_terms[i][c] / melting_time_scale;
else if (c == porosity_idx && this->get_timestep_number() > 0)
reaction_rate_out->reaction_rates[i][c] = melting_rate / melting_time_scale;
reaction_rate_out->reaction_rates[i][c] = porosity_change / melting_time_scale;
else
reaction_rate_out->reaction_rates[i][c] = 0.0;
}
Expand Down Expand Up @@ -531,11 +555,34 @@ namespace aspect
"melting should be used (if false), assuming that the melt fraction only "
"depends on temperature and pressure, and how much melt has already been "
"generated at a given point, but not considering movement of melt in "
"the melting parameterization.");
"the melting parameterization."
"\n\n"
"Note that melt does not freeze unless the 'Freezing rate' parameter is set "
"to a value larger than 0.");
prm.declare_entry ("Freezing rate", "0.0",
Patterns::Double (0),
"Freezing rate of melt when in subsolidus regions."
"Units: $1/yr$.");
"Freezing rate of melt when in subsolidus regions. "
"If this parameter is set to a number larger than 0.0, it specifies the "
"fraction of melt that will freeze per year, as soon as the porosity "
"exceeds the equilibrium melt fraction, and the equilibrium melt fraction "
"falls below the depletion. In this case, melt will freeze according to the "
"given rate until one of those conditions is not fulfilled anymore. The "
"reasoning behind this is that there should not be more melt present than "
"the equilibrium melt fraction, as melt production decreases with increasing "
"depletion, but the freezing process of melt also reduces the depletion by "
"the same amount, and as soon as the depletion falls below the equilibrium "
"melt fraction, we expect that material should melt again (no matter how "
"much melt is present). This is quite a simplification and not a realistic "
"freezing parameterization, but without tracking the melt composition, there "
"is no way to compute freezing rates accurately. "
"If this parameter is set to zero, no freezing will occur. "
"Note that if operator splitting is used, freezing can never be faster than "
"determined by the ``Melting time scale for operator splitting''. In this case "
"the product of the ``Freezing rate'' and the ``Melting time scale for operator "
"splitting'' defines how fast freezing occurs with respect to melting (if the "
"product is 0.5, melting will occur twice as fast as freezing). "
"Units: 1/yr or 1/s, depending on the ``Use years "
"in output instead of seconds'' parameter.");
prm.declare_entry ("Melting time scale for operator splitting", "1e3",
Patterns::Double (0),
"In case the operator splitting scheme is used, the porosity field can not "
Expand Down Expand Up @@ -702,6 +749,7 @@ namespace aspect
if (this->get_parameters().convert_to_years == true)
{
melting_time_scale *= year_in_seconds;
freezing_rate /= year_in_seconds;
}

if (this->get_parameters().use_operator_splitting)
Expand All @@ -714,6 +762,12 @@ namespace aspect
+ Utilities::to_string(melting_time_scale) + "."));
AssertThrow(melting_time_scale > 0,
ExcMessage("The Melting time scale for operator splitting must be larger than 0!"));
AssertThrow(freezing_rate * this->get_parameters().reaction_time_step <= 1.0,
ExcMessage("The reaction time step " + Utilities::to_string(this->get_parameters().reaction_time_step)
+ " in the operator splitting scheme is too large to compute freezing rates! "
"You have to choose it in such a way that it is smaller than the inverse of the "
"'Freezing rate' chosen in the material model, which is currently "
+ Utilities::to_string(1.0/freezing_rate) + "."));
}

A1 = prm.get_double ("A1");
Expand Down
1 change: 1 addition & 0 deletions tests/melting_rate_operator_splitting.prm
View file
Original file line number Diff line number Diff line change
Expand Up @@ -112,6 +112,7 @@ subsection Material model
set Exponential melt weakening factor = 0.0
set Melt extraction depth = 0.0
set Melting time scale for operator splitting = 1.111111111e3
set Freezing rate = 9e-4

# simplifiy the melting model so that porosity depends on temperature linearly
set A1 = 726.85 # this is in °C --> 1000 K
Expand Down
124 changes: 124 additions & 0 deletions tests/rising_melt_blob_freezing.prm
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,124 @@
#########################################################
# This is a test for the melting rate functionality,
# in particular the freezing of melt.
# It uses the melt simple model.

# The initial porosity has the form of a Gaussian function
# there is a fixed (constant) freezing rate.

# The amplitude of the Gaussian is 0.2, and the melting time
# scale is chosen in such a way that the fraction of melt
# should be cut in half every year (visible in the maximum
# value of the porosity field).

############### Global parameters

set Dimension = 2
set Start time = 0
set End time = 2
set Use years in output instead of seconds = true
set Use direct solver for Stokes system = false
set Adiabatic surface temperature = 1600
set Nonlinear solver scheme = iterated IMPES
set Max nonlinear iterations = 50
set Maximum time step = 1

set Pressure normalization = surface
set Surface pressure = 0

subsection Melt settings
set Melt transport threshold = 0.0
end

subsection Material model

set Model name = melt simple
subsection Melt simple
set Melt extraction depth = 0.0
set Freezing rate = 1.0
set Melting time scale for operator splitting = 1.442695041
end
end

set Use operator splitting = true
subsection Solver parameters
subsection Operator splitting parameters
set Reaction time step = 0.005
set Reaction time steps per advection step = 1
end
end

############### Parameters describing the model
subsection Geometry model
set Model name = box

subsection Box
set X extent = 100000
set Y extent = 100000
end
end

subsection Model settings
set Tangential velocity boundary indicators = 0, 1, 2, 3
set Include melt transport = true
end

subsection Compositional fields
set Number of fields = 2
set Names of fields = porosity, peridotite
end

subsection Gravity model
set Model name = vertical
subsection Vertical
set Magnitude = 10.0
end
end

subsection Initial temperature model
set Model name = function
subsection Function
set Function expression = 293
end
end

subsection Initial composition model
set Model name = function

subsection Function
set Variable names = x,y
set Function constants = pi=3.1415926,x0=50000,y0=50000,c=10000
set Function expression = 0.2 * exp(-((x-x0)*(x-x0)+(y-y0)*(y-y0))/(2*c*c)); 0.0 * exp(-((x-x0)*(x-x0)+(y-y0)*(y-y0))/(2*c*c))
end
end


# The final part of this input file describes how many times the
# mesh is refined and what to do with the solution once computed
subsection Mesh refinement
set Initial adaptive refinement = 0
set Initial global refinement = 5
set Time steps between mesh refinement = 0
end

subsection Postprocess

set List of postprocessors = composition statistics, visualization

subsection Visualization

set List of output variables = density, viscosity, thermal expansivity, melt material properties

subsection Melt material properties
set List of properties = permeability, fluid density, compaction viscosity, fluid viscosity
end

set Number of grouped files = 0

set Output format = vtu

set Time between graphical output = 0

end

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

Number of active cells: 1,024 (on 6 levels)
Number of degrees of freedom: 32,842 (8,450+2,178+8,450+1,089+4,225+4,225+4,225)

*** Timestep 0: t=0 years
Solving temperature system... 0 iterations.
Solving porosity system ... 0 iterations.
Skipping peridotite composition solve because RHS is zero.
Rebuilding Stokes preconditioner...
Solving Stokes system... 200+665 iterations.
Solving for u_f in 12 iterations.
Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.76497e-16, 1.44358e-16, 0, 0.168633
Relative nonlinear residual (total system) after nonlinear iteration 1: 0.168633

Solving temperature system... 0 iterations.
Solving porosity system ... 0 iterations.
Skipping peridotite composition solve because RHS is zero.
Rebuilding Stokes preconditioner...
Solving Stokes system... 0+0 iterations.
Solving for u_f in 12 iterations.
Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.76497e-16, 1.44358e-16, 0, 9.99026e-08
Relative nonlinear residual (total system) after nonlinear iteration 2: 9.99026e-08


Postprocessing:
Compositions min/max/mass: 2.778e-12/0.2/1.257e+08 // 0/0/0
Writing graphical output: output-rising_melt_blob_freezing/solution/solution-00000

*** Timestep 1: t=1 years
Solving composition reactions in 200 substep(s).
Solving temperature system... 0 iterations.
Solving porosity system ... 6 iterations.
Solving peridotite system ... 5 iterations.
Rebuilding Stokes preconditioner...
Solving Stokes system... 200+726 iterations.
Solving for u_f in 12 iterations.
Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.79631e-16, 3.76619e-06, 2.26327e-07, 0.00840578
Relative nonlinear residual (total system) after nonlinear iteration 1: 0.00840578

Solving temperature system... 0 iterations.
Solving porosity system ... 6 iterations.
Solving peridotite system ... 4 iterations.
Rebuilding Stokes preconditioner...
Solving Stokes system... 10+0 iterations.
Solving for u_f in 12 iterations.
Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.73728e-16, 1.88328e-06, 1.38447e-07, 1.08144e-07
Relative nonlinear residual (total system) after nonlinear iteration 2: 1.88328e-06


Postprocessing:
Compositions min/max/mass: -7.531e-11/0.09988/6.276e+07 // -0.1252/-1.39e-12/-7.019e+07
Writing graphical output: output-rising_melt_blob_freezing/solution/solution-00001

*** Timestep 2: t=2 years
Solving composition reactions in 200 substep(s).
Solving temperature system... 0 iterations.
Solving porosity system ... 6 iterations.
Solving peridotite system ... 5 iterations.
Rebuilding Stokes preconditioner...
Solving Stokes system... 200+601 iterations.
Solving for u_f in 12 iterations.
Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.77111e-16, 3.79756e-06, 7.12047e-08, 0.0313658
Relative nonlinear residual (total system) after nonlinear iteration 1: 0.0313658

Solving temperature system... 0 iterations.
Solving porosity system ... 6 iterations.
Solving peridotite system ... 3 iterations.
Rebuilding Stokes preconditioner...
Solving Stokes system... 8+0 iterations.
Solving for u_f in 12 iterations.
Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.75463e-16, 1.4052e-06, 2.21463e-08, 1.42281e-07
Relative nonlinear residual (total system) after nonlinear iteration 2: 1.4052e-06


Postprocessing:
Compositions min/max/mass: -8.803e-11/0.04988/3.134e+07 // -0.1877/3.631e-11/-1.052e+08
Writing graphical output: output-rising_melt_blob_freezing/solution/solution-00002

Termination requested by criterion: end time


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

0 comments on commit 04f5c1a

Please sign in to comment.