diff --git a/benchmarks/advection/rotate_shape_dg_bp.prm b/benchmarks/advection/rotate_shape_dg_bp.prm new file mode 100644 index 00000000000..b1c2b7309a1 --- /dev/null +++ b/benchmarks/advection/rotate_shape_dg_bp.prm @@ -0,0 +1,22 @@ +# Like rotate_shape_supg but with discontinuous Galerkin method +# and bound-preserving limiter + +set Dimension = 2 + +include rotate_shape_supg.prm + +set Output directory = output-rotate-shape-dg-bp + +subsection Discretization + set Use discontinuous temperature discretization = true + set Use discontinuous composition discretization = true + + subsection Stabilization parameters + set Limiter for discontinuous temperature solution = bound preserving + set Limiter for discontinuous composition solution = bound preserving + set Global temperature maximum = 1 + set Global temperature minimum = 0 + set Global composition maximum = 1 + set Global composition minimum = 0 + end +end diff --git a/benchmarks/advection/rotate_shape_dg_weno.prm b/benchmarks/advection/rotate_shape_dg_weno.prm new file mode 100644 index 00000000000..aae23e8b1e7 --- /dev/null +++ b/benchmarks/advection/rotate_shape_dg_weno.prm @@ -0,0 +1,28 @@ +# Like rotate_shape_supg but with discontinuous Galerkin method +# and WENO limiter + +set Dimension = 2 + +include rotate_shape_supg.prm + +set Output directory = output-rotate-shape-dg-weno + +subsection Discretization + set Use discontinuous temperature discretization = true + set Use discontinuous composition discretization = true + + subsection Stabilization parameters + set Limiter for discontinuous temperature solution = WENO + set Limiter for discontinuous composition solution = WENO + end +end + +subsection Postprocess + subsection Visualization + set List of output variables = kxrcf indicator + + subsection KXRCF indicator + set Name of advection field = C_1 + end + end +end diff --git a/benchmarks/buiter_et_al_2008_jgr/figure_6_1e20.prm b/benchmarks/buiter_et_al_2008_jgr/figure_6_1e20.prm index d3006d93625..e72349df709 100644 --- a/benchmarks/buiter_et_al_2008_jgr/figure_6_1e20.prm +++ b/benchmarks/buiter_et_al_2008_jgr/figure_6_1e20.prm @@ -100,7 +100,6 @@ subsection Discretization set Use discontinuous composition discretization = true subsection Stabilization parameters - set Use limiter for discontinuous composition solution = true # apply the limiter to the DG solutions set Global composition maximum = 100.0, 1.0, 1.0, 1.0 set Global composition minimum = 0.0, 0.0, 0.0, 1.0 end diff --git a/benchmarks/buiter_et_al_2016_jsg/doc/convergence.part.prm b/benchmarks/buiter_et_al_2016_jsg/doc/convergence.part.prm index 544ad88b0ac..6e2907fd476 100644 --- a/benchmarks/buiter_et_al_2016_jsg/doc/convergence.part.prm +++ b/benchmarks/buiter_et_al_2016_jsg/doc/convergence.part.prm @@ -21,9 +21,9 @@ subsection Discretization set Use discontinuous composition discretization = true subsection Stabilization parameters - set Use limiter for discontinuous composition solution = true set Global composition maximum = 1, 1, 1, 1, 1, 100 set Global composition minimum = 0, 0, 0, 0, 0, 0 + set Limiter for discontinuous composition solution = bound preserving end end diff --git a/benchmarks/buiter_et_al_2016_jsg/exp_1.prm b/benchmarks/buiter_et_al_2016_jsg/exp_1.prm index 51b2116626e..976b4f505b3 100644 --- a/benchmarks/buiter_et_al_2016_jsg/exp_1.prm +++ b/benchmarks/buiter_et_al_2016_jsg/exp_1.prm @@ -106,11 +106,10 @@ subsection Discretization set Use discontinuous composition discretization = true subsection Stabilization parameters - set Use limiter for discontinuous composition solution = true - # [sand, bound, block, sticky_air, total_strain] set Global composition maximum = 1., 1., 1., 1., 1e5 set Global composition minimum = 0., 0., 0., 0., 0. + set Limiter for discontinuous composition solution = bound preserving end end diff --git a/benchmarks/buiter_et_al_2016_jsg/exp_2_high_resolution.prm b/benchmarks/buiter_et_al_2016_jsg/exp_2_high_resolution.prm index b2496ed311c..dba24612726 100644 --- a/benchmarks/buiter_et_al_2016_jsg/exp_2_high_resolution.prm +++ b/benchmarks/buiter_et_al_2016_jsg/exp_2_high_resolution.prm @@ -48,9 +48,9 @@ subsection Discretization set Use discontinuous composition discretization = true subsection Stabilization parameters - set Use limiter for discontinuous composition solution = true set Global composition maximum = 1, 1, 1, 1, 1, 100 set Global composition minimum = 0, 0, 0, 0, 0, 0 + set Limiter for discontinuous composition solution = bound preserving end end diff --git a/benchmarks/infill_density/infill_ascii.prm b/benchmarks/infill_density/infill_ascii.prm index dbed3116c5c..aad1f865492 100644 --- a/benchmarks/infill_density/infill_ascii.prm +++ b/benchmarks/infill_density/infill_ascii.prm @@ -72,9 +72,9 @@ subsection Discretization set Use discontinuous composition discretization = true subsection Stabilization parameters - set Use limiter for discontinuous composition solution = true set Global composition maximum = 1.e13, 1.e13, 1.e13, 1.0 set Global composition minimum = -1.e13, -1.e13, -1.e13, 0.0 + set Limiter for discontinuous composition solution = bound preserving end end diff --git a/benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/input_v.prm b/benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/input_v.prm index 37534cc3107..7b80ff14209 100644 --- a/benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/input_v.prm +++ b/benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/input_v.prm @@ -2,7 +2,7 @@ # flow benchmark. This is used to test the velocity boundary conditions # of the Newton Stokes solver. -set Output directory = a2_itIMPES_ST1e-20_UFSfalse_NSP-SPD_NSA-none_C0_g6_ag0_AEWfalse_SF9e-1_NLT1e-14_ABT1e-2_LT1e-5_mLT9e-1_I150_P3_EW1_theta1_LS5_RSMfalse_AV-1_n5 +set Output directory = a2_itsingle Advection, single Stokes_ST1e-20_UFSfalse_NSP-SPD_NSA-none_C0_g6_ag0_AEWfalse_SF9e-1_NLT1e-14_ABT1e-2_LT1e-5_mLT9e-1_I150_P3_EW1_theta1_LS5_RSMfalse_AV-1_n5 set Dimension = 2 set CFL number = 1.0 set Maximum time step = 1 diff --git a/benchmarks/newton_solver_benchmark_set/tosi_et_al_2015/input.prm b/benchmarks/newton_solver_benchmark_set/tosi_et_al_2015/input.prm index b6cca66b629..86379ea1226 100644 --- a/benchmarks/newton_solver_benchmark_set/tosi_et_al_2015/input.prm +++ b/benchmarks/newton_solver_benchmark_set/tosi_et_al_2015/input.prm @@ -146,9 +146,9 @@ end # - Vrms over the whole domain (velocity statistics) # - Vrms along the surface (velocity boundary statistics) # - Vmax along the surface (velocity boundary statistics) -# - average rate of viscous dissipation (tosi viscous dissipation statistics) -# - average rate of work against gravity (tosi viscous dissipation statistics) -# - percentage error difference work and dissipation (tosi viscous dissipation statistics) +# - average rate of viscous dissipation (tosi heating statistics) +# - average rate of work against gravity (tosi heating statistics) +# - percentage error difference work and dissipation (tosi heating statistics) # # Tosi et al. (2015) also report depth averages; averages for viscosity # and temperature are produced every 0.05 units of time. diff --git a/benchmarks/nsinker/nsinker.prm b/benchmarks/nsinker/nsinker.prm index 3e7adb91f88..0716eaac90a 100644 --- a/benchmarks/nsinker/nsinker.prm +++ b/benchmarks/nsinker/nsinker.prm @@ -13,7 +13,6 @@ set Nonlinear solver scheme = no Advection, single Stokes set Max nonlinear iterations = 1 set Use years in output instead of seconds = false - # Follow as closely as possible the parameters from Rudi et al. (2017) subsection Solver parameters subsection Stokes solver parameters @@ -25,8 +24,9 @@ subsection Solver parameters set Use weighted BFBT for Schur complement = false set Krylov method for cheap solver steps = GMRES end + subsection AMG parameters - set AMG aggregation threshold = 0.02 + set AMG aggregation threshold = 0.02 end end diff --git a/benchmarks/nsinker/nsinker_bfbt.prm b/benchmarks/nsinker/nsinker_bfbt.prm index ed104447980..7de5ffda08f 100644 --- a/benchmarks/nsinker/nsinker_bfbt.prm +++ b/benchmarks/nsinker/nsinker_bfbt.prm @@ -14,10 +14,9 @@ set Nonlinear solver scheme = no Advection, single Stokes set Max nonlinear iterations = 1 set Use years in output instead of seconds = false - # Follow as closely as possible the parameters from Rudi et al. (2017) subsection Solver parameters - subsection Stokes solver parameters + subsection Stokes solver parameters set Use full A block as preconditioner = true set Number of cheap Stokes solver steps = 500 set Maximum number of expensive Stokes solver steps = 1000 @@ -26,8 +25,9 @@ subsection Solver parameters set Use weighted BFBT for Schur complement = true set Krylov method for cheap solver steps = GMRES end + subsection AMG parameters - set AMG aggregation threshold = 0.02 + set AMG aggregation threshold = 0.02 end end diff --git a/benchmarks/particle_integration_scheme/circle.prm b/benchmarks/particle_integration_scheme/circle.prm index 83ac384ab47..3b7a17da328 100644 --- a/benchmarks/particle_integration_scheme/circle.prm +++ b/benchmarks/particle_integration_scheme/circle.prm @@ -5,7 +5,7 @@ set Use years in output instead of seconds = false set CFL number = 1.0 set Output directory = circle_euler_1.0 set Timing output frequency = 100 -set Nonlinear solver scheme = Advection only +set Nonlinear solver scheme = single Advection, no Stokes subsection Geometry model set Model name = box diff --git a/benchmarks/rayleigh_taylor_instability_free_surface/kaus_rayleigh_taylor_instability.prm b/benchmarks/rayleigh_taylor_instability_free_surface/kaus_rayleigh_taylor_instability.prm index 08f66c0330c..96cbf4deadc 100644 --- a/benchmarks/rayleigh_taylor_instability_free_surface/kaus_rayleigh_taylor_instability.prm +++ b/benchmarks/rayleigh_taylor_instability_free_surface/kaus_rayleigh_taylor_instability.prm @@ -27,9 +27,9 @@ subsection Discretization set Use discontinuous composition discretization = true subsection Stabilization parameters - set Use limiter for discontinuous composition solution = true set Global composition maximum = 1 set Global composition minimum = 0 + set Limiter for discontinuous composition solution = bound preserving end end diff --git a/benchmarks/solcx/compositional_fields/solcx_compositional_fields.prm b/benchmarks/solcx/compositional_fields/solcx_compositional_fields.prm index 724762bcdfb..19f189b6c0a 100644 --- a/benchmarks/solcx/compositional_fields/solcx_compositional_fields.prm +++ b/benchmarks/solcx/compositional_fields/solcx_compositional_fields.prm @@ -70,9 +70,9 @@ subsection Discretization set Use locally conservative discretization = true subsection Stabilization parameters - set Use limiter for discontinuous composition solution = true # apply the limiter to the DG solutions set Global composition maximum = 1, 1e6 set Global composition minimum = -1, 1 + set Limiter for discontinuous temperature solution = bound preserving # apply the limiter to the DG solutions end end diff --git a/benchmarks/solcx/compositional_fields/solcx_particles.prm b/benchmarks/solcx/compositional_fields/solcx_particles.prm index 624d2616aaa..e742152af5c 100644 --- a/benchmarks/solcx/compositional_fields/solcx_particles.prm +++ b/benchmarks/solcx/compositional_fields/solcx_particles.prm @@ -55,7 +55,7 @@ subsection Discretization set Use discontinuous composition discretization = true subsection Stabilization parameters - set Use limiter for discontinuous composition solution = false # apply the limiter to the DG solutions + set Limiter for discontinuous temperature solution = none # apply the limiter to the DG solutions end end diff --git a/benchmarks/solubility/plugin/solubility.cc b/benchmarks/solubility/plugin/solubility.cc index c899cb5ed8c..896ed59ffa9 100644 --- a/benchmarks/solubility/plugin/solubility.cc +++ b/benchmarks/solubility/plugin/solubility.cc @@ -359,7 +359,7 @@ namespace aspect melt_fractions (const MaterialModel::MaterialModelInputs &in, std::vector &melt_fractions) const { - for (unsigned int q=0; qintrospection().compositional_index_for_name("water_content"); const unsigned int porosity_idx = this->introspection().compositional_index_for_name("porosity"); @@ -399,7 +399,7 @@ namespace aspect if (this->get_parameters().use_operator_splitting && out.template get_additional_output>() == nullptr) { - const unsigned int n_points = out.viscosities.size(); + const unsigned int n_points = out.n_evaluation_points(); out.additional_outputs.push_back( std::make_unique> (n_points, this->n_compositional_fields())); } diff --git a/benchmarks/viscoelastic_beam_modified/20km_opentopbot_weno.prm b/benchmarks/viscoelastic_beam_modified/20km_opentopbot_weno.prm new file mode 100644 index 00000000000..5f10ad30d36 --- /dev/null +++ b/benchmarks/viscoelastic_beam_modified/20km_opentopbot_weno.prm @@ -0,0 +1,35 @@ +# This parameter file modifies the benchmark 20kmh_05drho.prm to use the WENO limiter +# for viscoelastic stress components. + +include $ASPECT_SOURCE_DIR/benchmarks/viscoelastic_bending_beam/20km_opentopbot.prm + +set Output directory = 20kmh_05drho_weno + +subsection Discretization + subsection Stabilization parameters + # Use WENO limiter for viscoelastic stress components, while using the + # bound-preserving limiter for the chemical field + set Limiter for discontinuous composition solution = WENO, WENO, WENO, bound preserving + + # In this benchmark, the spurious oscillations at the boundaries of the beam + # are extremely serious, so we set the linear weight of neighbor cells one + # magnitude greater than the default value to make the reconstructed solution + # smoother. + set WENO linear weight of neighbor cells = 0.01 + + # The global maximum and minimum values are only used by the chemical field. + set Global composition maximum = 1.0 + set Global composition minimum = 0.0 + end +end + +# Visualize the KXRCF indicator for ve_stress_xy +subsection Postprocess + subsection Visualization + set List of output variables = material properties, strain rate, kxrcf indicator + + subsection KXRCF indicator + set Name of advection field = ve_stress_xy + end + end +end diff --git a/benchmarks/viscoelastic_bending_beam/viscoelastic_bending_beam.prm b/benchmarks/viscoelastic_bending_beam/viscoelastic_bending_beam.prm index 3d14a096e72..ca656ef8389 100644 --- a/benchmarks/viscoelastic_bending_beam/viscoelastic_bending_beam.prm +++ b/benchmarks/viscoelastic_bending_beam/viscoelastic_bending_beam.prm @@ -74,9 +74,9 @@ subsection Discretization set Use discontinuous composition discretization = true subsection Stabilization parameters - set Use limiter for discontinuous composition solution = true set Global composition maximum = 1.e11, 1.e11, 1.e11, 1.0 set Global composition minimum = -1.e11, -1.e11, -1.e11, 0.0 + set Limiter for discontinuous composition solution = bound preserving end end diff --git a/benchmarks/viscoelastic_bending_beam/viscoelastic_bending_beam_weno.prm b/benchmarks/viscoelastic_bending_beam/viscoelastic_bending_beam_weno.prm new file mode 100644 index 00000000000..34f163cfdd1 --- /dev/null +++ b/benchmarks/viscoelastic_bending_beam/viscoelastic_bending_beam_weno.prm @@ -0,0 +1,35 @@ +# This parameter file modifies the benchmark viscoealastic_bending_beam.prm +# to use the WENO limiter for viscoelastic stress components. + +include $ASPECT_SOURCE_DIR/benchmarks/viscoelastic_bending_beam/viscoelastic_bending_beam.prm + +set Output directory = output_viscoelastic_bending_beam_weno + +subsection Discretization + subsection Stabilization parameters + # Use WENO limiter for viscoelastic stress components, while using the + # bound-preserving limiter for the chemical field + set Limiter for discontinuous composition solution = WENO, WENO, WENO, bound preserving + + # In this benchmark, the spurious oscillations at the boundaries of the beam + # are extremely serious, so we set the linear weight of neighbor cells one + # magnitude greater than the default value to make the reconstructed solution + # smoother. + set WENO linear weight of neighbor cells = 0.01 + + # The global maximum and minimum values are only used by the chemical field. + set Global composition maximum = 1.0 + set Global composition minimum = 0.0 + end +end + +# Visualize the KXRCF indicator for ve_stress_xy +subsection Postprocess + subsection Visualization + set List of output variables = material properties, strain rate, kxrcf indicator + + subsection KXRCF indicator + set Name of advection field = ve_stress_xy + end + end +end diff --git a/benchmarks/viscoelastic_plate_flexure/viscoelastic_plate_flexure.prm b/benchmarks/viscoelastic_plate_flexure/viscoelastic_plate_flexure.prm index 8b4f09583d9..ae3a3bc71a8 100644 --- a/benchmarks/viscoelastic_plate_flexure/viscoelastic_plate_flexure.prm +++ b/benchmarks/viscoelastic_plate_flexure/viscoelastic_plate_flexure.prm @@ -61,9 +61,9 @@ subsection Discretization set Use discontinuous composition discretization = true subsection Stabilization parameters - set Use limiter for discontinuous composition solution = true set Global composition maximum = 1.e11, 1.e11, 1.e11, 1.0, 1.0 set Global composition minimum = -1.e11, -1.e11, -1.e11, 0.0, 0.0 + set Limiter for discontinuous composition solution = bound preserving end end diff --git a/contrib/python/scripts/__pycache__/aspect_input.cpython-312.pyc b/contrib/python/scripts/__pycache__/aspect_input.cpython-312.pyc new file mode 100644 index 00000000000..f20772791ed Binary files /dev/null and b/contrib/python/scripts/__pycache__/aspect_input.cpython-312.pyc differ diff --git a/contrib/utilities/update_scripts/prm_files/reformat_limiter_settings.py b/contrib/utilities/update_scripts/prm_files/reformat_limiter_settings.py new file mode 100644 index 00000000000..92f11886d15 --- /dev/null +++ b/contrib/utilities/update_scripts/prm_files/reformat_limiter_settings.py @@ -0,0 +1,99 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +""" This script changes 'Use limiter for discontinuous temperature/composition solution' +to 'limiter for discontinuous temperature/composition solution' in order to allow users +to specify different limiters for different fields +""" + + +import sys +import os +import re +import argparse + +__author__ = 'The authors of the ASPECT code' +__copyright__ = 'Copyright 2024, ASPECT' +__license__ = 'GNU GPL 2 or later' + +# Add the ASPECT root directory to the path so we can import from the aspect_data module +sys.path.append(os.path.join(os.path.dirname(__file__), '..', '..', '..')) +import python.scripts.aspect_input as aspect + + + +def reformat_limiter_settings(parameters): + """change 'Use limiter for discontinuous temperature/composition solution' + to 'limiter for discontinuous temperature/composition solution' in order to allow users + to specify different limiters for different fields + """ + + old_entries = ["Use limiter for discontinuous temperature solution", + "Use limiter for discontinuous composition solution"] + new_entries = ["Limiter for discontinuous temperature solution", + "Limiter for discontinuous composition solution"] + + # Go to the subsection + if "Discretization" in parameters: + if "Stabilization parameters" in parameters["Discretization"]["value"]: + subsection = parameters["Discretization"]["value"]["Stabilization parameters"]["value"] + + for i in range(0,2): + if old_entries[i] in subsection: + + # Get the parameter value. Notice that the value may be followed by some + # extra comments, in which case we need to split them. + comments = subsection[old_entries[i]]["comment"] + values_and_extra_comments = subsection[old_entries[i]]["value"].split('#', 1) + values = values_and_extra_comments[0].rstrip() + extra_comments = '' if len(values_and_extra_comments) == 1 else " #" + values_and_extra_comments[1] + + # Delete the old entry + del subsection[old_entries[i]] + + new_values = '' + if i == 0: + # For temperature field, the parameter value is a boolean + new_values = "bound preserving" if values == "true" else "none" + + else: + # For compositional fields, the parameter value may be a boolean or + # a list of booleans. + split_values = values.split(',') + for i in range(0, len(split_values)): + if split_values[i] == "true": + new_values += "bound preserving" + else: + new_values += "none" + + if i < len(split_values)-1: + new_values += ", " + + # append the new values with the extra comments + new_values += extra_comments + subsection[new_entries[i]] = {"comment" : comments, + "value" : new_values, + "type" : "parameter", + "alignment spaces": 1} + + return parameters + + + +def main(input_file, output_file): + """Echo the input arguments to standard output""" + parameters = aspect.read_parameter_file(input_file) + parameters = reformat_limiter_settings(parameters) + aspect.write_parameter_file(parameters, output_file) + + + +if __name__ == '__main__': + parser = argparse.ArgumentParser( + prog='ASPECT .prm file reformatter', + description='Reformats ASPECT .prm files to follow our general formatting guidelines. See the documentation of this script for details.') + parser.add_argument('input_file', type=str, help='The .prm file to reformat') + parser.add_argument('output_file', type=str, help='The .prm file to write the reformatted file to') + args = parser.parse_args() + + sys.exit(main(args.input_file, args.output_file)) diff --git a/cookbooks/allken_et_al_2012_rift_interaction/allken.prm b/cookbooks/allken_et_al_2012_rift_interaction/allken.prm index a42b73d30c3..fbe89f631e1 100644 --- a/cookbooks/allken_et_al_2012_rift_interaction/allken.prm +++ b/cookbooks/allken_et_al_2012_rift_interaction/allken.prm @@ -5,7 +5,7 @@ set Dimension = 3 set Start time = 0 set End time = 5e6 set Use years in output instead of seconds = true -set Nonlinear solver scheme = iterated Stokes +set Nonlinear solver scheme = single Advection, iterated Stokes set Nonlinear solver tolerance = 1e-4 set Max nonlinear iterations = 25 set Output directory = output-allken_etal_2012 diff --git a/cookbooks/convection-box/tutorial-onset-of-convection/model_input/tutorial.prm b/cookbooks/convection-box/tutorial-onset-of-convection/model_input/tutorial.prm index 5c8302c747c..18e7f499a7b 100644 --- a/cookbooks/convection-box/tutorial-onset-of-convection/model_input/tutorial.prm +++ b/cookbooks/convection-box/tutorial-onset-of-convection/model_input/tutorial.prm @@ -82,21 +82,6 @@ end # whereas all other parts of the boundary are insulated (i.e., # no heat flux through these boundaries; this is also often used # to specify symmetry boundaries). -# subsection Model settings OLD DEPRECATED -# set Fixed temperature boundary indicators = 2,3 OLD MOVED - -# The next parameters then describe on which parts of the -# boundary we prescribe a zero or nonzero velocity and -# on which parts the flow is allowed to be tangential. -# Here, all four sides of the box allow tangential -# unrestricted flow but with a zero normal component: -# set Zero velocity boundary indicators = OLD MOVED -# set Prescribed velocity boundary indicators = OLD MOVED -# set Tangential velocity boundary indicators = 0,1,2,3 OLD MOVED - -subsection Boundary velocity model - set Tangential velocity boundary indicators = left, right, bottom, top -end # The final part of this section describes whether we # want to include adiabatic heating (from a small @@ -118,7 +103,11 @@ end # at the left and right does not matter.) subsection Boundary temperature model set Fixed temperature boundary indicators = 2,3 - set Model name = box + set List of model names = box + + # subsection Model settings OLD DEPRECATED + # set Fixed temperature boundary indicators = 2,3 OLD MOVED + subsection Box set Bottom temperature = 3600 @@ -126,6 +115,20 @@ subsection Boundary temperature model end end +subsection Boundary velocity model + set Tangential velocity boundary indicators = left, right, bottom, top + + # The next parameters then describe on which parts of the + # boundary we prescribe a zero or nonzero velocity and + # on which parts the flow is allowed to be tangential. + # Here, all four sides of the box allow tangential + # unrestricted flow but with a zero normal component: + # set Zero velocity boundary indicators = OLD MOVED + # set Prescribed velocity boundary indicators = OLD MOVED + # set Tangential velocity boundary indicators = 0,1,2,3 OLD MOVED + +end + # The final part is to specify what ASPECT should do with the # solution once computed at the end of every time step. The # process of evaluating the solution is called `postprocessing' diff --git a/cookbooks/future/mantle_setup_start.prm b/cookbooks/future/mantle_setup_start.prm index 4b733a7e240..98a79045800 100644 --- a/cookbooks/future/mantle_setup_start.prm +++ b/cookbooks/future/mantle_setup_start.prm @@ -7,7 +7,7 @@ set Dimension = 2 set Use years in output instead of seconds = true set End time = 1.5e8 # run until 1.5e8, then restart with Melt scaling factor threshold = 1e-8 set Output directory = output_ULVZ_2 -set Nonlinear solver scheme = iterated IMPES +set Nonlinear solver scheme = iterated Advection and Stokes set Max nonlinear iterations = 20 #set Resume computation = true @@ -101,7 +101,7 @@ end subsection Boundary composition model set Allow fixed composition on outflow boundaries = false set Fixed composition boundary indicators = top - set Model name = initial composition + set List of model names = initial composition end subsection Heating model diff --git a/cookbooks/inner_core_convection/inner_core_traction.prm b/cookbooks/inner_core_convection/inner_core_traction.prm index 215f28bc066..01137135c4c 100644 --- a/cookbooks/inner_core_convection/inner_core_traction.prm +++ b/cookbooks/inner_core_convection/inner_core_traction.prm @@ -131,6 +131,6 @@ end subsection Solver parameters subsection Stokes solver parameters - set Stokes solver type = block AMG + set Stokes solver type = block AMG end end diff --git a/cookbooks/kinematically_driven_subduction_2d/doc/Case1_postprocessing.prm b/cookbooks/kinematically_driven_subduction_2d/doc/Case1_postprocessing.prm index 13f6305edf3..c894b14bc71 100644 --- a/cookbooks/kinematically_driven_subduction_2d/doc/Case1_postprocessing.prm +++ b/cookbooks/kinematically_driven_subduction_2d/doc/Case1_postprocessing.prm @@ -1,5 +1,5 @@ subsection Postprocess - set List of postprocessors = visualization, velocity statistics, heating statistics, maximum depth of field, composition velocity statistics, viscous dissipation statistics, trench location + set List of postprocessors = visualization, velocity statistics, heating statistics, trench location subsection Visualization set List of output variables = density, viscosity, strain rate, error indicator diff --git a/cookbooks/kinematically_driven_subduction_2d/doc/Case2a_postprocessing.prm b/cookbooks/kinematically_driven_subduction_2d/doc/Case2a_postprocessing.prm index 97f6121cceb..9305efc0827 100644 --- a/cookbooks/kinematically_driven_subduction_2d/doc/Case2a_postprocessing.prm +++ b/cookbooks/kinematically_driven_subduction_2d/doc/Case2a_postprocessing.prm @@ -1,5 +1,5 @@ subsection Postprocess - set List of postprocessors = visualization, velocity statistics, heating statistics, maximum depth of field, composition velocity statistics, viscous dissipation statistics, temperature statistics, trench location, isotherm depth, composition RMS temperature statistics + set List of postprocessors = visualization, velocity statistics, heating statistics, temperature statistics, trench location, isotherm depth, composition RMS temperature statistics subsection Composition velocity statistics #Sum the velocities of the fields that make up the subducting plate. diff --git a/cookbooks/kinematically_driven_subduction_2d/kinematically_driven_subduction_2d_case1.prm b/cookbooks/kinematically_driven_subduction_2d/kinematically_driven_subduction_2d_case1.prm index ba095e61c31..ba634ce5f0d 100644 --- a/cookbooks/kinematically_driven_subduction_2d/kinematically_driven_subduction_2d_case1.prm +++ b/cookbooks/kinematically_driven_subduction_2d/kinematically_driven_subduction_2d_case1.prm @@ -174,7 +174,7 @@ subsection Formulation end subsection Postprocess - set List of postprocessors = visualization, velocity statistics, heating statistics, maximum depth of field, composition velocity statistics, viscous dissipation statistics, trench location + set List of postprocessors = visualization, velocity statistics, heating statistics, trench location subsection Visualization set List of output variables = density, viscosity, strain rate, error indicator diff --git a/cookbooks/kinematically_driven_subduction_2d/kinematically_driven_subduction_2d_case2a.prm b/cookbooks/kinematically_driven_subduction_2d/kinematically_driven_subduction_2d_case2a.prm index b7c767ab9e4..121df4b9d5c 100644 --- a/cookbooks/kinematically_driven_subduction_2d/kinematically_driven_subduction_2d_case2a.prm +++ b/cookbooks/kinematically_driven_subduction_2d/kinematically_driven_subduction_2d_case2a.prm @@ -86,7 +86,7 @@ subsection Initial temperature model end subsection Postprocess - set List of postprocessors = visualization, velocity statistics, heating statistics, maximum depth of field, composition velocity statistics, viscous dissipation statistics, temperature statistics, trench location, isotherm depth, composition RMS temperature statistics + set List of postprocessors = visualization, velocity statistics, heating statistics, temperature statistics, trench location, isotherm depth, composition RMS temperature statistics subsection Composition velocity statistics #Sum the velocities of the fields that make up the subducting plate. diff --git a/cookbooks/magnetic_stripes/magnetic_stripes.cc b/cookbooks/magnetic_stripes/magnetic_stripes.cc index 5a3cd1340f2..0c13a2b6ef9 100644 --- a/cookbooks/magnetic_stripes/magnetic_stripes.cc +++ b/cookbooks/magnetic_stripes/magnetic_stripes.cc @@ -63,7 +63,7 @@ namespace aspect break; } - for (unsigned int i=0; i < in.position.size(); ++i) + for (unsigned int i=0; i < in.n_evaluation_points(); ++i) { const double depth = this->get_geometry_model().depth(in.position[i]); const double reaction_depth = 7000.0; diff --git a/cookbooks/magnetic_stripes/magnetic_stripes.prm b/cookbooks/magnetic_stripes/magnetic_stripes.prm index 4dc64baec60..feed74fb113 100644 --- a/cookbooks/magnetic_stripes/magnetic_stripes.prm +++ b/cookbooks/magnetic_stripes/magnetic_stripes.prm @@ -17,9 +17,9 @@ subsection Discretization set Use discontinuous composition discretization = true subsection Stabilization parameters - set Use limiter for discontinuous composition solution = true set Global composition maximum = 1 set Global composition minimum = -1 + set Limiter for discontinuous composition solution = bound preserving end end @@ -107,7 +107,7 @@ end subsection Boundary composition model set Fixed composition boundary indicators = bottom - set Model name = initial composition + set List of model names = initial composition end subsection Material model diff --git a/cookbooks/morency_doin_2004/morency_doin.prm b/cookbooks/morency_doin_2004/morency_doin.prm index 40880818d7c..951205a24ca 100644 --- a/cookbooks/morency_doin_2004/morency_doin.prm +++ b/cookbooks/morency_doin_2004/morency_doin.prm @@ -103,7 +103,7 @@ end subsection Solver parameters subsection Stokes solver parameters - set Stokes solver type = block AMG + set Stokes solver type = block AMG end end diff --git a/cookbooks/tomography_based_plate_motions/plugins/tomography_based_plate_motions.cc b/cookbooks/tomography_based_plate_motions/plugins/tomography_based_plate_motions.cc index d3c61e9a379..8a9fde17252 100644 --- a/cookbooks/tomography_based_plate_motions/plugins/tomography_based_plate_motions.cc +++ b/cookbooks/tomography_based_plate_motions/plugins/tomography_based_plate_motions.cc @@ -1635,7 +1635,7 @@ namespace aspect // reduce grain size. if (out.template get_additional_output>() == nullptr) { - const unsigned int n_points = out.viscosities.size(); + const unsigned int n_points = out.n_evaluation_points(); out.additional_outputs.push_back( std::make_unique> (n_points)); } @@ -1643,7 +1643,7 @@ namespace aspect // We need the prescribed field outputs to interpolate the grain size onto a compositional field. if (out.template get_additional_output>() == nullptr) { - const unsigned int n_points = out.viscosities.size(); + const unsigned int n_points = out.n_evaluation_points(); out.additional_outputs.push_back( std::make_unique> (n_points, this->n_compositional_fields())); } @@ -1652,7 +1652,7 @@ namespace aspect if (use_table_properties && out.template get_additional_output>() == nullptr) { - const unsigned int n_points = out.viscosities.size(); + const unsigned int n_points = out.n_evaluation_points(); out.additional_outputs.push_back( std::make_unique> (n_points)); } @@ -1660,7 +1660,7 @@ namespace aspect if (this->get_parameters().temperature_method == Parameters::AdvectionFieldMethod::prescribed_field && out.template get_additional_output>() == nullptr) { - const unsigned int n_points = out.viscosities.size(); + const unsigned int n_points = out.n_evaluation_points(); out.additional_outputs.push_back( std::make_unique> (n_points)); } @@ -1668,14 +1668,14 @@ namespace aspect // We need additional field outputs for the unscaled viscosity if (out.template get_additional_output>() == nullptr) { - const unsigned int n_points = out.viscosities.size(); + const unsigned int n_points = out.n_evaluation_points(); out.additional_outputs.push_back( std::make_unique> (n_points)); } if (out.template get_additional_output>() == nullptr) { - const unsigned int n_points = out.viscosities.size(); + const unsigned int n_points = out.n_evaluation_points(); out.additional_outputs.push_back( std::make_unique> (n_points)); } diff --git a/doc/manual/citing_aspect.bib b/doc/manual/citing_aspect.bib index 1ef12192db7..f332ea06521 100644 --- a/doc/manual/citing_aspect.bib +++ b/doc/manual/citing_aspect.bib @@ -2068,4 +2068,32 @@ @article{SAKI2024107212 keywords = {Mantle transition zone discontinuities, PP and SS precursors, Caribbean subduction, Plume-slab interaction} } +@article{zhong2013, +title = {A simple weighted essentially nonoscillatory limiter for {Runge}–{Kutta} discontinuous {Galerkin} methods}, +volume = {232}, +issn = {0021-9991}, +doi = {10.1016/j.jcp.2012.08.028}, +language = {en}, +number = {1}, +journal = {Journal of Computational Physics}, +author = {Zhong, Xinghui and Shu, Chi-Wang}, +month = jan, +year = {2013}, +keywords = {Discontinuous Galerkin method, WENO limiter}, +pages = {397--415} +} +@article{krivodonova2004, +series = {Workshop on {Innovative} {Time} {Integrators} for {PDEs}}, +title = {Shock detection and limiting with discontinuous {Galerkin} methods for hyperbolic conservation laws}, +volume = {48}, +issn = {0168-9274}, +doi = {10.1016/j.apnum.2003.11.002}, +language = {en}, +number = {3}, +journal = {Applied Numerical Mathematics}, +author = {Krivodonova, L. and Xin, J. and Remacle, J. -F. and Chevaugeon, N. and Flaherty, J. E.}, +month = mar, +year = {2004}, +pages = {323--338} +} diff --git a/include/aspect/parameters.h b/include/aspect/parameters.h index 201fb425efa..baabeb3865a 100644 --- a/include/aspect/parameters.h +++ b/include/aspect/parameters.h @@ -313,6 +313,41 @@ namespace aspect } }; + /** + * A struct to provide the different choices of limiters for + * the discontinuous Galerkin method. + */ + struct DGLimiterType + { + enum Kind + { + bound_preserving, + WENO, + none + }; + + static + Kind + parse(const std::string &input) + { + if (input == "bound preserving") + return bound_preserving; + else if (input == "WENO") + return WENO; + else if (input == "none") + return none; + else + AssertThrow(false, ExcNotImplemented()); + + return Kind(); + } + + static const std::string pattern() + { + return "bound preserving|WENO|none"; + } + }; + /** * This enum represents the different choices for the linear solver * for the Stoke system. See @p stokes_solver_type. @@ -639,7 +674,7 @@ namespace aspect std::vector additional_refinement_times; unsigned int adaptive_refinement_interval; bool skip_solvers_on_initial_refinement; - bool skip_setup_initial_conditions_on_initial_refinement; + bool skip_setup_initial_temperature_on_initial_refinement; bool run_postprocessors_on_initial_refinement; bool run_postprocessors_on_nonlinear_iterations; /** @@ -656,12 +691,15 @@ namespace aspect std::vector stabilization_beta; double stabilization_gamma; double discontinuous_penalty; - bool use_limiter_for_discontinuous_temperature_solution; - std::vector use_limiter_for_discontinuous_composition_solution; + typename DGLimiterType::Kind limiter_for_discontinuous_temperature_solution; + std::vector limiter_for_discontinuous_composition_solution; double global_temperature_max_preset; double global_temperature_min_preset; std::vector global_composition_max_preset; std::vector global_composition_min_preset; + double temperature_KXRCF_indicator_threshold; + std::vector composition_KXRCF_indicator_threshold; + double WENO_linear_weight; std::vector compositional_fields_with_disabled_boundary_entropy_viscosity; /** diff --git a/include/aspect/particle/property/crystal_preferred_orientation.h b/include/aspect/particle/property/crystal_preferred_orientation.h index 1c10dcc087b..29c038f87c8 100644 --- a/include/aspect/particle/property/crystal_preferred_orientation.h +++ b/include/aspect/particle/property/crystal_preferred_orientation.h @@ -121,7 +121,7 @@ namespace aspect * * We store the same number of grains for all minerals (e.g. olivine and enstatite * grains), although their volume fractions may not be the same. This is because we need a minimum number - * of grains per tracer to perform reliable statistics on it. This minimum should be the same for all + * of grains per particle to perform reliable statistics on it. This minimum should be the same for all * minerals. * * @ingroup ParticleProperties diff --git a/include/aspect/postprocess/visualization/kxrcf_indicator.h b/include/aspect/postprocess/visualization/kxrcf_indicator.h new file mode 100644 index 00000000000..bcd4022612d --- /dev/null +++ b/include/aspect/postprocess/visualization/kxrcf_indicator.h @@ -0,0 +1,80 @@ +/* + Copyright (C) 2011 - 2022 by the authors of the ASPECT code. + + This file is part of ASPECT. + + ASPECT is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2, or (at your option) + any later version. + + ASPECT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ASPECT; see the file LICENSE. If not see + . +*/ + + +#ifndef _aspect_postprocess_visualization_kxrcf_indicator_h +#define _aspect_postprocess_visualization_kxrcf_indicator_h + +#include +#include + + +namespace aspect +{ + namespace Postprocess + { + namespace VisualizationPostprocessors + { + /** + * A class derived that implements a function that provides the + * KXRCF indicator for a given advection field on each cell for + * graphical output. + */ + template + class KXRCFIndicator : public CellDataVectorCreator, + public SimulatorAccess + { + public: + /** + * Constructor. + */ + KXRCFIndicator(); + + /** + * @copydoc CellDataVectorCreator::execute() + */ + std::pair>> + execute () const override; + + /** + * Declare the parameters this class takes through input files. + */ + static + void + declare_parameters (ParameterHandler &prm); + + /** + * Read the parameters this class declares from the parameter file. + */ + void + parse_parameters (ParameterHandler &prm) override; + + private: + /** + * A parameter that tells us for which advection field the + * KXRCF indicator should be visualized. + */ + unsigned int field_index; + }; + } + } +} + +#endif diff --git a/include/aspect/simulator.h b/include/aspect/simulator.h index 3df4a6752a9..0c7fbd1630f 100644 --- a/include/aspect/simulator.h +++ b/include/aspect/simulator.h @@ -1373,10 +1373,37 @@ namespace aspect * the discontinuous Galerkin solution will stay in the prescribed bounds. * * This function is implemented in - * source/simulator/helper_functions.cc. + * source/simulator/limiters.cc. + */ + void apply_BP_limiter_to_dg_solutions (const AdvectionField &advection_field); + + /** + * Apply the WENO limiter to the discontinuous Galerkin solutions. + * WENO is short for Weighted Essentially Non-Oscillatory, which is + * a class of high-resolution schemes used in the numerical solution + * of hyperbolic conservation laws. The basic idea of WENO limiter is + * to replace the solution in troubled cells (cells with oscillation) + * by a polynomial reconstruction that takes the neighbor cells into + * account. The WENO scheme implemented in the program is a simple + * variant proposed by Zhong and Shu, 2013. + * + * This function is implemented in + * source/simulator/limiters.cc. */ - void apply_limiter_to_dg_solutions (const AdvectionField &advection_field); + void apply_WENO_limiter_to_dg_solutions (const AdvectionField &advection_field); + /** + * Fills a vector with the KXRCF indicator for a given advection field + * on each local cell. The KXRCF indicator is a metric of discontinuity + * for hyperbolic conservation laws. Cells with high KXRCF values are + * identified as "troubled cells" and will be smoothed by the WENO limiter. + * + * This function is implemented in + * source/simulator/limiters.cc. + */ + template + void compute_KXRCF_indicators(Vector &KXRCF_indicators, + const AdvectionField &advection_field) const; /** * Compute the reactions in case of operator splitting: diff --git a/include/aspect/simulator_access.h b/include/aspect/simulator_access.h index a65c0cf4879..999a22a6980 100644 --- a/include/aspect/simulator_access.h +++ b/include/aspect/simulator_access.h @@ -427,6 +427,14 @@ namespace aspect void get_artificial_viscosity_composition(Vector &viscosity_per_cell, const unsigned int compositional_variable) const; + + /** + * Compute the KXRCF indicator on each locally owned cell for a specific + * advection field. + */ + void + compute_KXRCF_indicators(Vector &KXRCF_indicators, + const unsigned int field_index) const; /** @} */ diff --git a/include/aspect/volume_of_fluid/assembly.h b/include/aspect/volume_of_fluid/assembly.h index bc212c798e2..fb46191c864 100644 --- a/include/aspect/volume_of_fluid/assembly.h +++ b/include/aspect/volume_of_fluid/assembly.h @@ -18,8 +18,8 @@ . */ -#ifndef _aspect_volume_of_fluid_assembly_h -#define _aspect_volume_of_fluid_assembly_h +#ifndef _aspect_volume_of_fluid_simulator/assemblers/interface.h +#define _aspect_volume_of_fluid_simulator/assemblers/interface.h #include #include diff --git a/include/aspect/volume_of_fluid/handler.h b/include/aspect/volume_of_fluid/handler.h index 56be01c21f6..8e741060777 100644 --- a/include/aspect/volume_of_fluid/handler.h +++ b/include/aspect/volume_of_fluid/handler.h @@ -24,7 +24,7 @@ #include #include #include -#include +#include #include diff --git a/source/material_model/melt_boukare.cc b/source/material_model/melt_boukare.cc index 2e0b0ab91c2..4022c8161ae 100644 --- a/source/material_model/melt_boukare.cc +++ b/source/material_model/melt_boukare.cc @@ -508,7 +508,7 @@ namespace aspect const unsigned int n_endmembers = endmember_names.size(); EndmemberProperties endmembers(n_endmembers); - for (unsigned int q=0; q endmember_mole_fractions_per_phase(n_endmembers); @@ -627,7 +627,7 @@ namespace aspect const unsigned int n_endmembers = endmember_names.size(); EndmemberProperties endmembers(n_endmembers); - for (unsigned int q=0; q endmember_mole_fractions_per_phase(n_endmembers); std::vector endmember_mole_fractions_in_composite(n_endmembers); diff --git a/source/material_model/reactive_fluid_transport.cc b/source/material_model/reactive_fluid_transport.cc index 84bd54dc09f..1f2d76610bb 100644 --- a/source/material_model/reactive_fluid_transport.cc +++ b/source/material_model/reactive_fluid_transport.cc @@ -123,7 +123,7 @@ namespace aspect melt_fractions (const MaterialModel::MaterialModelInputs &in, std::vector &melt_fractions) const { - for (unsigned int q=0; qintrospection().compositional_index_for_name("porosity"); switch (fluid_solid_reaction_scheme) diff --git a/source/particle/property/crystal_preferred_orientation.cc b/source/particle/property/crystal_preferred_orientation.cc index 3f849f0ca4d..58aab2dfdcb 100644 --- a/source/particle/property/crystal_preferred_orientation.cc +++ b/source/particle/property/crystal_preferred_orientation.cc @@ -139,7 +139,7 @@ namespace aspect // // Note that we store exactly the same number of grains of all minerals (e.g. olivine and enstatite // grains), although their volume fractions may not be the same. We need a minimum amount - // of grains per tracer to perform reliable statistics on it. This minimum is the same for all phases. + // of grains per particle to perform reliable statistics on it. This minimum is the same for all phases. // and enstatite. // // Furthermore, for this plugin the following dims are always 3. When using 2d an infinitely thin 3d domain is assumed. diff --git a/source/postprocess/visualization/kxrcf_indicator.cc b/source/postprocess/visualization/kxrcf_indicator.cc new file mode 100644 index 00000000000..b932f96d09f --- /dev/null +++ b/source/postprocess/visualization/kxrcf_indicator.cc @@ -0,0 +1,135 @@ +/* + Copyright (C) 2011 - 2022 by the authors of the ASPECT code. + + This file is part of ASPECT. + + ASPECT is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2, or (at your option) + any later version. + + ASPECT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ASPECT; see the file LICENSE. If not see + . +*/ + + +#include + +namespace aspect +{ + namespace Postprocess + { + namespace VisualizationPostprocessors + { + template + KXRCFIndicator::KXRCFIndicator() + : + CellDataVectorCreator("") + {} + + + template + std::pair>> + KXRCFIndicator::execute() const + { + std::pair>> + return_value ("KXRCF_indicator", + std::make_unique>(this->get_triangulation().n_active_cells())); + this->compute_KXRCF_indicators(*return_value.second, field_index); + return return_value; + } + + + template + void + KXRCFIndicator::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Postprocess"); + { + prm.enter_subsection("Visualization"); + { + prm.enter_subsection("KXRCF indicator"); + { + prm.declare_entry ("Name of advection field", "", + Patterns::Anything(), + "The name of the advection field whose output " + "should be visualized. "); + } + prm.leave_subsection(); + } + prm.leave_subsection(); + } + prm.leave_subsection(); + } + + + template + void + KXRCFIndicator::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Postprocess"); + { + prm.enter_subsection("Visualization"); + { + prm.enter_subsection("KXRCF indicator"); + { + const std::string field_name = prm.get("Name of advection field"); + + if (field_name == "temperature") + field_index = 0; + else + { + AssertThrow(this->introspection().compositional_name_exists(field_name), + ExcMessage("No compositional field with name <" + + field_name + + "> exists for which you want to visualize the KXRCF indicator.")); + + field_index = this->introspection().compositional_index_for_name(field_name) + 1; + } + } + prm.leave_subsection(); + } + prm.leave_subsection(); + } + prm.leave_subsection(); + } + } + } +} + + +// explicit instantiations +namespace aspect +{ + namespace Postprocess + { + namespace VisualizationPostprocessors + { + ASPECT_REGISTER_VISUALIZATION_POSTPROCESSOR(KXRCFIndicator, + "kxrcf indicator", + "A visualization output object that generates output " + "showing the value of the KXRCF indicator for a given " + "advection field (either temperature or compositional " + "field) on each cell. The KXRCF indicator is a metric " + "of discontinuity for hyperbolic conservation laws. " + "If the KXRCF value of a cell is higher than the value of " + "``Temperature/Composition KXRCF indicator threshold'' " + "of the corresponding field in the input parameter file, " + "the cell will be identified as a ``troubled-cell'' that " + "should be smoothed by the WENO limiter. For details, see " + "\\cite{Krivodonova:etal:2004}." + "\n\n" + "This postprocessor should only be used for discontinuous " + "advection fields. Otherwise, the postprocessor will produce " + "a meaningless visualization output." + "\n\n" + "Physical units: none.") + } + } +} diff --git a/source/simulator/core.cc b/source/simulator/core.cc index d74c564d1ea..9b1f7b535fb 100644 --- a/source/simulator/core.cc +++ b/source/simulator/core.cc @@ -485,14 +485,16 @@ namespace aspect if (MappingQCache *map = dynamic_cast*>(&(*mapping))) map->initialize(MappingQGeneric(4), triangulation); - bool dg_limiter_enabled = parameters.use_limiter_for_discontinuous_temperature_solution; + bool bp_limiter_enabled = (parameters.limiter_for_discontinuous_temperature_solution + == Parameters::DGLimiterType::bound_preserving); for (unsigned int c=0; c::DGLimiterType::bound_preserving); - // Check that DG limiters are only used with cartesian mapping - if (dg_limiter_enabled) + // Check that BP limiters are only used with cartesian mapping + if (bp_limiter_enabled) AssertThrow(geometry_model->natural_coordinate_system() == Utilities::Coordinates::CoordinateSystem::cartesian, - ExcMessage("The limiter for the discontinuous temperature and composition solutions " + ExcMessage("The bound preserving limiter for the discontinuous temperature and composition solutions " "has not been tested in non-Cartesian geometries and currently requires " "the use of a Cartesian geometry model.")); @@ -2049,7 +2051,7 @@ namespace aspect timestep_number = 0; time_step = old_time_step = 0; - if (! parameters.skip_setup_initial_conditions_on_initial_refinement + if (! parameters.skip_setup_initial_temperature_on_initial_refinement || ! (pre_refinement_step < parameters.initial_adaptive_refinement)) { diff --git a/source/simulator/helper_functions.cc b/source/simulator/helper_functions.cc index 0450a1de14d..bd84cabf160 100644 --- a/source/simulator/helper_functions.cc +++ b/source/simulator/helper_functions.cc @@ -1392,248 +1392,6 @@ namespace aspect - template - void Simulator::apply_limiter_to_dg_solutions (const AdvectionField &advection_field) - { - // TODO: Modify to more robust method - // Skip if this composition field is being set from the volume_of_fluid handler - if (!advection_field.is_temperature() && - parameters.volume_of_fluid_tracking_enabled) - if (volume_of_fluid_handler->field_index_for_name(introspection.name_for_compositional_index(advection_field.compositional_variable)) - != volume_of_fluid_handler->get_n_fields()) - return; - - /* - * First setup the quadrature points which are used to find the maximum and minimum solution values at those points. - * A quadrature formula that combines all quadrature points constructed as all tensor products of - * 1) one dimensional Gauss points; 2) one dimensional Gauss-Lobatto points. - * We require that the Gauss-Lobatto points (2) appear in only one direction. - * Therefore, possible combination - * in 2d: the combinations are 21, 12 - * in 3d: the combinations are 211, 121, 112 - */ - const QGauss<1> quadrature_formula_1 (advection_field.polynomial_degree(introspection)+1); - const QGaussLobatto<1> quadrature_formula_2 (advection_field.polynomial_degree(introspection)+1); - - const unsigned int n_q_points_1 = quadrature_formula_1.size(); - const unsigned int n_q_points_2 = quadrature_formula_2.size(); - const unsigned int n_q_points = dim * n_q_points_2 * static_cast(std::pow(n_q_points_1, dim-1)); - - std::vector> quadrature_points; - quadrature_points.reserve(n_q_points); - - switch (dim) - { - case 2: - { - // append quadrature points combination 12 - for (unsigned int i=0; i < n_q_points_1 ; ++i) - { - const double x = quadrature_formula_1.point(i)(0); - for (unsigned int j=0; j < n_q_points_2 ; ++j) - { - const double y = quadrature_formula_2.point(j)(0); - quadrature_points.push_back(Point(x,y)); - } - } - // append quadrature points combination 21 - for (unsigned int i=0; i < n_q_points_2 ; ++i) - { - const double x = quadrature_formula_2.point(i)(0); - for (unsigned int j=0; j < n_q_points_1 ; ++j) - { - const double y = quadrature_formula_1.point(j)(0); - quadrature_points.push_back(Point(x,y)); - } - } - break; - } - - case 3: - { - // append quadrature points combination 121 - for ( unsigned int i=0; i < n_q_points_1 ; ++i) - { - const double x = quadrature_formula_1.point(i)(0); - for ( unsigned int j=0; j < n_q_points_2 ; ++j) - { - const double y = quadrature_formula_2.point(j)(0); - for ( unsigned int k=0; k < n_q_points_1 ; ++k) - { - const double z = quadrature_formula_1.point(k)(0); - quadrature_points.push_back(Point(x,y,z)); - } - } - } - // append quadrature points combination 112 - for (unsigned int i=0; i < n_q_points_1 ; ++i) - { - const double x = quadrature_formula_1.point(i)(0); - for (unsigned int j=0; j < n_q_points_1 ; ++j) - { - const double y = quadrature_formula_1.point(j)(0); - for (unsigned int k=0; k < n_q_points_2 ; ++k) - { - const double z = quadrature_formula_2.point(k)(0); - quadrature_points.push_back(Point(x,y,z)); - } - } - } - // append quadrature points combination 211 - for (unsigned int i=0; i < n_q_points_2 ; ++i) - { - const double x = quadrature_formula_2.point(i)(0); - for ( unsigned int j=0; j < n_q_points_1 ; ++j) - { - const double y = quadrature_formula_1.point(j)(0); - for ( unsigned int k=0; k < n_q_points_1 ; ++k) - { - const double z = quadrature_formula_1.point(k)(0); - quadrature_points.push_back(Point(x,y,z)); - } - } - } - break; - } - - default: - Assert (false, ExcNotImplemented()); - } - - Assert (quadrature_points.size() == n_q_points, ExcInternalError()); - const Quadrature quadrature_formula(quadrature_points); - - // Quadrature rules only used for the numerical integration for better accuracy - const Quadrature &quadrature_formula_0 - = (advection_field.is_temperature() ? - introspection.quadratures.temperature : - introspection.quadratures.compositional_fields); - const unsigned int n_q_points_0 = quadrature_formula_0.size(); - - // fe values for points evaluation - FEValues fe_values (*mapping, - finite_element, - quadrature_formula, - update_values | - update_quadrature_points); - std::vector values (n_q_points); - // fe values for numerical integration, with a number of quadrature points - // that is equal to 1/dim times the number of total points above - FEValues fe_values_0 (*mapping, - finite_element, - quadrature_formula_0, - update_values | - update_quadrature_points | - update_JxW_values); - std::vector values_0 (n_q_points_0); - - const FEValuesExtractors::Scalar field - = (advection_field.is_temperature() - ? - introspection.extractors.temperature - : - introspection.extractors.compositional_fields[advection_field.compositional_variable] - ); - - const double max_solution_exact_global = (advection_field.is_temperature() - ? - parameters.global_temperature_max_preset - : - parameters.global_composition_max_preset[advection_field.compositional_variable] - ); - const double min_solution_exact_global = (advection_field.is_temperature() - ? - parameters.global_temperature_min_preset - : - parameters.global_composition_min_preset[advection_field.compositional_variable] - ); - - LinearAlgebra::BlockVector distributed_solution (introspection.index_sets.system_partitioning, - mpi_communicator); - const unsigned int block_idx = advection_field.block_index(introspection); - distributed_solution.block(block_idx) = solution.block(block_idx); - - std::vector local_dof_indices (finite_element.dofs_per_cell); - - for (const auto &cell : dof_handler.active_cell_iterators()) - if (cell->is_locally_owned()) - { - cell->get_dof_indices (local_dof_indices); - // used to find the maximum, minimum - fe_values.reinit (cell); - fe_values[field].get_function_values(solution, values); - // used for the numerical integration - fe_values_0.reinit (cell); - fe_values_0[field].get_function_values(solution, values_0); - - // Find the local max and local min - const double min_solution_local = *std::min_element (values.begin(), values.end()); - const double max_solution_local = *std::max_element (values.begin(), values.end()); - // Find the trouble cell - if (min_solution_local < min_solution_exact_global - || max_solution_local > max_solution_exact_global) - { - // Compute the cell area and cell solution average - double local_area = 0.0; - double local_solution_average = 0.0; - for (unsigned int q = 0; q < n_q_points_0; ++q) - { - local_area += fe_values_0.JxW(q); - local_solution_average += values_0[q]*fe_values_0.JxW(q); - } - local_solution_average /= local_area; - - /* - * Define theta: a scaling constant used to correct the old solution by the formula - * new_value = theta * (old_value-old_solution_cell_average)+old_solution_cell_average - * where theta \in [0,1] defined as below. - * After the correction, the new solution does not exceed the user-given - * exact global maximum/minimum values. Meanwhile, the new solution's cell average - * equals to the old solution's cell average. - */ - double theta = 1.0; - if (std::abs(max_solution_local-local_solution_average) > std::numeric_limits::min()) - { - theta = std::min(theta, std::abs((max_solution_exact_global-local_solution_average) - / (max_solution_local-local_solution_average))); - } - if (std::abs(min_solution_local-local_solution_average) > std::numeric_limits::min()) - { - theta = std::min(theta, std::abs((min_solution_exact_global-local_solution_average) - / (min_solution_local-local_solution_average))); - } - - /* Modify the advection degrees of freedom of the numerical solution. - * Note that we are using DG elements, so every DoF on a locally owned cell is locally owned; - * this means that we do not need to check whether the 'distributed_solution' vector actually - * stores the element we read from/write to here. - */ - for (unsigned int j = 0; - j < finite_element.base_element(advection_field.base_element(introspection)).dofs_per_cell; - ++j) - { - const unsigned int support_point_index = finite_element.component_to_system_index( - (advection_field.is_temperature() - ? - introspection.component_indices.temperature - : - introspection.component_indices.compositional_fields[advection_field.compositional_variable] - ), - /*dof index within component=*/ j); - const double solution_value = solution(local_dof_indices[support_point_index]); - const double limited_solution_value = theta * (solution_value-local_solution_average) + local_solution_average; - distributed_solution(local_dof_indices[support_point_index]) = limited_solution_value; - } - } - } - - distributed_solution.compress(VectorOperation::insert); - // now get back to the original vector - solution.block(block_idx) = distributed_solution.block(block_idx); - } - - - template void Simulator::update_solution_vectors_with_reaction_results (const unsigned int block_index, const LinearAlgebra::BlockVector &distributed_vector, @@ -2740,7 +2498,6 @@ namespace aspect template bool Simulator::stokes_matrix_depends_on_solution() const; \ template bool Simulator::stokes_A_block_is_symmetric() const; \ template void Simulator::interpolate_onto_velocity_system(const TensorFunction<1,dim> &func, LinearAlgebra::Vector &vec);\ - template void Simulator::apply_limiter_to_dg_solutions(const AdvectionField &advection_field); \ template void Simulator::compute_reactions(); \ template void Simulator::initialize_current_linearization_point (); \ template void Simulator::interpolate_material_output_into_advection_field(const AdvectionField &adv_field); \ diff --git a/source/simulator/limiters.cc b/source/simulator/limiters.cc new file mode 100644 index 00000000000..bce8fba5b4a --- /dev/null +++ b/source/simulator/limiters.cc @@ -0,0 +1,876 @@ +#include +#include +#include + +namespace aspect +{ + template + void + Simulator::apply_BP_limiter_to_dg_solutions(const AdvectionField &advection_field) + { + // TODO: Modify to more robust method + // Skip if this composition field is being set from the volume_of_fluid handler + if (!advection_field.is_temperature() && + parameters.volume_of_fluid_tracking_enabled) + if (volume_of_fluid_handler->field_index_for_name(introspection.name_for_compositional_index(advection_field.compositional_variable)) + != volume_of_fluid_handler->get_n_fields()) + return; + + /* + * First setup the quadrature points which are used to find the maximum and minimum solution values at those points. + * A quadrature formula that combines all quadrature points constructed as all tensor products of + * 1) one dimensional Gauss points; 2) one dimensional Gauss-Lobatto points. + * We require that the Gauss-Lobatto points (2) appear in only one direction. + * Therefore, possible combination + * in 2d: the combinations are 21, 12 + * in 3d: the combinations are 211, 121, 112 + */ + const QGauss<1> quadrature_formula_1 (advection_field.polynomial_degree(introspection)+1); + const QGaussLobatto<1> quadrature_formula_2 (advection_field.polynomial_degree(introspection)+1); + + const unsigned int n_q_points_1 = quadrature_formula_1.size(); + const unsigned int n_q_points_2 = quadrature_formula_2.size(); + const unsigned int n_q_points = dim * n_q_points_2 * static_cast(std::pow(n_q_points_1, dim-1)); + + std::vector> quadrature_points; + quadrature_points.reserve(n_q_points); + + switch (dim) + { + case 2: + { + // append quadrature points combination 12 + for (unsigned int i=0; i < n_q_points_1 ; ++i) + { + const double x = quadrature_formula_1.point(i)(0); + for (unsigned int j=0; j < n_q_points_2 ; ++j) + { + const double y = quadrature_formula_2.point(j)(0); + quadrature_points.push_back(Point(x,y)); + } + } + // append quadrature points combination 21 + for (unsigned int i=0; i < n_q_points_2 ; ++i) + { + const double x = quadrature_formula_2.point(i)(0); + for (unsigned int j=0; j < n_q_points_1 ; ++j) + { + const double y = quadrature_formula_1.point(j)(0); + quadrature_points.push_back(Point(x,y)); + } + } + break; + } + + case 3: + { + // append quadrature points combination 121 + for ( unsigned int i=0; i < n_q_points_1 ; ++i) + { + const double x = quadrature_formula_1.point(i)(0); + for ( unsigned int j=0; j < n_q_points_2 ; ++j) + { + const double y = quadrature_formula_2.point(j)(0); + for ( unsigned int k=0; k < n_q_points_1 ; ++k) + { + const double z = quadrature_formula_1.point(k)(0); + quadrature_points.push_back(Point(x,y,z)); + } + } + } + // append quadrature points combination 112 + for (unsigned int i=0; i < n_q_points_1 ; ++i) + { + const double x = quadrature_formula_1.point(i)(0); + for (unsigned int j=0; j < n_q_points_1 ; ++j) + { + const double y = quadrature_formula_1.point(j)(0); + for (unsigned int k=0; k < n_q_points_2 ; ++k) + { + const double z = quadrature_formula_2.point(k)(0); + quadrature_points.push_back(Point(x,y,z)); + } + } + } + // append quadrature points combination 211 + for (unsigned int i=0; i < n_q_points_2 ; ++i) + { + const double x = quadrature_formula_2.point(i)(0); + for ( unsigned int j=0; j < n_q_points_1 ; ++j) + { + const double y = quadrature_formula_1.point(j)(0); + for ( unsigned int k=0; k < n_q_points_1 ; ++k) + { + const double z = quadrature_formula_1.point(k)(0); + quadrature_points.push_back(Point(x,y,z)); + } + } + } + break; + } + + default: + Assert (false, ExcNotImplemented()); + } + + Assert (quadrature_points.size() == n_q_points, ExcInternalError()); + const Quadrature quadrature_formula(quadrature_points); + + // Quadrature rules only used for the numerical integration for better accuracy + const Quadrature &quadrature_formula_0 + = (advection_field.is_temperature() ? + introspection.quadratures.temperature : + introspection.quadratures.compositional_fields); + const unsigned int n_q_points_0 = quadrature_formula_0.size(); + + // fe values for points evaluation + FEValues fe_values (*mapping, + finite_element, + quadrature_formula, + update_values | + update_quadrature_points); + std::vector values (n_q_points); + // fe values for numerical integration, with a number of quadrature points + // that is equal to 1/dim times the number of total points above + FEValues fe_values_0 (*mapping, + finite_element, + quadrature_formula_0, + update_values | + update_quadrature_points | + update_JxW_values); + std::vector values_0 (n_q_points_0); + + const FEValuesExtractors::Scalar field + = (advection_field.is_temperature() + ? + introspection.extractors.temperature + : + introspection.extractors.compositional_fields[advection_field.compositional_variable] + ); + + const double max_solution_exact_global = (advection_field.is_temperature() + ? + parameters.global_temperature_max_preset + : + parameters.global_composition_max_preset[advection_field.compositional_variable] + ); + const double min_solution_exact_global = (advection_field.is_temperature() + ? + parameters.global_temperature_min_preset + : + parameters.global_composition_min_preset[advection_field.compositional_variable] + ); + + LinearAlgebra::BlockVector distributed_solution (introspection.index_sets.system_partitioning, + mpi_communicator); + const unsigned int block_idx = advection_field.block_index(introspection); + distributed_solution.block(block_idx) = solution.block(block_idx); + + std::vector local_dof_indices (finite_element.dofs_per_cell); + + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + { + cell->get_dof_indices (local_dof_indices); + // used to find the maximum, minimum + fe_values.reinit (cell); + fe_values[field].get_function_values(solution, values); + // used for the numerical integration + fe_values_0.reinit (cell); + fe_values_0[field].get_function_values(solution, values_0); + + // Find the local max and local min + const double min_solution_local = *std::min_element (values.begin(), values.end()); + const double max_solution_local = *std::max_element (values.begin(), values.end()); + // Find the trouble cell + if (min_solution_local < min_solution_exact_global + || max_solution_local > max_solution_exact_global) + { + // Compute the cell area and cell solution average + double local_area = 0.0; + double local_solution_average = 0.0; + for (unsigned int q = 0; q < n_q_points_0; ++q) + { + local_area += fe_values_0.JxW(q); + local_solution_average += values_0[q]*fe_values_0.JxW(q); + } + local_solution_average /= local_area; + + /* + * Define theta: a scaling constant used to correct the old solution by the formula + * new_value = theta * (old_value-old_solution_cell_average)+old_solution_cell_average + * where theta \in [0,1] defined as below. + * After the correction, the new solution does not exceed the user-given + * exact global maximum/minimum values. Meanwhile, the new solution's cell average + * equals to the old solution's cell average. + */ + double theta = 1.0; + if (std::abs(max_solution_local-local_solution_average) > std::numeric_limits::min()) + { + theta = std::min(theta, std::abs((max_solution_exact_global-local_solution_average) + / (max_solution_local-local_solution_average))); + } + if (std::abs(min_solution_local-local_solution_average) > std::numeric_limits::min()) + { + theta = std::min(theta, std::abs((min_solution_exact_global-local_solution_average) + / (min_solution_local-local_solution_average))); + } + + /* Modify the advection degrees of freedom of the numerical solution. + * Note that we are using DG elements, so every DoF on a locally owned cell is locally owned; + * this means that we do not need to check whether the 'distributed_solution' vector actually + * stores the element we read from/write to here. + */ + for (unsigned int j = 0; + j < finite_element.base_element(advection_field.base_element(introspection)).dofs_per_cell; + ++j) + { + const unsigned int support_point_index = finite_element.component_to_system_index( + (advection_field.is_temperature() + ? + introspection.component_indices.temperature + : + introspection.component_indices.compositional_fields[advection_field.compositional_variable] + ), + /*dof index within component=*/ j); + const double solution_value = solution(local_dof_indices[support_point_index]); + const double limited_solution_value = theta * (solution_value-local_solution_average) + local_solution_average; + distributed_solution(local_dof_indices[support_point_index]) = limited_solution_value; + } + } + } + + distributed_solution.compress(VectorOperation::insert); + // now get back to the original vector + solution.block(block_idx) = distributed_solution.block(block_idx); + } + + + + namespace internal + { + struct KXRCFCellData + { + double cell_norm; + double inflow_face_jump; + double inflow_face_area; + + KXRCFCellData() + : cell_norm(0.0) + , inflow_face_jump(0.0) + , inflow_face_area(0.0) + {} + }; + } + + + template + template + void + Simulator::compute_KXRCF_indicators(Vector &KXRCF_indicators, + const AdvectionField &advection_field) const + { + // stuff for computing the integrals in cells and cell faces + QGauss quadrature(2); + QGauss face_quadrature(2); + + FEValues fe_values(*mapping, + finite_element, + quadrature, + update_values); + + FEFaceValues fe_face_values(*mapping, + finite_element, + face_quadrature, + update_values | + update_normal_vectors | + update_JxW_values); + + FESubfaceValues fe_subface_values(*mapping, + finite_element, + face_quadrature, + update_values | + update_normal_vectors | + update_JxW_values); + + FEFaceValues neighbor_fe_face_values(*mapping, + finite_element, + face_quadrature, + update_values | + update_JxW_values); + + const FEValuesExtractors::Scalar &field_extractor = advection_field.scalar_extractor(introspection); + + const unsigned int n_q_points = fe_values.n_quadrature_points; + const unsigned int n_face_q_points = fe_face_values.n_quadrature_points; + + std::vector field_values(n_q_points); + std::vector face_field_values(n_face_q_points); + std::vector neighbor_face_field_values(n_face_q_points); + + std::vector> face_velocity_values(n_face_q_points); + std::vector> face_mesh_velocity_values(n_face_q_points); + + // KXRCF indicator requires the computation of (1) cell maximum norm; + // (2) integration of the jump over inflow faces of each field. To + // avoid doing integration twice on each cell face, we set up a map + // to store the cell data and then loop over the cells using the same + // strategy as in Assemblers::AdvectionSystemInteriorFace::execute(). + std::map data; + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + data.insert(std::make_pair(cell->active_cell_index(), + internal::KXRCFCellData())); + + // Now loop over locally owned cells and compute the maximum norm + // and the jump over inflow faces. + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + { + auto cell_data = data.find(cell->active_cell_index()); + + // Compute the maximum norm of field values in this cell. + fe_values.reinit(cell); + fe_values[field_extractor].get_function_values(solution, field_values); + + double cell_norm = 0.0; + for (unsigned int q = 0; q < n_q_points; ++q) + cell_norm = std::max(cell_norm, std::abs(field_values[q])); + + cell_data->second.cell_norm = cell_norm; + + // Compute the jump of field values over inflow faces. + for (const unsigned int face_no : cell->face_indices()) + { + if (!cell->at_boundary(face_no) || cell->has_periodic_neighbor(face_no)) + { + const typename DoFHandler::cell_iterator + neighbor = cell->neighbor_or_periodic_neighbor(face_no); + + const bool cell_has_periodic_neighbor = cell->has_periodic_neighbor(face_no); + + // 'neighbor' defined above is NOT active_cell_iterator, this includes cells + // that are refined + if (!neighbor->has_children()) + { + if (neighbor->level() == cell->level() && + neighbor->is_active() && + (((neighbor->is_locally_owned()) && (cell->index() < neighbor->index())) + || + ((!neighbor->is_locally_owned()) && (cell->subdomain_id() < neighbor->subdomain_id())))) + { + // cell and neighbor are equal-sized, and cell has been chosen to + // assemble this face, so calculate from cell + const unsigned int neighbor2 = + (cell_has_periodic_neighbor + ? + cell->periodic_neighbor_of_periodic_neighbor(face_no) + : + cell->neighbor_of_neighbor(face_no)); + + auto neighbor_data = data.find(neighbor->active_cell_index()); + + // Set up face values. + fe_face_values.reinit(cell, face_no); + fe_face_values[field_extractor].get_function_values(solution, face_field_values); + + fe_face_values[introspection.extractors.velocities].get_function_values( + solution, face_velocity_values); + if (parameters.mesh_deformation_enabled) + fe_face_values[introspection.extractors.velocities].get_function_values( + mesh_deformation->mesh_velocity, face_mesh_velocity_values); + + // Set up neighbor face values. + neighbor_fe_face_values.reinit(neighbor, neighbor2); + neighbor_fe_face_values[field_extractor].get_function_values(solution, neighbor_face_field_values); + + const double face_area = cell->face(face_no)->measure(); + + double v_n = 0.0; + double face_jump = 0.0; + for (unsigned int q = 0; q < n_face_q_points; ++q) + { + Tensor<1,dim> current_u = face_velocity_values[q]; + if (parameters.mesh_deformation_enabled) + current_u -= face_mesh_velocity_values[q]; + + v_n += (current_u * fe_face_values.normal_vector(q)) * fe_face_values.JxW(q); + + face_jump += (face_field_values[q] - neighbor_face_field_values[q]) * fe_face_values.JxW(q); + } + + if (v_n < 0.0) // inflow + { + cell_data->second.inflow_face_jump += face_jump; + cell_data->second.inflow_face_area += face_area; + } + else // outflow + { + if (neighbor_data != data.end()) + { + neighbor_data->second.inflow_face_jump -= face_jump; + neighbor_data->second.inflow_face_area += face_area; + } + } + + } + else + { + // Neighbor is taking responsibility for assembly of this face, because + // either (1) neighbor is coarser, or + // (2) neighbor is equally-sized and + // (a) neighbor is on a different subdomain, with lower subdomain_id(), or + // (b) neighbor is on the same subdomain and has lower index(). + } + } + else + { + // Neighbor has children, so always assemble from here. + const unsigned int neighbor2 = + (cell_has_periodic_neighbor + ? + cell->periodic_neighbor_face_no(face_no) + : + cell->neighbor_face_no(face_no)); + + // Loop over subfaces. We know that the neighbor is finer, so we could loop over the subfaces of the current + // face. but if we are at a periodic boundary, then the face of the current cell has no children, so instead use + // the children of the periodic neighbor's corresponding face since we know that the letter does indeed have + // children (because we know that the neighbor is refined). + typename DoFHandler::face_iterator neighbor_face = neighbor->face(neighbor2); + for (unsigned int subface_no = 0; subface_no < neighbor_face->n_children(); ++subface_no) + { + const typename DoFHandler::active_cell_iterator neighbor_child = + (cell_has_periodic_neighbor + ? + cell->periodic_neighbor_child_on_subface(face_no, subface_no) + : + cell->neighbor_child_on_subface(face_no, subface_no)); + + auto neighbor_data = data.find(neighbor_child->active_cell_index()); + + // Set up subface values. + fe_subface_values.reinit(cell, face_no, subface_no); + fe_subface_values[field_extractor].get_function_values(solution, face_field_values); + + fe_subface_values[introspection.extractors.velocities].get_function_values( + solution, face_velocity_values); + if (parameters.mesh_deformation_enabled) + fe_subface_values[introspection.extractors.velocities].get_function_values( + mesh_deformation->mesh_velocity, face_mesh_velocity_values); + + neighbor_fe_face_values.reinit(neighbor_child, neighbor2); + neighbor_fe_face_values[field_extractor].get_function_values(solution, neighbor_face_field_values); + + const double face_area = neighbor_child->face(neighbor2)->measure(); + + double v_n = 0.0; + double face_jump = 0.0; + for (unsigned int q = 0; q < n_face_q_points; ++q) + { + Tensor<1,dim> current_u = face_velocity_values[q]; + if (parameters.mesh_deformation_enabled) + current_u -= face_mesh_velocity_values[q]; + + v_n += (current_u * fe_subface_values.normal_vector(q)) * fe_subface_values.JxW(q); + + face_jump += (face_field_values[q] - neighbor_face_field_values[q]) * fe_subface_values.JxW(q); + } + + if (v_n < 0.0) // inflow + { + cell_data->second.inflow_face_jump += face_jump; + cell_data->second.inflow_face_area += face_area; + } + else // outflow + { + if (neighbor_data != data.end()) + { + neighbor_data->second.inflow_face_jump -= face_jump; + neighbor_data->second.inflow_face_area += face_area; + } + } + } + } + } + else + { + // The current face is a boundary face. + } + } + } + + // Finally, compute the KXRCF indicator and pick the largest one as + // troubled cell indicator. + const unsigned int degree = advection_field.polynomial_degree(introspection); + const double power = (degree + 1.0) / 2.0; + const double scaling_factor = 1. / global_Omega_diameter; + + const unsigned int block_idx = advection_field.block_index(introspection); + const double epsilon = std::max(1e-20, + 1e-3 * (solution.block(block_idx).max() - solution.block(block_idx).min())); + + KXRCF_indicators.reinit(triangulation.n_active_cells()); + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + { + // The characteristic mesh spacing is calculated by d/D, where + // d and D are the diameters of the present cell and the geometry + // model, respectively. + const double h_pow = std::pow(cell->diameter() * scaling_factor, power); + const auto cell_data = data.find(cell->active_cell_index()); + if (cell_data->second.inflow_face_area == 0.0) + KXRCF_indicators[cell->active_cell_index()] = 0.0; + else + KXRCF_indicators[cell->active_cell_index()] = + std::abs(cell_data->second.inflow_face_jump) / + (h_pow * cell_data->second.inflow_face_area * (cell_data->second.cell_norm + epsilon)); + } + } + + + template + void + Simulator::apply_WENO_limiter_to_dg_solutions(const AdvectionField &advection_field) + { + // Skip if this composition field is being set from the volume_of_fluid handler + if (!advection_field.is_temperature() && + parameters.volume_of_fluid_tracking_enabled) + if (volume_of_fluid_handler->field_index_for_name(introspection.name_for_compositional_index(advection_field.compositional_variable)) + != volume_of_fluid_handler->get_n_fields()) + return; + + const unsigned int degree = advection_field.polynomial_degree(introspection); + AssertThrow(degree < 3, + ExcMessage("Currently, WENO limiter has only been implemented for " + "advection fields with polynomial degree 1 or 2.")); + + // We need two objects of FEValues in the following computation: + // fe_values: computes the cell average and smoothness indicator of the target cell, + // and provides the coordinates and JxWs of the quadrature points of the target + // cell when computing the smoothness indicators of neighbor cells; + // neighbor_fe_values: computes the cell averages and field derivatives of + // neighbor cells at quadrature points of the target cell. With the derivatives + // at hand, we can use the JxWs provided by fe_values to compute the smoothness + // indicators of neighbor cells. + // The coordinates of quadrature points of neighbor_fe_values change from cell to + // cell since the mesh may not be Cartesian, hence neighbor_fe_values has to be + // reconstructed in each cell. + Quadrature quadrature(advection_field.is_temperature() ? + introspection.quadratures.temperature : + introspection.quadratures.compositional_fields); + + UpdateFlags update_flags = update_values | update_gradients | + update_quadrature_points | update_JxW_values | + (degree > 1 ? update_hessians : update_default); + + UpdateFlags neighbor_update_flags = update_values | update_gradients | + (degree > 1 ? update_hessians : update_default); + + FEValues fe_values(*mapping, finite_element, quadrature, update_flags); + + const unsigned int n_q_points = fe_values.n_quadrature_points; + const unsigned int field_dofs_per_cell = finite_element.base_element(advection_field.base_element(introspection)).dofs_per_cell; + const unsigned int field_component = advection_field.component_index(introspection); + + const FEValuesExtractors::Scalar &field_extractor = advection_field.scalar_extractor(introspection); + + std::vector field_values(n_q_points); + std::vector> field_gradients(n_q_points); + std::vector> field_hessians(n_q_points); + + std::vector> neighbor_quadrature_points(n_q_points); + + std::vector reconstructed_values(n_q_points); + + std::vector cell_averages; + std::vector smoothness_indicators; + std::vector linear_weights; + std::vector nonlinear_weights; + std::vector> field_values_in_target_cell; + + // stuff for projecting the reconstructed polynomial onto grid nodes + FullMatrix cell_matrix(field_dofs_per_cell, field_dofs_per_cell); + Vector cell_rhs(field_dofs_per_cell); + Vector cell_solution(field_dofs_per_cell); + std::vector phi_field(field_dofs_per_cell); + std::vector cell_dof_indices(finite_element.dofs_per_cell); + + // Detect troubled cells with KXRCF indicator. + Vector KXRCF_indicators; + compute_KXRCF_indicators(KXRCF_indicators, advection_field); + + TrilinosWrappers::MPI::BlockVector distributed_solution(introspection.index_sets.system_partitioning, + mpi_communicator); + + const unsigned int block_idx = advection_field.block_index(introspection); + distributed_solution.block(block_idx) = solution.block(block_idx); + + const double KXRCF_indicator_threshold = (advection_field.is_temperature() ? + parameters.temperature_KXRCF_indicator_threshold : + parameters.composition_KXRCF_indicator_threshold[advection_field.compositional_variable]); + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + { + // If the KXRCF indicator of the present cell is lower than the threshold, + // then we do not have to do WENO reconstruction for the cell. + if (KXRCF_indicators[cell->active_cell_index()] < KXRCF_indicator_threshold) + continue; + + cell_averages.clear(); + smoothness_indicators.clear(); + field_values_in_target_cell.clear(); + linear_weights.clear(); + + // Compute the field values and derivatives at quadrature points. + fe_values.reinit(cell); + fe_values[field_extractor].get_function_values(solution, field_values); + fe_values[field_extractor].get_function_gradients(solution, field_gradients); + if (degree > 1) + fe_values[field_extractor].get_function_hessians(solution, field_hessians); + + // Compute the cell average and the smoothness indicator. + double cell_measure = 0.0; + double field_integral = 0.0; + double field_gradient_square_integral = 0.0; + double field_hessian_square_integral = 0.0; + for (unsigned int q = 0; q < n_q_points; ++q) + { + cell_measure += fe_values.JxW(q); + field_integral += field_values[q] * fe_values.JxW(q); + + field_gradient_square_integral += field_gradients[q].norm_square() * fe_values.JxW(q); + if (degree > 1) + field_hessian_square_integral += field_hessians[q].norm_square() * fe_values.JxW(q); + } + + double beta = field_gradient_square_integral; + if (degree > 1) + beta += field_hessian_square_integral * cell_measure; + + cell_averages.push_back(field_integral / cell_measure); + smoothness_indicators.push_back(beta); + field_values_in_target_cell.push_back(field_values); + // Linear weight of the target cell cannot be determined at present time + // (because we don't know how many neighbors the target cell has). Fill in + // the position with an invalid value. + linear_weights.push_back(numbers::signaling_nan()); + + // Loop over the neighbor cells. + for (const unsigned int face_no : cell->face_indices()) + { + if (cell->face(face_no)->at_boundary()) + continue; + + const typename DoFHandler::cell_iterator + neighbor = cell->neighbor_or_periodic_neighbor(face_no); + + // 'neighbor' defined above is NOT active_cell_iterator, this includes cells + // that are refined + if (!neighbor->has_children()) + { + // Compute the field values at quadrature points of the neighbor. + + // Compute the field derivatives at quadrature points of the target cell. + // To do so, we need to transform the quadrature points of the target cell + // to the 'unit cell' defined by the neighbor. + mapping->transform_points_real_to_unit_cell(neighbor, + fe_values.get_quadrature_points(), + neighbor_quadrature_points); + + Quadrature neighbor_quadrature(neighbor_quadrature_points, + quadrature.get_weights()); + + FEValues neighbor_fe_values(*mapping, + finite_element, + neighbor_quadrature, + neighbor_update_flags); + + neighbor_fe_values.reinit(neighbor); + neighbor_fe_values[field_extractor].get_function_values(solution, field_values); + neighbor_fe_values[field_extractor].get_function_gradients(solution, field_gradients); + if (degree > 1) + neighbor_fe_values[field_extractor].get_function_hessians(solution, field_hessians); + + // Compute the cell average and the smoothness indicator. + double field_integral = 0.0; + double field_gradient_square_integral = 0.0; + double field_hessian_square_integral = 0.0; + for (unsigned int q = 0; q < n_q_points; ++q) + { + field_integral += field_values[q] * fe_values.JxW(q); + field_gradient_square_integral += field_gradients[q].norm_square() * fe_values.JxW(q); + if (degree > 1) + field_hessian_square_integral += field_hessians[q].norm_square() * fe_values.JxW(q); + } + + double beta = field_gradient_square_integral; + if (degree > 1) + beta += field_hessian_square_integral * cell_measure; + + cell_averages.push_back(field_integral / cell_measure); + smoothness_indicators.push_back(beta); + field_values_in_target_cell.push_back(field_values); + linear_weights.push_back(parameters.WENO_linear_weight); + } + else + { + // Neighbor has children. We use a half of the children that are adjacent to + // the target cell (i.e., children that share a face with the target cell) to + // compute the cell average and smoothness indicator. + const unsigned int neighbor_face_no = + (cell->has_periodic_neighbor(face_no) + ? + cell->periodic_neighbor_face_no(face_no) + : + cell->neighbor_face_no(face_no)); + + typename DoFHandler::face_iterator neighbor_face = neighbor->face(neighbor_face_no); + for (unsigned int subface_no = 0; subface_no < neighbor_face->n_children(); ++subface_no) + { + const unsigned int child_no = + GeometryInfo::child_cell_on_face(RefinementCase::isotropic_refinement, + neighbor_face_no, + subface_no, + neighbor->face_orientation(face_no), + neighbor->face_flip(face_no), + neighbor->face_rotation(face_no)); + + typename DoFHandler::active_cell_iterator + neighbor_child = neighbor->child(child_no); + + // Compute the field derivatives at quadrature points of the target cell. + // To do so, we need to transform the quadrature points of the target cell + // to the 'unit cell' defined by the neighbor. + mapping->transform_points_real_to_unit_cell(neighbor_child, + fe_values.get_quadrature_points(), + neighbor_quadrature_points); + + Quadrature neighbor_quadrature(neighbor_quadrature_points, + quadrature.get_weights()); + + FEValues neighbor_fe_values(*mapping, + finite_element, + neighbor_quadrature, + neighbor_update_flags); + + neighbor_fe_values.reinit(neighbor_child); + neighbor_fe_values[field_extractor].get_function_values(solution, field_values); + neighbor_fe_values[field_extractor].get_function_gradients(solution, field_gradients); + if (degree > 1) + neighbor_fe_values[field_extractor].get_function_hessians(solution, field_hessians); + + // Compute the cell average and the smoothness indicator. + double field_integral = 0.0; + double field_gradient_square_integral = 0.0; + double field_hessian_square_integral = 0.0; + for (unsigned int q = 0; q < n_q_points; ++q) + { + field_integral += field_values[q] * fe_values.JxW(q); + field_gradient_square_integral += field_gradients[q].norm_square() * fe_values.JxW(q); + if (degree > 1) + field_hessian_square_integral += field_hessians[q].norm_square() * fe_values.JxW(q); + } + + double beta = field_gradient_square_integral; + if (degree > 1) + beta += field_hessian_square_integral * cell_measure; + + cell_averages.push_back(field_integral / cell_measure); + smoothness_indicators.push_back(beta); + field_values_in_target_cell.push_back(field_values); + linear_weights.push_back(parameters.WENO_linear_weight / neighbor_face->n_children()); + } + } + } + + // Compute the linear weights of the target cell. + linear_weights[0] = 1.0; + for (unsigned int i = 1; i < linear_weights.size(); ++i) + linear_weights[0] -= linear_weights[i]; + + // Compute the nonlinear weights. + const double epsilon = 1e-6 * std::accumulate(smoothness_indicators.begin(), + smoothness_indicators.end(), + 0.0); + nonlinear_weights.clear(); + for (unsigned int i = 0; i < linear_weights.size(); ++i) + nonlinear_weights.push_back(linear_weights[i] / + Utilities::fixed_power<2,double>(smoothness_indicators[i] + epsilon)); + + // Normalize the nonlinear weights. + const double sum = std::accumulate(nonlinear_weights.begin(), + nonlinear_weights.end(), + 0.0); + + for (unsigned int i = 0; i < nonlinear_weights.size(); ++i) + nonlinear_weights[i] /= sum; + + // Compute the values of the reconstructed polynomial at quadrature points. + std::fill(reconstructed_values.begin(), reconstructed_values.end(), 0.0); + for (unsigned int q = 0; q < n_q_points; ++q) + for (unsigned int i = 0; i < nonlinear_weights.size(); ++i) + reconstructed_values[q] += (field_values_in_target_cell[i][q] - cell_averages[i] + cell_averages[0]) * nonlinear_weights[i]; + + // Project the reconstructed polynomial onto grid nodes. + cell_matrix = 0; + cell_rhs = 0; + for (unsigned int q = 0; q < n_q_points; ++q) + { + for (unsigned int i = 0, i_field = 0; i_field < field_dofs_per_cell; /*increment at end of loop*/) + { + if (finite_element.system_to_component_index(i).first == field_component) + { + phi_field[i_field] = fe_values[field_extractor].value(i, q); + ++i_field; + } + ++i; + } + + for (unsigned int i = 0; i < field_dofs_per_cell; ++i) + { + cell_rhs(i) += phi_field[i] * reconstructed_values[q] * fe_values.JxW(q); + + for (unsigned int j = 0; j < field_dofs_per_cell; ++j) + cell_matrix(i, j) += phi_field[i] * phi_field[j] * fe_values.JxW(q); + } + } + + cell_matrix.gauss_jordan(); + cell_matrix.vmult(cell_solution, cell_rhs); + + // Copy the results to the distributed solution vector. + cell->get_dof_indices(cell_dof_indices); + for (unsigned int i = 0, i_field = 0; i_field < field_dofs_per_cell; /*increment at end of loop*/) + { + if (finite_element.system_to_component_index(i).first == field_component) + { + distributed_solution[cell_dof_indices[i]] = cell_solution(i_field); + ++i_field; + } + ++i; + } + } + + // Update the solution vector. + distributed_solution.compress(VectorOperation::insert); + solution.block(block_idx) = distributed_solution.block(block_idx); + } +} + +// explicit instantiations +namespace aspect +{ + +#define INSTANTIATE(dim) \ + template void Simulator::apply_BP_limiter_to_dg_solutions(const AdvectionField &advection_field); \ + template void Simulator::compute_KXRCF_indicators(Vector &, \ + const AdvectionField &) const; \ + template void Simulator::compute_KXRCF_indicators(Vector &, \ + const AdvectionField &) const; \ + template void Simulator::apply_WENO_limiter_to_dg_solutions(const AdvectionField &advection_field); + + ASPECT_INSTANTIATE(INSTANTIATE) + +#undef INSTANTIATE +} diff --git a/source/simulator/parameters.cc b/source/simulator/parameters.cc index 630301c8fce..3fd32e22c08 100644 --- a/source/simulator/parameters.cc +++ b/source/simulator/parameters.cc @@ -1139,47 +1139,53 @@ namespace aspect "is largely empirically decided -- it must be large enough to ensure " "the bilinear form is coercive, but not so large as to penalize " "discontinuity at all costs."); - prm.declare_entry ("Use limiter for discontinuous temperature solution", "false", - Patterns::Bool (), - "Whether to apply the bound preserving limiter as a correction after computing " - "the discontinuous temperature solution. The limiter will only have an " - "effect if the 'Global temperature maximum' and " - "'Global temperature minimum' parameters are defined in the .prm file. " - "This limiter keeps the discontinuous solution in the range given by " - "'Global temperature maximum' and 'Global temperature minimum'. " - "Because this limiter modifies the solution it no longer " - "satisfies the assembled equation. Therefore, " - "the nonlinear residual for this field is meaningless, and in nonlinear " - "solvers we will ignore the residual for this field to evaluate " + prm.declare_entry ("Limiter for discontinuous temperature solution", "none", + Patterns::Selection(DGLimiterType::pattern()), + "The type of limiter to apply to the discontinuous temperature solution. " + "Available options are:\n" + " * `bound preserving': a limiter that keeps the discontinuous solution " + "in the range given by `Global temperature maximum' and `Global temperature " + "minimum'.\n" + " * `WENO': a limiter that eliminates spurious oscillations by reconstructing " + "a polynomial function in places where shocks are detected, as described in " + "\\cite{zhong:etal:2013}.\n" + " * `none`: if chosen, no limiter is applied to the discontinuous solution.\n" + "Note that if this entry is not set to `none', then the limiter will modify " + "the solution, so the nonlinear residual for this field is meaningless, and " + "in nonlinear solvers we will ignore the residual for this field to evaluate " "if the nonlinear solver has converged."); - prm.declare_entry ("Use limiter for discontinuous composition solution", "false", - Patterns::List(Patterns::Bool()), - "Whether to apply the bound preserving limiter as a correction after having " - "the discontinuous composition solution. The limiter will only have an " - "effect if the 'Global composition maximum' and " - "'Global composition minimum' parameters are defined in the .prm file. " - "This limiter keeps the discontinuous solution in the range given by " - "Global composition maximum' and 'Global composition minimum'. " - "The number of input values in this parameter separated by ',' has to be " - "one or the number of the compositional fields. When only one value " - "is supplied, this same value is assumed for all compositional fields, otherwise " - "each value represents if the limiter should be applied to the respective " - "compositional field. " - "Because this limiter modifies the solution it no longer " - "satisfies the assembled equation. Therefore, " - "the nonlinear residual for this field is meaningless, and in nonlinear " - "solvers we will ignore the residual for this field to evaluate " + prm.declare_entry ("Limiter for discontinuous composition solution", "none", + Patterns::List(Patterns::Selection(DGLimiterType::pattern())), + "The type of limiter to apply to the discontinuous composition solution. The " + "number of the input values separated by ',' has to be one or the same as the " + "number of the compositional fields. When only one value is supplied, this " + "same value is assumed for all compositional fields.\n" + "Available options are:\n" + " * `bound preserving': a limiter that keeps the discontinuous solution " + "in the range given by `Global temperature maximum' and `Global temperature " + "minimum'.\n" + " * `WENO': a limiter that eliminates spurious oscillations by reconstructing " + "a polynomial function in places where shocks are detected, as described in " + "\\cite{zhong:etal:2013}.\n" + " * `none`: if chosen, no limiter is applied to the discontinuous solution.\n" + "Note that if this entry is not set to `none', then the limiter will modify " + "the solution, so the nonlinear residual for this field is meaning less, and " + "in nonlinear solvers we will ignore the residual for this field to evaluate " "if the nonlinear solver has converged."); prm.declare_entry ("Global temperature maximum", boost::lexical_cast(std::numeric_limits::max()), Patterns::Double (), "The maximum global temperature value that will be used in the bound preserving " - "limiter for the discontinuous solutions from temperature advection fields."); + "limiter for the discontinuous solutions from temperature advection fields. " + "The input value is active only when `Limiter for discontinuous temperature " + "solution' is set to `bound preserving'."); prm.declare_entry ("Global temperature minimum", boost::lexical_cast(std::numeric_limits::lowest()), Patterns::Double (), "The minimum global temperature value that will be used in the bound preserving " - "limiter for the discontinuous solutions from temperature advection fields."); + "limiter for the discontinuous solutions from temperature advection fields." + "The input value is active only when `Limiter for discontinuous temperature " + "solution' is set to `bound preserving'."); prm.declare_entry ("Global composition maximum", boost::lexical_cast(std::numeric_limits::max()), Patterns::List(Patterns::Double ()), @@ -1187,7 +1193,9 @@ namespace aspect "limiter for the discontinuous solutions from composition advection fields. " "The number of the input 'Global composition maximum' values separated by ',' has to be " "one or the same as the number of the compositional fields. When only one value " - "is supplied, this same value is assumed for all compositional fields."); + "is supplied, this same value is assumed for all compositional fields. " + "The input value is active only when `Limiter for discontinuous composition " + "solution' is set to `bound preserving'."); prm.declare_entry ("Global composition minimum", boost::lexical_cast(std::numeric_limits::lowest()), Patterns::List(Patterns::Double ()), @@ -1195,7 +1203,32 @@ namespace aspect "limiter for the discontinuous solutions from composition advection fields. " "The number of the input 'Global composition minimum' values separated by ',' has to be " "one or the same as the number of the compositional fields. When only one value " - "is supplied, this same value is assumed for all compositional fields."); + "is supplied, this same value is assumed for all compositional fields.\n" + "The input value is active only when `Limiter for discontinuous composition " + "solution' is set to `bound preserving'."); + prm.declare_entry ("Temperature KXRCF indicator threshold", "1.0", + Patterns::Double(0.0), + "The threshold of KXRCF indicator for the temperature field, as described in " + "\\cite{Krivodonova:etal:2004}. If the KXRCF indicator of a cell is greater than " + "the threshold, then the cell is marked as 'troubled' and will be smoothed by " + "the WENO limiter.\n" + "The input value is active only when `Limiter for discontinuous temperature " + "solution' is set to `WENO'."); + prm.declare_entry ("Composition KXRCF indicator threshold", "1.0", + Patterns::List(Double(0.0)), + "The threshold of KXRCF indicator for the temperature field, as described in " + "\\cite{Krivodonova:etal:2004}. If the KXRCF indicator of a cell is greater than " + "the threshold, then the cell is marked as 'troubled' and will be smoothed by " + "the WENO limiter. The number of the input values separated by ',' has to be one " + "or the same as the number of the compositional fields. When only one value is " + "supplied, this same value is assumed for all compositional fields.\n" + "The input value is active only when `Limiter for discontinuous composition " + "solution' is set to `WENO'."); + prm.declare_entry ("WENO linear weight of neighbor cells", "0.001", + Patterns::Double(0.0, 1.0), + "The linear weight of neighbor cells in WENO scheme. The larger this value is, " + "the more weight is given to information from neighbor cells when reconstructing the polynomial " + "reconstruction. See \\cite{zhong:etal:2013} for detail."); } prm.leave_subsection (); } @@ -1607,15 +1640,15 @@ namespace aspect additional_refinement_time *= year_in_seconds; skip_solvers_on_initial_refinement = prm.get_bool("Skip solvers on initial refinement"); - skip_setup_initial_conditions_on_initial_refinement = prm.get_bool("Skip setup initial conditions on initial refinement"); + skip_setup_initial_temperature_on_initial_refinement = prm.get_bool("Skip setup initial conditions on initial refinement"); - if (skip_setup_initial_conditions_on_initial_refinement == true && skip_solvers_on_initial_refinement == false) + if (skip_setup_initial_temperature_on_initial_refinement == true && skip_solvers_on_initial_refinement == false) AssertThrow(false, ExcMessage("Cannot execute solvers if no initial conditions are set up. " "You must set skip_solvers_on_initial_refinement to true.")); run_postprocessors_on_initial_refinement = prm.get_bool("Run postprocessors on initial refinement"); - if (skip_setup_initial_conditions_on_initial_refinement == true && run_postprocessors_on_initial_refinement == true) + if (skip_setup_initial_temperature_on_initial_refinement == true && run_postprocessors_on_initial_refinement == true) AssertThrow(false, ExcMessage("Cannot run postprocessors if no initial conditions are set up. " "You must set run_postprocessors_on_initial_refinement to false.")); } @@ -1787,13 +1820,16 @@ namespace aspect stabilization_gamma = prm.get_double ("gamma"); discontinuous_penalty = prm.get_double ("Discontinuous penalty"); - use_limiter_for_discontinuous_temperature_solution - = prm.get_bool("Use limiter for discontinuous temperature solution"); - use_limiter_for_discontinuous_composition_solution - = Utilities::possibly_extend_from_1_to_N(Utilities::string_to_bool - (Utilities::split_string_list(prm.get("Use limiter for discontinuous composition solution"))), - n_compositional_fields, - "Use limiter for discontinuous composition solution"); + + limiter_for_discontinuous_temperature_solution = DGLimiterType::parse(prm.get("Limiter for discontinuous temperature solution")); + std::vector x_limiter_for_discontinuous_composition_solution = + Utilities::possibly_extend_from_1_to_N(Utilities::split_string_list(prm.get("Limiter for discontinuous composition solution")), + n_compositional_fields, + "Limiter for discontinuous composition solution"); + limiter_for_discontinuous_composition_solution.resize(n_compositional_fields); + for (unsigned int c = 0; c < n_compositional_fields; ++c) + limiter_for_discontinuous_composition_solution[c] = DGLimiterType::parse(x_limiter_for_discontinuous_composition_solution[c]); + global_temperature_max_preset = prm.get_double ("Global temperature maximum"); global_temperature_min_preset = prm.get_double ("Global temperature minimum"); global_composition_max_preset = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double @@ -1804,6 +1840,14 @@ namespace aspect (Utilities::split_string_list(prm.get ("Global composition minimum"))), n_compositional_fields, "Global composition minimum"); + + temperature_KXRCF_indicator_threshold = prm.get_double ("Temperature KXRCF indicator threshold"); + composition_KXRCF_indicator_threshold = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double + (Utilities::split_string_list(prm.get("Composition KXRCF indicator threshold"))), + n_compositional_fields, + "Composition KXRCF indicator threshold"); + WENO_linear_weight = prm.get_double("WENO linear weight of neighbor cells"); + compositional_fields_with_disabled_boundary_entropy_viscosity = Utilities::split_string_list(prm.get("List of compositional fields with disabled boundary entropy viscosity")); diff --git a/source/simulator/simulator_access.cc b/source/simulator/simulator_access.cc index 990c55a8564..12b4a3affeb 100644 --- a/source/simulator/simulator_access.cc +++ b/source/simulator/simulator_access.cc @@ -321,6 +321,19 @@ namespace aspect + template + void + SimulatorAccess::compute_KXRCF_indicators(Vector &KXRCF_indicators, + const unsigned int field_index) const + { + const typename Simulator::AdvectionField advection_field = + (field_index == 0 ? Simulator::AdvectionField::temperature() + : Simulator::AdvectionField::composition(field_index-1)); + simulator->compute_KXRCF_indicators(KXRCF_indicators, advection_field); + } + + + template const LinearAlgebra::BlockVector & SimulatorAccess::get_current_linearization_point () const diff --git a/source/simulator/solver.cc b/source/simulator/solver.cc index bb9ba0227da..c91409716bb 100644 --- a/source/simulator/solver.cc +++ b/source/simulator/solver.cc @@ -736,20 +736,35 @@ namespace aspect pcout << solver_control.last_step() << " iterations." << std::endl; - if ((advection_field.is_discontinuous(introspection) - && - ( - (advection_field.is_temperature() && parameters.use_limiter_for_discontinuous_temperature_solution) - || - (!advection_field.is_temperature() && parameters.use_limiter_for_discontinuous_composition_solution[advection_field.compositional_variable]) - ))) + if (advection_field.is_discontinuous(introspection)) { - apply_limiter_to_dg_solutions(advection_field); - // by applying the limiter we have modified the solution to no longer - // satisfy the equation. Therefore the residual is meaningless and cannot - // converge to zero in nonlinear iterations. Disable residual computation - // for this field. - return 0.0; + const typename Parameters::DGLimiterType::Kind + limiter = (advection_field.is_temperature() ? + parameters.limiter_for_discontinuous_temperature_solution : + parameters.limiter_for_discontinuous_composition_solution[advection_field.compositional_variable]); + + switch (limiter) + { + case Parameters::DGLimiterType::bound_preserving: + apply_BP_limiter_to_dg_solutions(advection_field); + break; + case Parameters::DGLimiterType::WENO: + apply_WENO_limiter_to_dg_solutions(advection_field); + break; + case Parameters::DGLimiterType::none: + break; + default: + AssertThrow(false, ExcNotImplemented()); + } + + if (limiter != Parameters::DGLimiterType::none) + // by applying the limiter we have modified the solution to no longer + // satisfy the equation. Therefore the residual is meaningless and cannot + // converge to zero in nonlinear iterations. Disable residual computation + // for this field. + return 0.0; + else + return initial_residual; } return initial_residual; diff --git a/source/volume_of_fluid/assembler.cc b/source/volume_of_fluid/assembler.cc index 1d92fb49bfc..a8b734510b2 100644 --- a/source/volume_of_fluid/assembler.cc +++ b/source/volume_of_fluid/assembler.cc @@ -24,7 +24,7 @@ #include #include #include -#include +#include namespace aspect diff --git a/tests/boukare_bulk_composition.prm b/tests/boukare_bulk_composition.prm index f27f707aa8c..1bf3548a74f 100644 --- a/tests/boukare_bulk_composition.prm +++ b/tests/boukare_bulk_composition.prm @@ -13,7 +13,7 @@ set CFL number = 1.0 set Use years in output instead of seconds = true set End time = 6e5 set Output directory = output_bulk_3 -set Nonlinear solver scheme = iterated IMPES +set Nonlinear solver scheme = iterated Advection and Stokes set Max nonlinear iterations = 10 set Nonlinear solver tolerance = 1e-4 set Pressure normalization = surface @@ -114,7 +114,7 @@ end subsection Boundary composition model set Allow fixed composition on outflow boundaries = false set Fixed composition boundary indicators = top - set Model name = initial composition + set List of model names = initial composition end subsection Heating model diff --git a/tests/boukare_conserve_mass.prm b/tests/boukare_conserve_mass.prm index 992bdf8b545..51fe5e3a3de 100644 --- a/tests/boukare_conserve_mass.prm +++ b/tests/boukare_conserve_mass.prm @@ -12,7 +12,7 @@ set Dimension = 2 set CFL number = 1.0 set Use years in output instead of seconds = true set End time = 1e5 -set Nonlinear solver scheme = iterated IMPES +set Nonlinear solver scheme = iterated Advection and Stokes set Max nonlinear iterations = 20 set Nonlinear solver tolerance = 1e-4 set Pressure normalization = surface @@ -87,7 +87,7 @@ end subsection Boundary composition model set Allow fixed composition on outflow boundaries = false set Fixed composition boundary indicators = top - set Model name = initial composition + set List of model names = initial composition end subsection Heating model diff --git a/tests/boukare_latent_heat.prm b/tests/boukare_latent_heat.prm index ce06e480c7f..9d758dca743 100644 --- a/tests/boukare_latent_heat.prm +++ b/tests/boukare_latent_heat.prm @@ -105,7 +105,7 @@ end subsection Boundary composition model set Allow fixed composition on outflow boundaries = false - set Model name = initial composition + set List of model names = initial composition end subsection Heating model diff --git a/tests/box_simple_surface_strain_rate_residual.prm b/tests/box_simple_surface_strain_rate_residual.prm index 058a948ec39..6d44326dd22 100644 --- a/tests/box_simple_surface_strain_rate_residual.prm +++ b/tests/box_simple_surface_strain_rate_residual.prm @@ -23,7 +23,7 @@ end subsection Boundary temperature model set List of model names = initial temperature - set Fixed temperature boundary indicators = 0,1,2,3 # default: + set Fixed temperature boundary indicators = 0,1,2,3 subsection Initial temperature set Maximal temperature = 101.0 @@ -36,7 +36,7 @@ subsection Boundary composition model end subsection Geometry model - set Model name = box # default: + set Model name = box subsection Box set X extent = 100e3 @@ -64,7 +64,7 @@ subsection Boundary velocity model end subsection Initial temperature model - set Model name = function # default: + set Model name = function subsection Function set Function expression = 100.0 diff --git a/tests/box_surface_base_variables.prm b/tests/box_surface_base_variables.prm index ec646ea17f2..dfe51aa1cf9 100644 --- a/tests/box_surface_base_variables.prm +++ b/tests/box_surface_base_variables.prm @@ -53,10 +53,6 @@ subsection Initial temperature model end end -# The parameters below this comment were created by the update script -# as replacement for the old 'Model settings' subsection. They can be -# safely merged with any existing subsections with the same name. - subsection Boundary temperature model set List of model names = function set Fixed temperature boundary indicators = 5 diff --git a/tests/dc_Picard_compressible.prm b/tests/dc_Picard_compressible.prm index d7cbe095d27..49e199d58ad 100644 --- a/tests/dc_Picard_compressible.prm +++ b/tests/dc_Picard_compressible.prm @@ -12,8 +12,6 @@ set Max nonlinear iterations = 150 set Max nonlinear iterations in pre-refinement = 0 set Pressure normalization = surface set Maximum time step = 1e6 - -#set Number of cheap Stokes solver steps = 200 # aspect 2.2 set Adiabatic surface temperature = 1600 set CFL number = 0.5 @@ -23,6 +21,8 @@ subsection Solver parameters set Linear solver tolerance = 1e-2 set Maximum number of expensive Stokes solver steps = 1000 set Number of cheap Stokes solver steps = 200 + + # set Number of cheap Stokes solver steps = 200 # aspect 2.2 end end diff --git a/tests/discontinuous_composition_bound_preserving_limiter.prm b/tests/discontinuous_composition_bound_preserving_limiter.prm index 916efe1e4a8..b6c742644b7 100644 --- a/tests/discontinuous_composition_bound_preserving_limiter.prm +++ b/tests/discontinuous_composition_bound_preserving_limiter.prm @@ -93,8 +93,8 @@ subsection Discretization set Use discontinuous composition discretization = true subsection Stabilization parameters - set Use limiter for discontinuous composition solution = true # apply the limiter to the DG solutions set Global composition maximum = 1.0, 2.0 set Global composition minimum = 0.0, 0.0 + set Limiter for discontinuous temperature solution = bound preserving # apply the limiter to the DG solutions end end diff --git a/tests/discontinuous_composition_bound_preserving_limiter_3d.prm b/tests/discontinuous_composition_bound_preserving_limiter_3d.prm index c381f0011ad..ec38cd80076 100644 --- a/tests/discontinuous_composition_bound_preserving_limiter_3d.prm +++ b/tests/discontinuous_composition_bound_preserving_limiter_3d.prm @@ -94,8 +94,8 @@ subsection Discretization set Use discontinuous composition discretization = true subsection Stabilization parameters - set Use limiter for discontinuous composition solution = true # apply the limiter to the DG solutions set Global composition maximum = 1.0 set Global composition minimum = 0.0 + set Limiter for discontinuous temperature solution = bound preserving # apply the limiter to the DG solutions end end diff --git a/tests/discontinuous_composition_bound_preserving_limiter_selected_fields.prm b/tests/discontinuous_composition_bound_preserving_limiter_selected_fields.prm index 5377bc15ed8..46f333e7aff 100644 --- a/tests/discontinuous_composition_bound_preserving_limiter_selected_fields.prm +++ b/tests/discontinuous_composition_bound_preserving_limiter_selected_fields.prm @@ -80,9 +80,10 @@ subsection Discretization set Use discontinuous composition discretization = true subsection Stabilization parameters - # apply the limiter to the DG solution of field 1, but not field 2 - set Use limiter for discontinuous composition solution = true, false set Global composition maximum = 1.0, 1.0 set Global composition minimum = 0.0, 0.0 + + # apply the limiter to the DG solution of field 1, but not field 2 + set Limiter for discontinuous composition solution = bound preserving, none end end diff --git a/tests/free_surface_iterated_stokes_gmg.prm b/tests/free_surface_iterated_stokes_gmg.prm index b8dd281ac99..26f1f1b6a62 100644 --- a/tests/free_surface_iterated_stokes_gmg.prm +++ b/tests/free_surface_iterated_stokes_gmg.prm @@ -5,7 +5,7 @@ set Dimension = 2 set Output directory = output_free_surface_iterated_stokes_gmg subsection Mesh refinement - set Initial global refinement = 3 # default: 2 + set Initial global refinement = 3 end subsection Solver parameters diff --git a/tests/interpolor_particles_harmonic_average.prm b/tests/interpolor_particles_harmonic_average.prm index a2edd2fa67d..192e4b57515 100644 --- a/tests/interpolor_particles_harmonic_average.prm +++ b/tests/interpolor_particles_harmonic_average.prm @@ -15,7 +15,6 @@ subsection Compositional fields set Mapped particle properties = end - subsection Initial composition model set Model name = unspecified end diff --git a/tests/iterated_advection_and_stokes_DG.prm b/tests/iterated_advection_and_stokes_DG.prm index 656f6126b6b..beca90c7156 100644 --- a/tests/iterated_advection_and_stokes_DG.prm +++ b/tests/iterated_advection_and_stokes_DG.prm @@ -26,8 +26,8 @@ subsection Discretization set Use discontinuous composition discretization = true subsection Stabilization parameters - set Use limiter for discontinuous composition solution = false set beta = 0.5 + set Limiter for discontinuous temperature solution = none end end diff --git a/tests/iterated_advection_and_stokes_DG_limiter.prm b/tests/iterated_advection_and_stokes_DG_limiter.prm index f07acf0e2d0..475859bf86b 100644 --- a/tests/iterated_advection_and_stokes_DG_limiter.prm +++ b/tests/iterated_advection_and_stokes_DG_limiter.prm @@ -28,13 +28,12 @@ subsection Discretization set Use discontinuous composition discretization = true subsection Stabilization parameters - set Use limiter for discontinuous composition solution = true - # porosity will not exceed 0.21 during set model time. # test will have to stay between 0 and 1. set Global composition maximum = 1.0, 1.0 set Global composition minimum = 0.0, 0.0 set beta = 0.5 + set Limiter for discontinuous composition solution = bound preserving end end diff --git a/tests/kinematically_driven_subduction_2d_case1.prm b/tests/kinematically_driven_subduction_2d_case1.prm index 650c7bad14d..fbb9e7219e6 100644 --- a/tests/kinematically_driven_subduction_2d_case1.prm +++ b/tests/kinematically_driven_subduction_2d_case1.prm @@ -154,7 +154,7 @@ subsection Formulation end subsection Postprocess - set List of postprocessors = visualization, velocity statistics, heating statistics, maximum depth of field, composition velocity statistics, viscous dissipation statistics + set List of postprocessors = visualization, velocity statistics, heating statistics subsection Visualization set List of output variables = density, viscosity, strain rate, error indicator diff --git a/tests/kinematically_driven_subduction_2d_case2a.prm b/tests/kinematically_driven_subduction_2d_case2a.prm index bf65782432d..a8e878d1a1e 100644 --- a/tests/kinematically_driven_subduction_2d_case2a.prm +++ b/tests/kinematically_driven_subduction_2d_case2a.prm @@ -155,7 +155,7 @@ subsection Initial temperature model end subsection Postprocess - set List of postprocessors = visualization, velocity statistics, heating statistics, maximum depth of field, composition velocity statistics, viscous dissipation statistics, temperature statistics, trench location, isotherm depth, composition RMS temperature statistics + set List of postprocessors = visualization, velocity statistics, heating statistics, temperature statistics, trench location, isotherm depth, composition RMS temperature statistics subsection Composition velocity statistics # Sum the velocities of the fields that make up the subducting plate. diff --git a/tests/kinematically_driven_subduction_2d_case2b.prm b/tests/kinematically_driven_subduction_2d_case2b.prm index 17af696dd8c..75080d37cfc 100644 --- a/tests/kinematically_driven_subduction_2d_case2b.prm +++ b/tests/kinematically_driven_subduction_2d_case2b.prm @@ -142,7 +142,7 @@ subsection Initial temperature model end subsection Postprocess - set List of postprocessors = visualization, velocity statistics, heating statistics, maximum depth of field, composition velocity statistics, viscous dissipation statistics, temperature statistics, trench location, isotherm depth, composition RMS temperature statistics + set List of postprocessors = visualization, velocity statistics, heating statistics, temperature statistics, trench location, isotherm depth, composition RMS temperature statistics subsection Composition velocity statistics # Sum the velocities of the fields that make up the subducting plate. diff --git a/tests/melt_boukare_EOS.prm b/tests/melt_boukare_EOS.prm index a66742a46dd..fc41a7d9818 100644 --- a/tests/melt_boukare_EOS.prm +++ b/tests/melt_boukare_EOS.prm @@ -89,7 +89,7 @@ end subsection Boundary composition model set Allow fixed composition on outflow boundaries = false set Fixed composition boundary indicators = top - set Model name = initial composition + set List of model names = initial composition end subsection Heating model diff --git a/tests/melt_boukare_averaging.prm b/tests/melt_boukare_averaging.prm index 6f88da030c5..35b70658038 100644 --- a/tests/melt_boukare_averaging.prm +++ b/tests/melt_boukare_averaging.prm @@ -89,7 +89,7 @@ end subsection Boundary composition model set Allow fixed composition on outflow boundaries = false set Fixed composition boundary indicators = top - set Model name = initial composition + set List of model names = initial composition end subsection Heating model diff --git a/tests/melt_visco_plastic_mohr_circle_compaction.prm b/tests/melt_visco_plastic_mohr_circle_compaction.prm index d8ed2f4bc8d..c2b990e5858 100644 --- a/tests/melt_visco_plastic_mohr_circle_compaction.prm +++ b/tests/melt_visco_plastic_mohr_circle_compaction.prm @@ -46,7 +46,6 @@ subsection Prescribed Stokes solution end end - subsection Material model set Model name = EXPERIMENTAL melt visco plastic diff --git a/tests/melt_visco_plastic_mohr_circle_shear.prm b/tests/melt_visco_plastic_mohr_circle_shear.prm index 6a9a756ae1f..12010e10714 100644 --- a/tests/melt_visco_plastic_mohr_circle_shear.prm +++ b/tests/melt_visco_plastic_mohr_circle_shear.prm @@ -45,7 +45,6 @@ subsection Prescribed Stokes solution end end - subsection Material model set Model name = EXPERIMENTAL melt visco plastic diff --git a/tests/no_dirichlet_on_outflow_refine.prm b/tests/no_dirichlet_on_outflow_refine.prm index 74cf957cc47..679c6360e4e 100644 --- a/tests/no_dirichlet_on_outflow_refine.prm +++ b/tests/no_dirichlet_on_outflow_refine.prm @@ -40,7 +40,7 @@ end subsection Boundary temperature model set Fixed temperature boundary indicators = bottom, top set Allow fixed temperature on outflow boundaries = false - set Model name = function + set List of model names = function subsection Function set Variable names = x,y diff --git a/tests/original_prm/original.prm b/tests/original_prm/original.prm index 701798a8f69..ff3555b1b19 100644 --- a/tests/original_prm/original.prm +++ b/tests/original_prm/original.prm @@ -10,4 +10,3 @@ include $ASPECT_SOURCE_DIR/tests/no_flow.prm # directory below this line: set Output directory = output-original_prm - diff --git a/tests/periodic_box_gmg.prm b/tests/periodic_box_gmg.prm index 1cc7f1860ed..941cda5467f 100644 --- a/tests/periodic_box_gmg.prm +++ b/tests/periodic_box_gmg.prm @@ -100,7 +100,7 @@ end subsection Mesh refinement set Additional refinement times = set Initial adaptive refinement = 2 - set Initial global refinement = 5 # default: 2 + set Initial global refinement = 5 set Minimum refinement level = 5 set Refinement fraction = 0.3 set Coarsening fraction = 0.03 diff --git a/tests/prmbackslash.prm b/tests/prmbackslash.prm index 63063c53811..735593109e8 100644 --- a/tests/prmbackslash.prm +++ b/tests/prmbackslash.prm @@ -51,11 +51,6 @@ subsection Mesh refinement set Initial global refinement = 5 end -# The parameters below this comment were created by the update script -# as replacement for the old 'Model settings' subsection. They can be -# safely merged with any existing subsections with the same name. - - subsection Boundary velocity model set Tangential velocity boundary indicators = 1 set Zero velocity boundary indicators = 0, 2, 3 diff --git a/tests/prmbackslash_2.prm b/tests/prmbackslash_2.prm index 7a211a97e07..8495edefa4e 100644 --- a/tests/prmbackslash_2.prm +++ b/tests/prmbackslash_2.prm @@ -1,156 +1,96 @@ -# Based on the sol_c\ -x_2.prm test input, split\ - every input file after 20 -# characters with a \ -backslash. In some cases,\ - split it more than once. -# The output should \ -still be the same. +# Based on the sol_cx_2.prm test input, split every input file after 20 +# characters with a backslash. In some cases, split it more than once. +# The output should still be the same. set Dimension = 2 - -set CFL number \ - = \ +set CFL number = \ 1.0 -set End time \ - = \ +set End time = \ 0 - - -set Resume computati\ -on = \ +set Resume computation = \ false -set Start time \ - = \ +set Start time = \ 0 -set Adiabatic surfac\ -e temperature = \ +set Adiabatic surface temperature = \ 0 -set Surface pressure\ - = \ +set Surface pressure = \ 0 -set Use years in out\ -put instead of seconds = \ +set Use years in output instead of seconds = \ false -set Nonlinear solver \ - scheme = \ +set Nonlinear solver scheme = \ no Advection, iterated Stokes -subsection Boundary \ -temperature model +subsection Boundary temperature model set Model name = b\ ox + set Fixed temperature boundary indicators = 0, 1 end +subsection Discretization + set Stokes velocity polynomial degree = 2 + set Temperature polynomial degree = 2 + set Use locally conservative discretization = false -subsection Discretiz\ -ation - set Stokes velocit\ -y polynomial degree \ - = 2 - set Temperature po\ -lynomial degree \ - = 2 - set Use locally co\ -nservative discretization\ - = false - subsection Stabili\ -zation parameters + subsection Stabilization parameters set alpha = 2 set beta = 0.07\ 8 set cR = 0.5 \ - end -end - + end -subsection Geometry \ -model - set Model name = b\ + subsection Geometry model + set Model name = b\ ox - subsection Box - set X extent = 1\ - + subsection Box + set X extent = 1\ set Y extent = 1\ - set Z extent = 1\ - end -end - + end -subsection Gravity m\ -odel - set Model name = v\ + subsection Gravity model + set Model name = v\ ertical -end + end - -subsection Initial t\ -emperature model - set Model name = p\ + subsection Initial temperature model + set Model name = p\ erturbed box -end + end - -subsection Material \ -model - set Model name = S\ + subsection Material model + set Model name = S\ olCxMaterial -end - - -subsection Mesh refi\ -nement - set Initial adapti\ -ve refinement = 0 \ - - set Initial global \ - refinement = 2 \ + end - - set Strategy \ - = de\ + subsection Mesh refinement + set Initial adaptive refinement = 0 \ + set Initial global refinement = 2 \ + set Strategy = de\ nsity, temperature -end - + end -subsection Boundary temp\ -erature model - set Fixed temperat\ -ure boundary indicators \ - = 0, 1 -end -subsection Boundary velo\ -city model - set Tangential vel\ -ocity boundary indicators\ - = 0,1,2,3 -end + subsection Boundary velocity model + set Tangential velocity boundary indicators = 0,1,2,3 + end + subsection Postprocess + set List of postprocessors = SolCxPostprocessor -subsection Postproce\ -ss - set List of postpr\ -ocessors = SolCxPostprocessor - subsection Depth a\ -verage - set Time between \ - graphical output = 1e8 - end + subsection Depth average + set Time between graphical output = 1e8 + end - subsection Visuali\ -zation - set Number of gr\ -ouped files = 0 - set Output forma\ -t = gnupl\ + subsection Visualization + set Number of grouped files = 0 + set Output format = gnupl\ ot - set Time between \ - graphical output = 0 #\ + set Time between graphical output = 0 #\ default: 1e8 + end + end end end diff --git a/tests/projected_density_vertical_pipe.prm b/tests/projected_density_vertical_pipe.prm index 608164a1ba8..c3ef1d902ac 100644 --- a/tests/projected_density_vertical_pipe.prm +++ b/tests/projected_density_vertical_pipe.prm @@ -13,7 +13,7 @@ set Nonlinear solver tolerance = 1e-5 subsection Boundary composition model set Allow fixed composition on outflow boundaries = true set Fixed composition boundary indicators = top, bottom - set Model name = box + set List of model names = box # The test is to make sure these values are ignored subsection Box diff --git a/tests/shearbox_particle_strain_rate.prm b/tests/shearbox_particle_strain_rate.prm index 725fb4c788a..8e29130e79f 100644 --- a/tests/shearbox_particle_strain_rate.prm +++ b/tests/shearbox_particle_strain_rate.prm @@ -110,7 +110,7 @@ subsection Gravity model end subsection Material model - set Model name = simple # default: + set Model name = simple subsection Simple model set Reference density = 1 diff --git a/tests/shell_surface_base_variables.prm b/tests/shell_surface_base_variables.prm index d47555aad0f..9d65bed9446 100644 --- a/tests/shell_surface_base_variables.prm +++ b/tests/shell_surface_base_variables.prm @@ -8,13 +8,9 @@ set End time = 1e300 set Start time = 0 set Adiabatic surface temperature = 1613.0 set Surface pressure = 0 -set Use years in output instead of seconds = false # default: true +set Use years in output instead of seconds = false set Nonlinear solver scheme = single Advection, single Stokes -# The parameters below this comment were created by the update script -# as replacement for the old 'Model settings' subsection. They can be -# safely merged with any existing subsections with the same name. - subsection Boundary temperature model set List of model names = initial temperature set Fixed temperature boundary indicators = 0, 1 @@ -49,10 +45,10 @@ subsection Material model subsection Simple model set Reference density = 3300 set Reference specific heat = 1250 - set Reference temperature = 1613 # default: 293 - set Thermal conductivity = 1e-6 # default: 4.7 + set Reference temperature = 1613 + set Thermal conductivity = 1e-6 set Thermal expansion coefficient = 2e-5 - set Viscosity = 1e22 # default: 5e24 + set Viscosity = 1e22 end end diff --git a/tests/slab_tip_depth.prm b/tests/slab_tip_depth.prm index 0e1532e6b75..35558060eba 100644 --- a/tests/slab_tip_depth.prm +++ b/tests/slab_tip_depth.prm @@ -5,7 +5,7 @@ # With increasing resolution, it should return the following # maximum depths of the composition fields: # 82, 110 and 82 km, resp. -# The postprocessors viscous dissipation statistics and +# The postprocessors heating statistics and # viscous dissipation per compositional field return # the viscous dissipation over the whole domain and # per area of each field, resp. With enough resolution, @@ -146,7 +146,7 @@ subsection Formulation end subsection Postprocess - set List of postprocessors = velocity statistics, maximum depth of field, viscous dissipation statistics + set List of postprocessors = velocity statistics, maximum depth of field, heating statistics end subsection Solver parameters diff --git a/tests/stokes_solver_fail.prm b/tests/stokes_solver_fail.prm index ad8853c1ca5..0f46b4734dd 100644 --- a/tests/stokes_solver_fail.prm +++ b/tests/stokes_solver_fail.prm @@ -13,5 +13,5 @@ subsection Solver parameters set Maximum number of expensive Stokes solver steps = 0 set Linear solver tolerance = 1e-35 set Stokes solver type = block AMG - end + end end diff --git a/tests/visco_elastic_top_beam.prm b/tests/visco_elastic_top_beam.prm index e27fc9f5fe3..7833dcac8b0 100644 --- a/tests/visco_elastic_top_beam.prm +++ b/tests/visco_elastic_top_beam.prm @@ -43,9 +43,9 @@ subsection Discretization set Use discontinuous composition discretization = true subsection Stabilization parameters - set Use limiter for discontinuous composition solution = true set Global composition maximum = 1.e11, 1.e11, 1.e11, 1.0 set Global composition minimum = -1.e11, -1.e11, -1.e11, 0.0 + set Limiter for discontinuous composition solution = bound preserving end end diff --git a/tests/viscous_top_beam.prm b/tests/viscous_top_beam.prm index 63b291dda6f..e4a1185cb57 100644 --- a/tests/viscous_top_beam.prm +++ b/tests/viscous_top_beam.prm @@ -51,9 +51,9 @@ subsection Discretization set Use discontinuous composition discretization = true subsection Stabilization parameters - set Use limiter for discontinuous composition solution = true set Global composition maximum = 1.e11, 1.e11, 1.e11, 1.0 set Global composition minimum = -1.e11, -1.e11, -1.e11, 0.0 + set Limiter for discontinuous composition solution = bound preserving end end