From 291e56f59a2605b937f26874143c3b5995414e4f Mon Sep 17 00:00:00 2001 From: Juliane Dannberg Date: Wed, 4 Apr 2018 18:31:26 -0700 Subject: [PATCH] address comments --- doc/update_prm_files_to_2.0.0.sed | 3 + include/aspect/compat.h | 4 +- include/aspect/melt.h | 152 ++++--- instant_melt.prm | 228 ---------- source/particle/property/melt_particle.cc | 4 +- source/postprocess/visualization/melt.cc | 5 +- source/simulator/core.cc | 4 +- source/simulator/melt.cc | 408 ++++++++++-------- source/simulator/parameters.cc | 16 +- source/simulator/solver.cc | 5 +- tests/compaction_length_refinement.prm | 4 - tests/composition_active_with_melt.prm | 1 - tests/free_surface_blob_nonzero_melt.prm | 4 - tests/global_melt.prm | 4 - tests/global_melt_parallel.prm | 4 - tests/initial_porosity.prm | 4 - tests/melt_and_traction.prm | 4 - tests/melt_force_vector.prm | 1 - tests/melt_material_4.prm | 4 - tests/melt_track.prm | 4 - tests/melt_transport.prm | 4 - tests/melt_velocity_boundary_conditions.prm | 4 - .../solver_history.txt | 12 +- tests/particle_melt_advection.prm | 6 +- tests/particle_property_composition.prm | 4 - tests/rising_melt_blob.prm | 25 +- tests/rising_melt_blob_freezing.prm | 4 - tests/shear_bands.prm | 4 - tests/shear_heating_with_melt.prm | 7 +- tests/solitary_wave.prm | 4 - tests/visco_plastic_derivatives.prm | 2 +- 31 files changed, 355 insertions(+), 584 deletions(-) delete mode 100644 instant_melt.prm diff --git a/doc/update_prm_files_to_2.0.0.sed b/doc/update_prm_files_to_2.0.0.sed index a25e8735a2e..2c42e91625a 100644 --- a/doc/update_prm_files_to_2.0.0.sed +++ b/doc/update_prm_files_to_2.0.0.sed @@ -66,3 +66,6 @@ s/boussinesq/Boussinesq/g # Rename Dynamic topography subsection s/subsection Dynamic Topography/subsection Dynamic topography/g + +# Remove all instances of `melt transport threshold` +/Melt transport threshold/d diff --git a/include/aspect/compat.h b/include/aspect/compat.h index 52e6e1ac295..e4d0060cd2e 100644 --- a/include/aspect/compat.h +++ b/include/aspect/compat.h @@ -59,8 +59,8 @@ namespace aspect { dealii::SolverControl::State return_value = dealii::SolverControl::check(step, check_value); - if (history_data_enabled) - history_data.resize(maxsteps+1); + if (step == 0) + history_data.resize(history_data.size()+1); return return_value; } diff --git a/include/aspect/melt.h b/include/aspect/melt.h index 0572af65ccd..470efeaa8ea 100644 --- a/include/aspect/melt.h +++ b/include/aspect/melt.h @@ -132,20 +132,23 @@ namespace aspect * Returns the cell-averaged and cut-off value of p_c_scale, * the factor we use to rescale the compaction pressure and to * decide if a cell is a melt cell. - * The last input argument consider_is_melt_cell determines if - * this computation takes into account if a cell is a "melt cell" - * (cells where we solve the melt transport equations, as - * indicated by the entries stored in the is_melt_cell vector of - * the melt handler) and return a value of zero if the cell is - * not a melt cell (if true), or whether the computation should - * disregard the information about which cells are melt cells, - * which is required when we want to update the is_melt_cell - * vector and find out which cells should be melt cells (if false). + * The last input argument @p consider_is_melt_cell determines if + * this computation takes into account if a cell is a "melt cell". + * Melt cells are cells where we solve the melt transport equations, + * as indicated by the entries stored in the is_melt_cell vector of + * the melt handler. In case @p consider_is_melt_cell is set to true, + * this function returns a value of zero if the cell is not a melt cell. + * If @p consider_is_melt_cell is set to false the computation + * disregards the information about which cells are melt cells, + * and computes p_c_scale from the cell-averaged Darcy coefficient + * for all cells. This is needed for example when we want to update + * the is_melt_cell vector and need to find out which cells should be + * marked as melt cells. */ double p_c_scale (const MaterialModel::MaterialModelInputs &inputs, const MaterialModel::MaterialModelOutputs &outputs, const MeltHandler &handler, - bool consider_is_melt_cell) const; + const bool consider_is_melt_cell) const; }; @@ -271,27 +274,20 @@ namespace aspect } - /** - * Class that contains all runtime parameters and other helper functions - * related to melt transport. A global instance can be retrieved with - * SimulatorAccess::get_melt_handler(), but keep in mind that it only - * exists if parameters.include_melt_transport is true. - */ - template - class MeltHandler: public SimulatorAccess + namespace Melt { - public: - MeltHandler(ParameterHandler &prm); - + template + struct Parameters + { /** * Declare additional parameters that are needed in models with - * melt transport (including the fluid pressure boundary conditions). + * melt transport. */ static void declare_parameters (ParameterHandler &prm); /** * Parse additional parameters that are needed in models with - * melt transport (including the fluid pressure boundary conditions). + * melt transport. * * This has to be called before edit_finite_element_variables, * so that the finite elements that are used for the additional melt @@ -300,6 +296,58 @@ namespace aspect */ void parse_parameters (ParameterHandler &prm); + /** + * The factor by how much the Darcy coefficient K_D in a cell can be smaller than + * the reference Darcy coefficient for this cell still to be considered a melt cell + * (for which the melt transport equations are solved). If the Darcy coefficient + * is smaller than the product of this value and the reference Dracy coefficient, + * the cell is not considered a melt cell and the Stokes system without melt + * transport is solved instead. In practice, this means that all terms in the + * assembly related to the migration of melt are set to zero, and the compaction + * pressure degrees of freedom are constrained to be zero. + */ + double K_D_variation_threshold; + + /** + * Whether to use a porosity weighted average of the melt and solid velocity + * to advect heat in the temperature equation or not. If this is set to true, + * additional terms are assembled on the left-hand side of the temperature + * advection equation in models with melt migration. + * If this is set to false, only the solid velocity is used (as in models + * without melt migration). + */ + bool heat_advection_by_melt; + + /** + * Whether to use a discontinuous element for the compaction pressure or not. + */ + bool use_discontinuous_p_c; + + /** + * Whether to cell-wise average the material properties that are used to + * compute the melt velocity or not. Note that the melt velocity is computed + * as the sum of the solid velocity and the phase separation flux (difference + * between melt and solid velocity). If this parameter is set to true, + * material properties in the computation of the phase separation flux will + * be averaged cell-wise. + */ + bool average_melt_velocity; + }; + } + + + /** + * Class that contains all runtime parameters and other helper functions + * related to melt transport. A global instance can be retrieved with + * SimulatorAccess::get_melt_handler(), but keep in mind that it only + * exists if parameters.include_melt_transport is true. + */ + template + class MeltHandler: public SimulatorAccess + { + public: + MeltHandler(ParameterHandler &prm); + /** * Create an additional material model output object that contains * the additional output variables needed in simulation with melt transport, @@ -323,6 +371,15 @@ namespace aspect */ void set_assemblers (Assemblers::Manager &assemblers) const; + + /** + * Initialize function. This is mainly to check that the melt transport + * parameters chosen in the input file are consistent with the rest of + * the options. We can not do this in the parse_parameters function, + * as we do not have simulator access at that point. + */ + void initialize() const; + /** * Setup SimulatorAccess for the plugins related to melt transport. */ @@ -380,58 +437,34 @@ namespace aspect bool is_melt_cell(const typename DoFHandler::active_cell_iterator &cell) const; /** - * Given the Darcy coefficient as computed by the material model, limit the - * coefficient to a minimum value based on the reference Darcy coefficient - * in melt cells, and set it to zero in cells that are not melt cells, and - * return this value. + * Given the Darcy coefficient @p K_D as computed by the material model, limit the + * coefficient to a minimum value (computed as the K_D variation threshold given in + * the input file times the reference Darcy coefficient) in melt cells and return + * this value. If @p is_melt_cell is false, return zero. */ double limited_darcy_coefficient(const double K_D, const bool is_melt_cell) const; /** - * The porosity limit for melt migration. For smaller porosities, the equations - * reduce to the Stokes equations and neglect melt transport. In practice, this - * means that all terms in the assembly related to the migration of melt are set - * to zero for porosities smaller than this threshold. - * This does not include the compaction term $p_c/\xi$, which is necessary for the - * solvability of the linear system, but does not influence the solution variables - * of the Stokes problem (in the absence of porosity). + * Return a pointer to the boundary fluid pressure. */ - double melt_transport_threshold; + const BoundaryFluidPressure::Interface & + get_boundary_fluid_pressure () const; /** - * Whether to use a porosity weighted average of the melt and solid velocity - * to advect heat in the temperature equation or not. If this is set to true, - * additional terms are assembled on the left-hand side of the temperature - * advection equation in models with melt migration. - * If this is set to false, only the solid velocity is used (as in models - * without melt migration). + * The object that stores the run-time parameters that control the how the + * melt transport equations are solved. */ - bool heat_advection_by_melt; + Melt::Parameters melt_parameters; + private: /** * Store a pointer to the fluid pressure boundary plugin, so that the * initialization can be done together with the other objects related to melt * transport. */ - std_cxx11::unique_ptr > boundary_fluid_pressure; + const std_cxx11::unique_ptr > boundary_fluid_pressure; - /** - * Whether to use a discontinuous element for the compaction pressure or not. - */ - bool use_discontinuous_p_c; - - /** - * Whether to cell-wise average the material properties that are used to - * compute the melt velocity or not. Note that the melt velocity is computed - * as the sum of the solid velocity and the phase separation flux (difference - * between melt and solid velocity). If this parameter is set to true, - * material properties in the computation of the phase separation flux will - * be averaged cell-wise. - */ - bool average_melt_velocity; - - private: /** * is_melt_cell_vector[cell->active_cell_index()] says whether we want to * solve the melt transport equations (as opposed to the Stokes equations without @@ -448,6 +481,7 @@ namespace aspect * we have computed this solution. */ ConstraintMatrix current_constraints; + }; } diff --git a/instant_melt.prm b/instant_melt.prm deleted file mode 100644 index e62162543e8..00000000000 --- a/instant_melt.prm +++ /dev/null @@ -1,228 +0,0 @@ -# Listing of Parameters -# --------------------- -# Global-scale mantle convection model with melt migration - -set Adiabatic surface temperature = 1600 # default: 0 -set CFL number = 1.0 -set Maximum time step = 1e6 -set Composition solver tolerance = 1e-14 -set Temperature solver tolerance = 1e-14 -set Nonlinear solver scheme = iterated IMPES -set Output directory = output-instant5 -set Max nonlinear iterations = 20 -set Linear solver tolerance = 1e-8 -set Nonlinear solver tolerance = 1e-5 - -# The number of space dimensions you want to run this program in. -set Dimension = 2 - -# The end time of the simulation. Units: years if the 'Use years in output -# instead of seconds' parameter is set; seconds otherwise. -# This end time is chosen in such a way that the solitary wave travels -# approximately 5 times its wavelength during the model time. -set End time = 2000 #70863 #9e8 - -set Pressure normalization = surface -set Surface pressure = 0 -set Resume computation = false -set Start time = 0 - -set Use years in output instead of seconds = true -set Use direct solver for Stokes system = false -set Number of cheap Stokes solver steps = 5000 #0 - -subsection Discretization - set Stokes velocity polynomial degree = 2 - set Composition polynomial degree = 1 - subsection Stabilization parameters - set beta = 0.2 - end -end - -subsection Compositional fields - set Number of fields = 2 - set Names of fields = porosity, peridotite -end - - -subsection Boundary temperature model - set Model name = initial temperature - - subsection Initial temperature - set Minimal temperature = 293 # default: 6000 - set Maximal temperature = 3700 # default: 0 - end -end - -subsection Boundary composition model - set Model name = initial composition -end - -subsection Boundary velocity model - subsection Function - set Function constants = b=100000, c=200000 - set Variable names = x,y - set Function expression = 0.0; -0.024995 + 0.1 * exp(-((x-b)*(x-b)+y*y)/(2*c*c)) - end -end - -subsection Geometry model - set Model name = box - - subsection Box - set X extent = 8700000 - set Y extent = 2900000 -# set X periodic = true - set X repetitions = 3 - end - -end - -subsection Gravity model - set Model name = vertical - - subsection Vertical - set Magnitude = 9.81 - end -end - - -subsection Initial temperature model - set Model name = adiabatic - subsection Adiabatic - set Age bottom boundary layer = 5e8 - set Age top boundary layer = 3e8 - set Amplitude = 50 - set Position = center - set Radius = 350000 - - subsection Function - set Function expression = 0;0 - end - end - - subsection Harmonic perturbation - set Magnitude = 50 - end -end - -subsection Initial composition model - set Model name = function - subsection Function - set Function constants = pi=3.1415926,a = 0.10, b = 2500000, c = 150000, d=1450000 - set Function expression = a * exp(-((y-b)*(y-b)+(0.5*(x-d))*(0.5*(x-d)))/(2*c*c)); 0.0*a * exp(-((y-b)*(y-b)+(0.2*(x-d))*(0.2*(x-d)))/(2*c*c)) - set Variable names = x,y - end -end - - -subsection Material model - - set Model name = melt global - subsection Melt global - set Thermal conductivity = 4.7 - set Reference solid density = 3400 - set Reference melt density = 3000 - set Thermal expansion coefficient = 2e-5 - set Reference permeability = 1e-8 - set Reference shear viscosity = 5e21 - set Reference bulk viscosity = 1e19 - set Exponential melt weakening factor = 10 - set Thermal viscosity exponent = 7 - set Thermal bulk viscosity exponent = 7 - set Reference temperature = 1600 - set Solid compressibility = 4.2e-12 - set Melt compressibility = 1.25e-11 - set Melt bulk modulus derivative = 0.0 - set Reference melt viscosity = 10 - - # change the density contrast of depleted material (in kg/m^3) - set Depletion density change = -200.0 # -100.0 # 0.0 - end -end - - -subsection Mesh refinement - set Coarsening fraction = 0.05 - set Refinement fraction = 0.8 - - set Initial adaptive refinement = 0 #2 # default: 2 - set Initial global refinement = 5 # default: - set Strategy = composition threshold, minimum refinement function #, nonadiabatic temperature - set Time steps between mesh refinement = 0 #5 - - subsection Minimum refinement function - set Coordinate system = depth - set Function expression = if (depth>1500000,5,5) - set Variable names = depth,phi - end - - subsection Composition threshold - set Compositional field thresholds = 1e-4,1.0 - end -end - - -subsection Boundary fluid pressure model - set Plugin name = density - subsection Density - set Density formulation = fluid density - end -end - -subsection Heating model - set List of model names = adiabatic heating -end - -subsection Model settings - set Fixed temperature boundary indicators = 2,3 - set Fixed composition boundary indicators = #2,3 - set Prescribed velocity boundary indicators = - - set Tangential velocity boundary indicators = 0,1,2,3 - set Zero velocity boundary indicators = -# set Remove nullspace = net x translation - - set Include melt transport = true -end - -subsection Melt settings - set Melt transport threshold = 1e-4 -end - - -subsection Postprocess - - set List of postprocessors = visualization,composition statistics,velocity statistics, temperature statistics, depth average - - subsection Visualization - - set List of output variables = melt material properties, material properties, nonadiabatic temperature, melt fraction, strain rate - - subsection Material properties - set List of material properties = density, viscosity, thermal expansivity, reaction terms - end - - subsection Melt material properties - set List of properties = fluid density, permeability, fluid viscosity, compaction viscosity - end - - set Number of grouped files = 0 - set Interpolate output = false - set Output format = vtu - set Time between graphical output = 6e6 - set Interpolate output = true - end - - subsection Depth average - set Number of zones = 12 - set Time between graphical output = 6e6 - end - -end - -subsection Checkpointing - set Time between checkpoint = 1700 -end - - diff --git a/source/particle/property/melt_particle.cc b/source/particle/property/melt_particle.cc index 4aad2b208e6..e3452137c09 100644 --- a/source/particle/property/melt_particle.cc +++ b/source/particle/property/melt_particle.cc @@ -111,7 +111,9 @@ namespace aspect prm.enter_subsection("Particles"); { prm.enter_subsection("Melt particle"); - threshold_for_melt_presence = prm.get_double ("Threshold for melt presence"); + { + threshold_for_melt_presence = prm.get_double ("Threshold for melt presence"); + } prm.leave_subsection(); } prm.leave_subsection(); diff --git a/source/postprocess/visualization/melt.cc b/source/postprocess/visualization/melt.cc index a4962bc8601..e40ba3919b8 100644 --- a/source/postprocess/visualization/melt.cc +++ b/source/postprocess/visualization/melt.cc @@ -115,7 +115,10 @@ namespace aspect AssertThrow(melt_outputs != NULL, ExcMessage("Need MeltOutputs from the material model for computing the melt properties.")); - const double p_c_scale = dynamic_cast*>(&this->get_material_model())->p_c_scale(in, out, this->get_melt_handler(), true); + const double p_c_scale = dynamic_cast*>(&this->get_material_model())->p_c_scale(in, + out, + this->get_melt_handler(), + true); for (unsigned int q=0; qinitialize_simulator (*this); + melt_handler->initialize(); } // If the solver type is a Newton type of solver, we need to set make sure diff --git a/source/simulator/melt.cc b/source/simulator/melt.cc index abc0fb85591..9e7d7c1ad49 100644 --- a/source/simulator/melt.cc +++ b/source/simulator/melt.cc @@ -66,12 +66,15 @@ namespace aspect double MeltInterface::p_c_scale (const MaterialModel::MaterialModelInputs &inputs, const MaterialModel::MaterialModelOutputs &outputs, - const MeltHandler &handler, - bool consider_is_melt_cell) const + const MeltHandler &melt_handler, + const bool consider_is_melt_cell) const { - Assert(inputs.current_cell.state() == IteratorState::valid, ExcMessage("sorry, I need a cell.")); + Assert(inputs.current_cell.state() == IteratorState::valid, + ExcMessage("You are trying to access the property of a cell (if it is " + "a melt cell or not), but the cell iterator in the provided material " + "model inputs does not point to a cell.")); - if (consider_is_melt_cell && !handler.is_melt_cell(inputs.current_cell)) + if (consider_is_melt_cell && !melt_handler.is_melt_cell(inputs.current_cell)) return 0.0; const MaterialModel::MeltOutputs *melt_outputs = outputs.template get_additional_output >(); @@ -92,10 +95,11 @@ namespace aspect // The same threshold is used when computing which cells are melt cells (the else branch), // with the difference that we return a p_c_scale of zero to indicate that the cell is not // a melt cell if the maximum Darcy coefficient of the cell is below the threshold. + const double K_D_threshold = melt_handler.melt_parameters.K_D_variation_threshold; if (consider_is_melt_cell) - K_D = std::max(K_D, 1e-3 * ref_K_D); + K_D = std::max(K_D, K_D_threshold * ref_K_D); else - K_D = (max_K_D < 1e-3 * ref_K_D) ? 0.0 : K_D; + K_D = (max_K_D < K_D_threshold * ref_K_D) ? 0.0 : K_D; // If the reference permeability is set to zero, there is no melt transport in the whole model and we return zero. return (ref_K_D > 0 ? std::sqrt(K_D / ref_K_D) : 0.0); @@ -165,7 +169,10 @@ namespace aspect const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points; const double pressure_scaling = this->get_pressure_scaling(); - const double p_c_scale = dynamic_cast*>(&this->get_material_model())->p_c_scale(scratch.material_model_inputs, scratch.material_model_outputs, this->get_melt_handler(), true); + const double p_c_scale = dynamic_cast*>(&this->get_material_model())->p_c_scale(scratch.material_model_inputs, + scratch.material_model_outputs, + this->get_melt_handler(), + true); MaterialModel::MeltOutputs *melt_outputs = scratch.material_model_outputs.template get_additional_output >(); @@ -212,7 +219,8 @@ namespace aspect const double eta = scratch.material_model_outputs.viscosities[q]; const double one_over_eta = 1. / eta; const double eta_two_thirds = scratch.material_model_outputs.viscosities[q] * 2.0 / 3.0; - const double K_D = this->get_melt_handler().limited_darcy_coefficient(melt_outputs->permeabilities[q] / melt_outputs->fluid_viscosities[q], p_c_scale > 0); + const double K_D = this->get_melt_handler().limited_darcy_coefficient(melt_outputs->permeabilities[q] / melt_outputs->fluid_viscosities[q], + p_c_scale > 0); /* - R = 1/eta M_p + K_D L_p for p @@ -257,21 +265,21 @@ namespace aspect (scratch.phi_p_c[i] * scratch.phi_p_c[j]) ) * JxW; + // add S between p_c and p_f - if (true) - data.local_matrix(i,j) += - ( - (p_c_scale * one_over_eta * - pressure_scaling * - pressure_scaling) - * scratch.phi_p[i] * scratch.phi_p_c[j] - + - (p_c_scale * one_over_eta * - pressure_scaling * - pressure_scaling) - * scratch.phi_p_c[i] * scratch.phi_p[j] - ) - * JxW; + data.local_matrix(i,j) += + ( + (p_c_scale * one_over_eta * + pressure_scaling * + pressure_scaling) + * scratch.phi_p[i] * scratch.phi_p_c[j] + + + (p_c_scale * one_over_eta * + pressure_scaling * + pressure_scaling) + * scratch.phi_p_c[i] * scratch.phi_p[j] + ) + * JxW; } } } @@ -355,7 +363,10 @@ namespace aspect Assert(dynamic_cast*>(&this->get_material_model()) != NULL, ExcMessage("Error: The current material model needs to be derived from MeltInterface to use melt transport.")); - const double p_c_scale = dynamic_cast*>(&this->get_material_model())->p_c_scale(scratch.material_model_inputs, scratch.material_model_outputs, this->get_melt_handler(), true); + const double p_c_scale = dynamic_cast*>(&this->get_material_model())->p_c_scale(scratch.material_model_inputs, + scratch.material_model_outputs, + this->get_melt_handler(), + true); const FEValuesExtractors::Scalar extractor_pressure = introspection.variable("fluid pressure").extractor_scalar(); const FEValuesExtractors::Scalar ex_p_c = introspection.variable("compaction pressure").extractor_scalar(); @@ -433,7 +444,8 @@ namespace aspect const unsigned int porosity_index = introspection.compositional_index_for_name("porosity"); const double porosity = std::max(scratch.material_model_inputs.composition[q][porosity_index],0.000); - const double K_D = this->get_melt_handler().limited_darcy_coefficient(melt_outputs->permeabilities[q] / melt_outputs->fluid_viscosities[q], p_c_scale > 0); + const double K_D = this->get_melt_handler().limited_darcy_coefficient(melt_outputs->permeabilities[q] / melt_outputs->fluid_viscosities[q], + p_c_scale > 0); const double viscosity_c = melt_outputs->compaction_viscosities[q]; const Tensor<1,dim> density_gradient_f = melt_outputs->fluid_density_gradients[q]; const double density_f = melt_outputs->fluid_densities[q]; @@ -535,14 +547,14 @@ namespace aspect const FEValuesExtractors::Scalar ex_p_f = introspection.variable("fluid pressure").extractor_scalar(); const unsigned int p_f_component_index = introspection.variable("fluid pressure").first_component_index; const unsigned int p_c_component_index = introspection.variable("compaction pressure").first_component_index; - const double p_c_scale = dynamic_cast*>(&this->get_material_model())->p_c_scale(scratch.material_model_inputs, scratch.material_model_outputs, this->get_melt_handler(), true); + const bool is_melt_cell = this->get_melt_handler().is_melt_cell(scratch.material_model_inputs.current_cell); const typename DoFHandler::face_iterator face = scratch.face_material_model_inputs.current_cell->face(scratch.face_number); MaterialModel::MeltOutputs *melt_outputs = scratch.face_material_model_outputs.template get_additional_output >(); std::vector grad_p_f(n_face_q_points); - this->get_melt_handler().boundary_fluid_pressure->fluid_pressure_gradient( + this->get_melt_handler().get_boundary_fluid_pressure().fluid_pressure_gradient( face->boundary_id(), scratch.face_material_model_inputs, scratch.face_material_model_outputs, @@ -558,7 +570,8 @@ namespace aspect const Tensor<1,dim> gravity = this->get_gravity_model().gravity_vector (scratch.face_finite_element_values.quadrature_point(q)); const double density_f = melt_outputs->fluid_densities[q]; - const double K_D = this->get_melt_handler().limited_darcy_coefficient(melt_outputs->permeabilities[q] / melt_outputs->fluid_viscosities[q], p_c_scale > 0); + const double K_D = this->get_melt_handler().limited_darcy_coefficient(melt_outputs->permeabilities[q] / melt_outputs->fluid_viscosities[q], + is_melt_cell); for (unsigned int i=0, i_stokes=0; i_stokesis_temperature() && this->get_melt_handler().is_melt_cell(scratch.cell) - && this->get_melt_handler().heat_advection_by_melt) + && this->get_melt_handler().melt_parameters.heat_advection_by_melt) { density_c_P_solid = (1.0 - porosity) * scratch.material_model_outputs.densities[q] * scratch.material_model_outputs.specific_heat[q]; density_c_P_melt = porosity * melt_outputs->fluid_densities[q] * scratch.material_model_outputs.specific_heat[q]; @@ -1150,7 +1163,7 @@ namespace aspect this->get_material_model().evaluate(in, out); - const double p_c_scale = dynamic_cast*>(&this->get_material_model())->p_c_scale(in,out,this->get_melt_handler(), true); + const bool is_melt_cell = this->get_melt_handler().is_melt_cell(in.current_cell); MaterialModel::MeltOutputs *melt_outputs = out.template get_additional_output >(); Assert(melt_outputs != NULL, ExcMessage("Need MeltOutputs from the material model for computing the melt variables.")); @@ -1160,7 +1173,7 @@ namespace aspect double K_D_over_phi = 1.0; - if (average_melt_velocity) + if (melt_parameters.average_melt_velocity) { // average the K_D and the porosity cell-wise (geometric mean) for (unsigned int q=0; q 0) + if (is_melt_cell) { const Tensor<1,dim> gravity = this->get_gravity_model().gravity_vector(in.position[q]); @@ -1287,15 +1300,18 @@ namespace aspect solution, p_c_values); fe_values[this->introspection().variable("fluid pressure").extractor_scalar()].get_function_values ( solution, p_f_values); - this->compute_material_model_input_values (solution, - fe_values, - cell, - true, - in); + in.reinit(fe_values, + cell, + this->introspection(), + solution, + false); this->get_material_model().evaluate(in, out); - const double p_c_scale = dynamic_cast*>(&this->get_material_model())->p_c_scale(in,out,this->get_melt_handler(), true); + const double p_c_scale = dynamic_cast*>(&this->get_material_model())->p_c_scale(in, + out, + this->get_melt_handler(), + true); for (unsigned int j=0; jget_fe().base_element(this->introspection().base_elements.pressure).dofs_per_cell; ++j) { @@ -1310,7 +1326,7 @@ namespace aspect const double phi = std::max(0.0, porosity_values[j]); double p = p_f_values[j]; - if (p_c_scale > 0) + if (p_c_scale > 0 && (1.0-phi) > std::numeric_limits::min()) p = (p_c_scale*p_c_values[j] - (phi-1.0) * p_f_values[j]) / (1.0-phi); distributed_vector(local_dof_indices[pressure_idx]) = p; @@ -1335,15 +1351,19 @@ namespace aspect namespace { + // This is a scratch object for the setting the compaction pressure constraints + // in cells without melt (where we do not solve the melt transport equations, + // so we set the compaction pressure to zero). template - struct PcAssembleData + struct PcConstraintsAssembleData { - PcAssembleData (const FiniteElement &finite_element, - const Quadrature &quadrature, - const Mapping &mapping, - const UpdateFlags update_flags, - const unsigned int n_compositional_fields, - const unsigned int stokes_dofs_per_cell) + // Standard constructor + PcConstraintsAssembleData (const FiniteElement &finite_element, + const Quadrature &quadrature, + const Mapping &mapping, + const UpdateFlags update_flags, + const unsigned int n_compositional_fields, + const unsigned int stokes_dofs_per_cell) : finite_element_values (mapping, finite_element, quadrature, update_flags), @@ -1354,7 +1374,8 @@ namespace aspect {} - PcAssembleData (const PcAssembleData &scratch) + // Copy constructor + PcConstraintsAssembleData (const PcConstraintsAssembleData &scratch) : finite_element_values (scratch.finite_element_values.get_mapping(), scratch.finite_element_values.get_fe(), @@ -1364,11 +1385,9 @@ namespace aspect dof_component_indices( scratch.dof_component_indices), material_model_inputs(scratch.material_model_inputs), material_model_outputs(scratch.material_model_outputs) - { - - } + {} - virtual ~PcAssembleData () + virtual ~PcConstraintsAssembleData () {} FEValues finite_element_values; @@ -1382,9 +1401,9 @@ namespace aspect }; template - struct PcCopyData + struct PcConstraintsCopyData { - explicit PcCopyData (const unsigned int stokes_dofs_per_cell) + explicit PcConstraintsCopyData (const unsigned int stokes_dofs_per_cell) { nonzero_dof_indices.reserve(stokes_dofs_per_cell); } @@ -1392,74 +1411,64 @@ namespace aspect std::vector nonzero_dof_indices; }; + + // Assembler for setting the compaction pressure constraints in cells without melt. + // Note that the assembler does not actually assemble anything, it just uses the + // structure of an assembler to loop over all compaction pressure dofs and to insert + // the ones that are in melt cells into an IndexSet (nonzero_dof_indices) that is + // later needed for setting the constraints. template - class PcNonZeroAssembler : public SimulatorAccess + class PcNonZeroDofsAssembler : public SimulatorAccess { public: - PcNonZeroAssembler(IndexSet &nonzero_entries, - std::vector &is_melt_cell) - : nonzero_entries (nonzero_entries), - is_melt_cell (is_melt_cell) + PcNonZeroDofsAssembler(const std::vector &is_melt_cell, + IndexSet &nonzero_entries) + : is_melt_cell (is_melt_cell), + nonzero_entries (nonzero_entries) {} void - local_assemble (const typename DoFHandler::active_cell_iterator &cell, - PcAssembleData &scratch, - PcCopyData &data) + local_save_nonzero_pc_dofs (const typename DoFHandler::active_cell_iterator &cell, + PcConstraintsAssembleData &scratch, + PcConstraintsCopyData &data) { - data.nonzero_dof_indices.clear(); cell->get_dof_indices (scratch.local_dof_indices); - scratch.finite_element_values.reinit (cell); - - MeltHandler::create_material_model_outputs(scratch.material_model_outputs); - const unsigned int pc_component_index = this->introspection().variable("compaction pressure").first_component_index; const unsigned int pc_base_index = this->introspection().variable("compaction pressure").base_index; - Assert(dynamic_cast*>(&this->get_material_model()) != - NULL, ExcMessage("Error: The current material model needs to be derived from MeltInterface to use melt transport.")); - - // Prepare the data structures for assembly - scratch.finite_element_values.reinit (cell); - scratch.material_model_inputs.reinit(scratch.finite_element_values, - cell, - this->introspection(), - this->get_current_linearization_point(), - true); - - this->get_material_model().evaluate(scratch.material_model_inputs, - scratch.material_model_outputs); - - const double p_c_scale = dynamic_cast*>(&this->get_material_model())->p_c_scale(scratch.material_model_inputs,scratch.material_model_outputs,this->get_melt_handler(), true); - const bool is_melt_cell = (p_c_scale > 0.0); + Assert(cell.state() == IteratorState::valid, + ExcMessage("You are trying to access the property of a cell (if it is " + "a melt cell or not), but the cell iterator provided does not " + "point to a cell.")); + const bool is_melt_cell = this->get_melt_handler().is_melt_cell(cell); if (is_melt_cell) { + data.nonzero_dof_indices.resize(this->get_fe().base_element(pc_base_index).dofs_per_cell); for (unsigned int j=0; jget_fe().base_element(pc_base_index).dofs_per_cell; ++j) { const unsigned int pressure_idx = this->get_fe().component_to_system_index(pc_component_index, /*dof index within component=*/ j); - data.nonzero_dof_indices.push_back(scratch.local_dof_indices[pressure_idx]); + data.nonzero_dof_indices[j] = scratch.local_dof_indices[pressure_idx]; } } } void - copy_local_to_global (const PcCopyData &data) + copy_local_to_global (const PcConstraintsCopyData &data) { nonzero_entries.add_indices(data.nonzero_dof_indices.begin(), data.nonzero_dof_indices.end()); } private: - - IndexSet &nonzero_entries; const std::vector &is_melt_cell; + IndexSet &nonzero_entries; }; @@ -1470,7 +1479,7 @@ namespace aspect MeltHandler:: add_current_constraints(ConstraintMatrix &constraints) { - // First, fill the constrains matrix with the current constraints from the beginning + // First, fill the constraints matrix with the current constraints from the beginning // of the time step (that do not include the "melt cell" contributions). constraints.clear (); constraints.reinit (this->introspection().index_sets.system_relevant_set); @@ -1485,7 +1494,7 @@ namespace aspect is_melt_cell_vector.resize(this->get_dof_handler().get_triangulation().n_active_cells()); { - // find our "melt cells" by looking at K_D + // find the "melt cells" by looking at p_c_scale const unsigned int n_compositional_fields = this->introspection().n_compositional_fields; FEValues finite_element_values(this->get_mapping(), fe, quadrature_formula, cell_update_flags); @@ -1504,54 +1513,59 @@ namespace aspect { finite_element_values.reinit (cell); - this->compute_material_model_input_values (this->get_current_linearization_point(), - finite_element_values, - cell, - true, - material_model_inputs); + material_model_inputs.reinit(finite_element_values, + cell, + this->introspection(), + this->get_current_linearization_point(), + false); this->get_material_model().evaluate(material_model_inputs, material_model_outputs); - const double p_c_scale = dynamic_cast*>(&this->get_material_model())->p_c_scale(material_model_inputs,material_model_outputs,this->get_melt_handler(), false /*=consider_is_melt_cell*/); + const double p_c_scale = dynamic_cast*>(&this->get_material_model())->p_c_scale(material_model_inputs, + material_model_outputs, + this->get_melt_handler(), + false /*=consider_is_melt_cell*/); const bool is_melt_cell = (p_c_scale > 0.0); is_melt_cell_vector[cell->active_cell_index()] = is_melt_cell; } } - PcNonZeroAssembler assembler(nonzero_pc_dofs, is_melt_cell_vector); + PcNonZeroDofsAssembler assembler(is_melt_cell_vector, nonzero_pc_dofs); assembler.initialize_simulator(this->get_simulator()); typedef FilteredIterator::active_cell_iterator> CellFilter; - unsigned int stokes_dofs_per_cell = dim * fe.base_element(this->introspection().base_elements.velocities).dofs_per_cell - + fe.base_element(this->introspection().base_elements.pressure).dofs_per_cell - + fe.base_element(this->introspection().base_elements.pressure).dofs_per_cell; + const unsigned int stokes_dofs_per_cell = dim * fe.base_element(this->introspection().base_elements.velocities).dofs_per_cell + + fe.base_element(this->introspection().base_elements.pressure).dofs_per_cell + + fe.base_element(this->introspection().base_elements.pressure).dofs_per_cell; + // Here we call the assembler in order to save all compaction pressure dofs + // that are in melt cells (and are nonzero) in the nonzero_pc_dofs index set. WorkStream:: run (CellFilter (IteratorFilters::LocallyOwnedCell(), this->get_dof_handler().begin_active()), CellFilter (IteratorFilters::LocallyOwnedCell(), this->get_dof_handler().end()), - std_cxx11::bind (&PcNonZeroAssembler::local_assemble, + std_cxx11::bind (&PcNonZeroDofsAssembler::local_save_nonzero_pc_dofs, &assembler, std_cxx11::_1, std_cxx11::_2, std_cxx11::_3), - std_cxx11::bind (&PcNonZeroAssembler::copy_local_to_global, + std_cxx11::bind (&PcNonZeroDofsAssembler::copy_local_to_global, &assembler, std_cxx11::_1), - PcAssembleData (fe, quadrature_formula, - this->get_mapping(), - cell_update_flags, - this->introspection().n_compositional_fields, - stokes_dofs_per_cell), - PcCopyData (stokes_dofs_per_cell)); - - // first pick all relevant p_c's: + PcConstraintsAssembleData (fe, quadrature_formula, + this->get_mapping(), + cell_update_flags, + this->introspection().n_compositional_fields, + stokes_dofs_per_cell), + PcConstraintsCopyData (stokes_dofs_per_cell)); + + // For the constraints, first pick all relevant p_c dofs: IndexSet for_constraints = this->introspection().index_sets.system_relevant_set & Utilities::extract_locally_active_dofs_with_component(this->get_dof_handler(), this->introspection().variable("compaction pressure").component_mask); @@ -1559,9 +1573,7 @@ namespace aspect // now subtract the ones that are nonzero as computed above: for_constraints.subtract_set(nonzero_pc_dofs); - // and constrain those. Note that we do not constrain p_c dofs that - // are on the interface between a cell "with melt" and a cell "without melt" - // (using the definition "at least one quadrature point with p_c_scale>0"). + // and constrain those. constraints.add_lines(for_constraints); constraints.close(); } @@ -1595,7 +1607,16 @@ namespace aspect const bool is_melt_cell) const { const double ref_K_D = dynamic_cast*>(&this->get_material_model())->reference_darcy_coefficient(); - return is_melt_cell ? std::max(K_D, 1.e-3*ref_K_D) : 0; + return is_melt_cell ? std::max(K_D, melt_parameters.K_D_variation_threshold*ref_K_D) : 0; + } + + + template + const BoundaryFluidPressure::Interface & + MeltHandler:: + get_boundary_fluid_pressure () const + { + return *this->boundary_fluid_pressure.get(); } @@ -1612,10 +1633,6 @@ namespace aspect variables.insert(variables.begin()+1, VariableDeclaration( "fluid pressure", - (parameters.use_locally_conservative_discretization) - ? - std_cxx11::shared_ptr >(new FE_DGP(parameters.stokes_velocity_degree-1)) - : std_cxx11::shared_ptr >(new FE_Q(parameters.stokes_velocity_degree-1)), 1, 0)); // same block as p_c even without a direct solver! @@ -1623,12 +1640,11 @@ namespace aspect variables.insert(variables.begin()+2, VariableDeclaration( "compaction pressure", - (parameters.use_locally_conservative_discretization || use_discontinuous_p_c) + melt_parameters.use_discontinuous_p_c ? std_cxx11::shared_ptr >(new FE_DGP(parameters.stokes_velocity_degree-1)) : - std_cxx11::shared_ptr >(new FE_Q(parameters.stokes_velocity_degree-1)) - , + std_cxx11::shared_ptr >(new FE_Q(parameters.stokes_velocity_degree-1)), 1, 1)); @@ -1689,92 +1705,114 @@ namespace aspect } - template - void - MeltHandler:: - declare_parameters (ParameterHandler &prm) + namespace Melt { - prm.enter_subsection ("Melt settings"); + template + void + Parameters:: + declare_parameters (ParameterHandler &prm) { - prm.declare_entry ("Melt transport threshold", "1e-3", - Patterns::Double (), - "The porosity limit for melt migration. For smaller porosities, the equations " - "reduce to the Stokes equations and do neglect melt transport. Only used " - "if Include melt transport is true. "); - prm.declare_entry ("Heat advection by melt", "false", - Patterns::Bool (), - "Whether to use a porosity weighted average of the melt and solid velocity " - "to advect heat in the temperature equation or not. If this is set to true, " - "additional terms are assembled on the left-hand side of the temperature " - "advection equation. Only used if Include melt transport is true. " - "If this is set to false, only the solid velocity is used (as in models " - "without melt migration)."); - prm.declare_entry ("Use discontinuous compaction pressure", "true", - Patterns::Bool (), - "Whether to use a discontinuous element for the compaction pressure or not. " - "From our preliminary tests, continuous elements seem to work better in models " - "where the porosity is > 0 everywhere in the domain, and discontinuous elements " - "work better in models where in parts of the domain the porosity = 0."); - prm.declare_entry ("Average melt velocity", "true", - Patterns::Bool (), - "Whether to cell-wise average the material properties that are used to " - "compute the melt velocity or not. The melt velocity is computed as the " - "sum of the solid velocity and the phase separation flux " - "$ - K_D / \\phi (\\nabla p_f - \\rho_f \\mathbf g)$. " - "If this parameter is set to true, $K_D$ and $\\phi$ will be averaged " - "cell-wise in the computation of the phase separation flux. " - "This is useful because in some models the melt velocity can have spikes " - "close to the interface between regions of melt and no melt, as both $K_D$ " - "and $\\phi$ go to zero for vanishing melt fraction. As the melt velocity is " - "used for computing the time step size, and in models that use heat " - "transport by melt or shear heating of melt, setting this parameter to true " - "can speed up the model and make it mode stable. In computations where " - "accuracy and convergence behavior of the melt velocity is important " - "(like in benchmark cases with an analytical solution), this parameter " - "should probably be set to 'false'."); + prm.enter_subsection ("Melt settings"); + { + prm.declare_entry ("Darcy coefficient variation threshold", "1e-3", + Patterns::Double (), + "The factor by how much the Darcy coefficient K_D in a cell can be smaller than " + "the reference Darcy coefficient for this cell still to be considered a melt cell " + "(for which the melt transport equations are solved). For smaller Darcy coefficients, " + "the Stokes equations (without melt) are solved instead. Only used if " + "``Include melt transport'' is true. "); + prm.declare_entry ("Heat advection by melt", "false", + Patterns::Bool (), + "Whether to use a porosity weighted average of the melt and solid velocity " + "to advect heat in the temperature equation or not. If this is set to true, " + "additional terms are assembled on the left-hand side of the temperature " + "advection equation. Only used if Include melt transport is true. " + "If this is set to false, only the solid velocity is used (as in models " + "without melt migration)."); + prm.declare_entry ("Use discontinuous compaction pressure", "true", + Patterns::Bool (), + "Whether to use a discontinuous element for the compaction pressure or not. " + "From our preliminary tests, continuous elements seem to work better in models " + "where the porosity is > 0 everywhere in the domain, and discontinuous elements " + "work better in models where in parts of the domain the porosity = 0."); + prm.declare_entry ("Average melt velocity", "true", + Patterns::Bool (), + "Whether to cell-wise average the material properties that are used to " + "compute the melt velocity or not. The melt velocity is computed as the " + "sum of the solid velocity and the phase separation flux " + "$ - K_D / \\phi (\\nabla p_f - \\rho_f \\mathbf g)$. " + "If this parameter is set to true, $K_D$ and $\\phi$ will be averaged " + "cell-wise in the computation of the phase separation flux. " + "This is useful because in some models the melt velocity can have spikes " + "close to the interface between regions of melt and no melt, as both $K_D$ " + "and $\\phi$ go to zero for vanishing melt fraction. As the melt velocity is " + "used for computing the time step size, and in models that use heat " + "transport by melt or shear heating of melt, setting this parameter to true " + "can speed up the model and make it mode stable. In computations where " + "accuracy and convergence behavior of the melt velocity is important " + "(like in benchmark cases with an analytical solution), this parameter " + "should probably be set to 'false'."); + } + prm.leave_subsection(); + + BoundaryFluidPressure::declare_parameters (prm); } - prm.leave_subsection(); - BoundaryFluidPressure::declare_parameters (prm); + template + void + Parameters:: + parse_parameters (ParameterHandler &prm) + { + prm.enter_subsection ("Melt settings"); + { + K_D_variation_threshold = prm.get_double("Darcy coefficient variation threshold"); + heat_advection_by_melt = prm.get_bool("Heat advection by melt"); + use_discontinuous_p_c = prm.get_bool("Use discontinuous compaction pressure"); + average_melt_velocity = prm.get_bool("Average melt velocity"); + } + prm.leave_subsection(); + } } template MeltHandler::MeltHandler (ParameterHandler &prm) + : + boundary_fluid_pressure(BoundaryFluidPressure::create_boundary_fluid_pressure(prm)) { - parse_parameters(prm); + melt_parameters.parse_parameters(prm); + boundary_fluid_pressure->parse_parameters(prm); } + template void - MeltHandler::initialize_simulator (const Simulator &simulator_object) + MeltHandler::initialize () const { - this->SimulatorAccess::initialize_simulator(simulator_object); - if (SimulatorAccess *sim = dynamic_cast*>(boundary_fluid_pressure.get())) - sim->initialize_simulator (simulator_object); - boundary_fluid_pressure->initialize (); - - if (use_discontinuous_p_c) + // The additional terms in the temperature systems have not been ported + // to the DG formulation: + AssertThrow(!this->get_parameters().use_discontinuous_temperature_discretization && + !this->get_parameters().use_discontinuous_composition_discretization, + ExcMessage("Using discontinuous elements for temperature " + "or composition in models with melt transport is currently not implemented.") ); + if (melt_parameters.use_discontinuous_p_c) AssertThrow(!this->model_has_prescribed_stokes_solution(), ExcMessage("You can not use a discontinuous p_c in a model " "with a presribed Stokes solution.")); + // We can not have a DG p_f. + AssertThrow(!this->get_parameters().use_locally_conservative_discretization, + ExcMessage ("Discontinuous elements for the fluid pressure " + "are not supported in models with melt transport.")); } + template void - MeltHandler:: - parse_parameters (ParameterHandler &prm) + MeltHandler::initialize_simulator (const Simulator &simulator_object) { - prm.enter_subsection ("Melt settings"); - { - melt_transport_threshold = prm.get_double("Melt transport threshold"); - heat_advection_by_melt = prm.get_bool("Heat advection by melt"); - use_discontinuous_p_c = prm.get_bool("Use discontinuous compaction pressure"); - average_melt_velocity = prm.get_bool("Average melt velocity"); - } - prm.leave_subsection(); - - boundary_fluid_pressure.reset(BoundaryFluidPressure::create_boundary_fluid_pressure(prm)); - boundary_fluid_pressure->parse_parameters (prm); + this->SimulatorAccess::initialize_simulator(simulator_object); + if (SimulatorAccess *sim = dynamic_cast*>(boundary_fluid_pressure.get())) + sim->initialize_simulator (simulator_object); + boundary_fluid_pressure->initialize (); } @@ -1889,6 +1927,10 @@ namespace aspect class \ MeltHandler; \ \ + namespace Melt \ + { \ + template struct Parameters; \ + } \ namespace Assemblers \ { \ template class MeltInterface; \ diff --git a/source/simulator/parameters.cc b/source/simulator/parameters.cc index 067ff2ae753..7a449aa01ff 100644 --- a/source/simulator/parameters.cc +++ b/source/simulator/parameters.cc @@ -1278,20 +1278,6 @@ namespace aspect "is of one degree lower and continuous, and if you selected " "a linear element for the velocity, you'd need a continuous " "element of degree zero for the pressure, which does not exist.")) - - if (include_melt_transport) - { - // The additional terms in the temperature systems have not been ported - // to the DG formulation: - AssertThrow(!use_discontinuous_temperature_discretization - && !use_discontinuous_composition_discretization, - ExcMessage ("Using discontinuous elements for temperature " - "or composition in models with melt transport is currently not implemented.")); - // We can not have a DG p_f. - AssertThrow(!use_locally_conservative_discretization, - ExcMessage ("Discontinuous elements for the pressure " - "in models with melt transport are not supported")); - } } prm.leave_subsection (); @@ -1684,7 +1670,7 @@ namespace aspect void Simulator::declare_parameters (ParameterHandler &prm) { Parameters::declare_parameters (prm); - MeltHandler::declare_parameters (prm); + Melt::Parameters::declare_parameters (prm); Newton::Parameters::declare_parameters (prm); Postprocess::Manager::declare_parameters (prm); MeshRefinement::Manager::declare_parameters (prm); diff --git a/source/simulator/solver.cc b/source/simulator/solver.cc index 22077868542..3c93297ca80 100644 --- a/source/simulator/solver.cc +++ b/source/simulator/solver.cc @@ -280,10 +280,10 @@ namespace aspect // first solve with the bottom left block, which we have built // as a mass matrix with the inverse of the viscosity { - aspect::SolverControl solver_control(1000, src.block(1).l2_norm() * S_block_tolerance); + SolverControl solver_control(1000, src.block(1).l2_norm() * S_block_tolerance); #ifdef ASPECT_USE_PETSC - SolverGMRES solver(solver_control); + SolverCG solver(solver_control); #else TrilinosWrappers::SolverCG solver(solver_control); #endif @@ -836,7 +836,6 @@ namespace aspect { // output solver history std::ofstream f((parameters.output_directory+"solver_history.txt").c_str()); - f << std::setprecision(16); // Only request the solver history if a history has actually been created if (parameters.n_cheap_stokes_solver_steps > 0) diff --git a/tests/compaction_length_refinement.prm b/tests/compaction_length_refinement.prm index 8144621e685..4cba33234ee 100644 --- a/tests/compaction_length_refinement.prm +++ b/tests/compaction_length_refinement.prm @@ -110,10 +110,6 @@ subsection Model settings set Include melt transport = true end -subsection Melt settings - set Melt transport threshold = 0 -end - subsection Postprocess diff --git a/tests/composition_active_with_melt.prm b/tests/composition_active_with_melt.prm index 591c8240bfc..929c52ee0d9 100644 --- a/tests/composition_active_with_melt.prm +++ b/tests/composition_active_with_melt.prm @@ -36,7 +36,6 @@ subsection Model settings end subsection Melt settings - set Melt transport threshold = 0.0 set Use discontinuous compaction pressure = false end diff --git a/tests/free_surface_blob_nonzero_melt.prm b/tests/free_surface_blob_nonzero_melt.prm index 4d482f8d7d0..54b13b1964d 100644 --- a/tests/free_surface_blob_nonzero_melt.prm +++ b/tests/free_surface_blob_nonzero_melt.prm @@ -113,10 +113,6 @@ subsection Model settings set Include melt transport = true end -subsection Melt settings - set Melt transport threshold = 0 -end - subsection Free surface set Free surface stabilization theta = 0.5 end diff --git a/tests/global_melt.prm b/tests/global_melt.prm index 95f2288b792..e6b49a40966 100644 --- a/tests/global_melt.prm +++ b/tests/global_melt.prm @@ -183,10 +183,6 @@ subsection Model settings set Include melt transport = true end -subsection Melt settings - set Melt transport threshold = 1e-4 -end - subsection Postprocess diff --git a/tests/global_melt_parallel.prm b/tests/global_melt_parallel.prm index 014ebf5a44e..06517361435 100644 --- a/tests/global_melt_parallel.prm +++ b/tests/global_melt_parallel.prm @@ -185,10 +185,6 @@ subsection Model settings set Include melt transport = true end -subsection Melt settings - set Melt transport threshold = 1e-4 -end - subsection Postprocess diff --git a/tests/initial_porosity.prm b/tests/initial_porosity.prm index 366aa56b7a9..185e95ea9d7 100644 --- a/tests/initial_porosity.prm +++ b/tests/initial_porosity.prm @@ -184,10 +184,6 @@ subsection Model settings set Include melt transport = true end -subsection Melt settings - set Melt transport threshold = 1e-4 -end - subsection Postprocess diff --git a/tests/melt_and_traction.prm b/tests/melt_and_traction.prm index 2cf465c77e0..e1cb35c0a55 100644 --- a/tests/melt_and_traction.prm +++ b/tests/melt_and_traction.prm @@ -60,10 +60,6 @@ end ##################### Settings for melt transport ######################## -subsection Melt settings - set Melt transport threshold = 0.0 -end - subsection Compositional fields set Number of fields = 2 set Names of fields = porosity, peridotite diff --git a/tests/melt_force_vector.prm b/tests/melt_force_vector.prm index 864021b103a..a5fcba31ec6 100644 --- a/tests/melt_force_vector.prm +++ b/tests/melt_force_vector.prm @@ -137,7 +137,6 @@ subsection Model settings set Zero velocity boundary indicators = # default: set Include melt transport = true -# set Melt transport threshold = 0 set Enable additional Stokes RHS = true end diff --git a/tests/melt_material_4.prm b/tests/melt_material_4.prm index e895a920b7b..4f7088969db 100644 --- a/tests/melt_material_4.prm +++ b/tests/melt_material_4.prm @@ -147,10 +147,6 @@ subsection Model settings set Include melt transport = true end -subsection Melt settings - set Melt transport threshold = 0 -end - subsection Postprocess diff --git a/tests/melt_track.prm b/tests/melt_track.prm index ba9a9b26aaf..74da2630b69 100644 --- a/tests/melt_track.prm +++ b/tests/melt_track.prm @@ -147,10 +147,6 @@ subsection Model settings set Include melt transport = true end -subsection Melt settings - set Melt transport threshold = 5e-3 -end - subsection Postprocess diff --git a/tests/melt_transport.prm b/tests/melt_transport.prm index 5a76f8e249a..ed94ec3080d 100644 --- a/tests/melt_transport.prm +++ b/tests/melt_transport.prm @@ -139,10 +139,6 @@ subsection Model settings set Include melt transport = true end -subsection Melt settings - set Melt transport threshold = 1e-4 -end - subsection Postprocess diff --git a/tests/melt_velocity_boundary_conditions.prm b/tests/melt_velocity_boundary_conditions.prm index 591ea72b54a..ccc267241fa 100644 --- a/tests/melt_velocity_boundary_conditions.prm +++ b/tests/melt_velocity_boundary_conditions.prm @@ -115,10 +115,6 @@ subsection Model settings set Include melt transport = true end -subsection Melt settings - set Melt transport threshold = 0 -end - subsection Free surface set Free surface stabilization theta = 0.5 end diff --git a/tests/n_expensive_stokes_solver_steps/solver_history.txt b/tests/n_expensive_stokes_solver_steps/solver_history.txt index bc67b2f2b7b..1d719b577bf 100644 --- a/tests/n_expensive_stokes_solver_steps/solver_history.txt +++ b/tests/n_expensive_stokes_solver_steps/solver_history.txt @@ -1,6 +1,6 @@ -0 686983212.7251304 -1 686318864.7549431 -2 16323811.30518775 -3 961925.3081440757 -4 196210.4232824568 -5 6802.435628751216 +0 6.86983e+08 +1 6.86319e+08 +2 1.63238e+07 +3 961925 +4 196210 +5 6802.44 diff --git a/tests/particle_melt_advection.prm b/tests/particle_melt_advection.prm index 818c2fd02e3..6a136ff0a31 100644 --- a/tests/particle_melt_advection.prm +++ b/tests/particle_melt_advection.prm @@ -36,7 +36,7 @@ subsection Compositional fields end subsection Boundary temperature model - set Model name = box + set List of model names = box subsection Box set Bottom temperature = 1600 @@ -119,10 +119,6 @@ subsection Model settings set Include melt transport = true end -subsection Melt settings - set Melt transport threshold = 1e-4 -end - subsection Postprocess diff --git a/tests/particle_property_composition.prm b/tests/particle_property_composition.prm index 1f668ca5607..33c38cc665c 100644 --- a/tests/particle_property_composition.prm +++ b/tests/particle_property_composition.prm @@ -147,10 +147,6 @@ subsection Model settings set Include melt transport = true end -subsection Melt settings - set Melt transport threshold = 5e-3 -end - subsection Postprocess diff --git a/tests/rising_melt_blob.prm b/tests/rising_melt_blob.prm index be771683e05..83e1ebe475b 100644 --- a/tests/rising_melt_blob.prm +++ b/tests/rising_melt_blob.prm @@ -1,17 +1,14 @@ ######################################################### -# This is a test for the melting rate functionality -# It uses a new material model similar to the melt -# simple model, but with a different melting rate. +# This is a test for the new melt solver. It checks if +# melt transport works correctly in a model that has both +# regions with and without melt (so only some of the cells +# are "melt cells"). -# The melting rate has the form of a Gaussian function -# and every timestep a fixed amount of material is melting. +# The initial porosity is zero everywhere, except for a +# Gaussian in the center of the model (the "melt blob"). +# As there is no melting or freezing, the blob should +# just rise, conserving the total porosity. -# The amplitude of the Gaussian is 0.0001, so both the -# porosity and the peridotite field should increase by -# this amount (visible in the maximum value). - -# Note that the maximum porosity can deviate slightly, -# as additional melt is created by compaction. ############### Global parameters set Dimension = 2 @@ -27,10 +24,6 @@ set CFL number = 0.2 set Pressure normalization = surface set Surface pressure = 0 -subsection Melt settings - set Melt transport threshold = 0.0 -end - ############### Parameters describing the model # Let us here choose again a box domain @@ -71,7 +64,7 @@ end subsection Boundary temperature model - set Model name = box + set List of model names = box subsection Box set Bottom temperature = 1600 diff --git a/tests/rising_melt_blob_freezing.prm b/tests/rising_melt_blob_freezing.prm index 3df47876375..25aab0e057d 100644 --- a/tests/rising_melt_blob_freezing.prm +++ b/tests/rising_melt_blob_freezing.prm @@ -26,10 +26,6 @@ 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 diff --git a/tests/shear_bands.prm b/tests/shear_bands.prm index 5829ec275e8..be9ba114bfb 100644 --- a/tests/shear_bands.prm +++ b/tests/shear_bands.prm @@ -131,10 +131,6 @@ subsection Model settings set Include melt transport = true end -subsection Melt settings - set Melt transport threshold = 1e-4 -end - subsection Postprocess set List of postprocessors = visualization,composition statistics,velocity statistics diff --git a/tests/shear_heating_with_melt.prm b/tests/shear_heating_with_melt.prm index 7d42b3798eb..feb92014490 100644 --- a/tests/shear_heating_with_melt.prm +++ b/tests/shear_heating_with_melt.prm @@ -4,8 +4,10 @@ # from the bottom boundary and melting adiabatically. This # leads to dissipation heating due to the relative movement of # melt and solid. As these terms scale with 1/porosity, this -# can lead to instabilities. This should be fixed by using a -# Melt transport threshold > 0, as done in this test. +# can lead to instabilities. However, this can be avoided by +# choosing a resolution that is high enough and choosing the +# Darcy coefficient variation threshold that is not too small +# (here, the default value is sufficient). # The test is based on the model of 1D upwelling melting in a # plume with melt-related heating in the Supplementary material @@ -195,7 +197,6 @@ subsection Model settings end subsection Melt settings - set Melt transport threshold = 1e-3 set Heat advection by melt = true end diff --git a/tests/solitary_wave.prm b/tests/solitary_wave.prm index 7dd6acee493..b6a1928c357 100644 --- a/tests/solitary_wave.prm +++ b/tests/solitary_wave.prm @@ -145,10 +145,6 @@ subsection Model settings set Include melt transport = true end -subsection Melt settings - set Melt transport threshold = 1e-4 -end - subsection Postprocess diff --git a/tests/visco_plastic_derivatives.prm b/tests/visco_plastic_derivatives.prm index 3c34e4556ed..6153c844692 100644 --- a/tests/visco_plastic_derivatives.prm +++ b/tests/visco_plastic_derivatives.prm @@ -77,7 +77,7 @@ end # Boundary composition specification subsection Boundary composition model - set Model name = initial composition + set List of model names = initial composition end # Material model (values for background material)