Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

implement WENO limiter for DG method #5929

Open
wants to merge 22 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from 13 commits
Commits
Show all changes
22 commits
Select commit Hold shift + click to select a range
ce09083
implement WENO limiter for DG method
Jun 21, 2024
b6aaaca
correct spelling mistakes
Jun 22, 2024
68e633c
add necessary headers for compilation without unity build
Jul 4, 2024
4bfb4a8
reformat the .prm file for limiter settings
Jul 4, 2024
13ed3c4
Update source/simulator/parameters.cc
YiminJin Jul 4, 2024
95edb8a
change boundary_preserving to bound_preserving
Jul 4, 2024
b0c4e8f
replace 'boundary preserving' by 'bound preserving'
Jul 5, 2024
a53c298
add citations for WENO limiter and KXRCF indicator
Jul 8, 2024
9d46a62
scale the characteristic mesh spacing by the diameter of Omega
Jul 8, 2024
06e7cb8
add viscoelastic_beam_modified benchmark with WENO limiter
Jul 8, 2024
36ff891
add viscoelastic_bending_beam benchmark with WENO limiter
Jul 8, 2024
0a25a71
add rotate_shape benchmark with WENO limiter
Jul 8, 2024
7103a42
add rotate_shape benchmark with bound-preserving limiter
Jul 8, 2024
59f361b
replace Patterns::Double() by Patterns::List(Patterns::Double()) for …
Jul 11, 2024
d6a1a0e
supplement comments for KXRCF indicator and WENO limiter
Jul 11, 2024
92b5d90
supplement the documentation for KXRCF indicator
Jul 11, 2024
e305776
Update source/simulator/parameters.cc
YiminJin Jul 11, 2024
f5cbc1c
Update source/simulator/parameters.cc
YiminJin Jul 11, 2024
4667dbc
reformat the .prm file for limiter settings
Jul 11, 2024
aea41a4
remove the .prm.bak files
Jul 27, 2024
fa78f1e
fix bug: missing case: comments behind parameter value
Jul 27, 2024
3c8173e
correct the parameter entries related with DG limiter
Jul 27, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
21 changes: 21 additions & 0 deletions benchmarks/advection/rotate_shape_dg_bp.prm
Original file line number Diff line number Diff line change
@@ -0,0 +1,21 @@
# 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
26 changes: 26 additions & 0 deletions benchmarks/advection/rotate_shape_dg_weno.prm
Original file line number Diff line number Diff line change
@@ -0,0 +1,26 @@
# 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
33 changes: 33 additions & 0 deletions benchmarks/viscoelastic_beam_modified/20km_opentopbot_weno.prm
Original file line number Diff line number Diff line change
@@ -0,0 +1,33 @@
# 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
Original file line number Diff line number Diff line change
@@ -0,0 +1,32 @@
# 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
Original file line number Diff line number Diff line change
@@ -0,0 +1,75 @@
#!/usr/bin/env python3
# -*- coding: utf-8 -*-

""" This script changes 'Use limiter for discontinuous temperature/composition solution'
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Have you tested this script? You will likely need to apply it anyway to all files that currently use the BP limiter in benchmarks/ cookbooks/ and tests/. An easy way to apply the script is to run the file contrib/utilities/update_scripts/update_all_files.sh.

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:

# Delete the old entry
use_limiter = subsection[old_entries[i]]["value"]
del subsection[old_entries[i]]

# If limiter is used, then set the value of the new entry to "bound preserving"
if use_limiter == "true":
subsection[new_entries[i]] = {"comment" : "",
"value" : "bound preserving",
"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))
28 changes: 28 additions & 0 deletions doc/manual/citing_aspect.bib
Original file line number Diff line number Diff line change
Expand Up @@ -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}
}
42 changes: 40 additions & 2 deletions include/aspect/parameters.h
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -656,12 +691,15 @@ namespace aspect
std::vector<double> stabilization_beta;
double stabilization_gamma;
double discontinuous_penalty;
bool use_limiter_for_discontinuous_temperature_solution;
std::vector<bool> use_limiter_for_discontinuous_composition_solution;
typename DGLimiterType::Kind limiter_for_discontinuous_temperature_solution;
std::vector<typename DGLimiterType::Kind> limiter_for_discontinuous_composition_solution;
double global_temperature_max_preset;
double global_temperature_min_preset;
std::vector<double> global_composition_max_preset;
std::vector<double> global_composition_min_preset;
double temperature_KXRCF_indicator_threshold;
std::vector<double> composition_KXRCF_indicator_threshold;
double WENO_linear_weight;

std::vector<std::string> compositional_fields_with_disabled_boundary_entropy_viscosity;
/**
Expand Down
80 changes: 80 additions & 0 deletions include/aspect/postprocess/visualization/kxrcf_indicator.h
Original file line number Diff line number Diff line change
@@ -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
<http://www.gnu.org/licenses/>.
*/


#ifndef _aspect_postprocess_visualization_kxrcf_indicator_h
#define _aspect_postprocess_visualization_kxrcf_indicator_h

#include <aspect/postprocess/visualization.h>
#include <aspect/simulator_access.h>


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.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please add an explanation what the KXRCF indicator is. What does it stand for and what does it mean to have a high KXRCF in a cell?

*/
template <int dim>
class KXRCFIndicator : public CellDataVectorCreator<dim>,
public SimulatorAccess<dim>
{
public:
/**
* Constructor.
*/
KXRCFIndicator();

/**
* @copydoc CellDataVectorCreator<dim>::execute()
*/
std::pair<std::string, std::unique_ptr<Vector<float>>>
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
23 changes: 21 additions & 2 deletions include/aspect/simulator.h
Original file line number Diff line number Diff line change
Expand Up @@ -1373,10 +1373,29 @@ namespace aspect
* the discontinuous Galerkin solution will stay in the prescribed bounds.
*
* This function is implemented in
* <code>source/simulator/helper_functions.cc</code>.
* <code>source/simulator/limiters.cc</code>.
*/
void apply_BP_limiter_to_dg_solutions (const AdvectionField &advection_field);

/**
* Apply the WENO limiter to the discontinuous Galerkin solutions.
* The WENO scheme implemented in this function is proposed by Zhong and Shu, 2013.
*
* This function is implemented in
* <code>source/simulator/limiters.cc</code>.
*/
void apply_limiter_to_dg_solutions (const AdvectionField &advection_field);
void apply_WENO_limiter_to_dg_solutions (const AdvectionField &advection_field);

/**
* Compute the KXRCF indicator that detect shocks in the advection field.
* The KXRCF indicator is used by the WENO limiter.
*
* This function is implemented in
* <code>source/simulator/limiters.cc</code>.
*/
template <typename T>
void compute_KXRCF_indicators(Vector<T> &KXRCF_indicators,
const AdvectionField &advection_field) const;

/**
* Compute the reactions in case of operator splitting:
Expand Down
Loading