implement WENO limiter for DG method

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

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

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

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
 

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
 

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
 

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
 

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

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.

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
 

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
 

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

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
 

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
 

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
 

benchmarks/solubility/plugin/solubility.cc

@@ -359,7 +359,7 @@ namespace aspect
     melt_fractions (const MaterialModel::MaterialModelInputs<dim> &in,
                     std::vector<double> &melt_fractions) const
     {
-      for (unsigned int q=0; q<in.temperature.size(); ++q)
+      for (unsigned int q=0; q<in.n_evaluation_points(); ++q)
         {
           const unsigned int water_idx = this->introspection().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<ReactionRateOutputs<dim>>() == 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<MaterialModel::ReactionRateOutputs<dim>> (n_points, this->n_compositional_fields()));
         }

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

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
 

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

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