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 1 commit
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
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
{
boundary_preserving,
WENO,
none
};

static
Kind
parse(const std::string &input)
{
if (input == "boundary preserving")
Copy link
Member

Choose a reason for hiding this comment

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

I think it would need to be bound preserving. boundary preserving refers to a geometrical boundary while bound preserving mean preserving the upper and lower bounds (values) of a solution.

return boundary_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 "boundary 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
8 changes: 8 additions & 0 deletions include/aspect/simulator_access.h
Original file line number Diff line number Diff line change
Expand Up @@ -427,6 +427,14 @@ namespace aspect
void
get_artificial_viscosity_composition(Vector<float> &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<float> &KXRCF_indicators,
const unsigned int field_index) const;
/** @} */


Expand Down
123 changes: 123 additions & 0 deletions source/postprocess/visualization/kxrcf_indicator.cc
Original file line number Diff line number Diff line change
@@ -0,0 +1,123 @@
/*
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/>.
*/


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

namespace aspect
{
namespace Postprocess
{
namespace VisualizationPostprocessors
{
template <int dim>
KXRCFIndicator<dim>::KXRCFIndicator()
:
CellDataVectorCreator<dim>("")
{}


template <int dim>
std::pair<std::string, std::unique_ptr<Vector<float>>>
KXRCFIndicator<dim>::execute() const
{
std::pair<std::string, std::unique_ptr<Vector<float>>>
return_value ("KXRCF_indicator",
std::make_unique<Vector<float>>(this->get_triangulation().n_active_cells()));
this->compute_KXRCF_indicators(*return_value.second, field_index);
return return_value;
}


template <int dim>
void
KXRCFIndicator<dim>::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 <int dim>
void
KXRCFIndicator<dim>::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 on each "
"cell."
Copy link
Member

Choose a reason for hiding this comment

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

As above, please extend this documentation by what the KXRCF indicator means and under what circumstances it makes sense to look at it.

"\n\n"
"Physical units: none.")
}
}
}
12 changes: 7 additions & 5 deletions source/simulator/core.cc
Original file line number Diff line number Diff line change
Expand Up @@ -485,14 +485,16 @@ namespace aspect
if (MappingQCache<dim> *map = dynamic_cast<MappingQCache<dim>*>(&(*mapping)))
map->initialize(MappingQGeneric<dim>(4), triangulation);

bool dg_limiter_enabled = parameters.use_limiter_for_discontinuous_temperature_solution;
bool bp_limiter_enabled = (parameters.limiter_for_discontinuous_temperature_solution
== Parameters<dim>::DGLimiterType::boundary_preserving);
for (unsigned int c=0; c<introspection.n_compositional_fields; ++c)
dg_limiter_enabled = dg_limiter_enabled || parameters.use_limiter_for_discontinuous_composition_solution[c];
bp_limiter_enabled = bp_limiter_enabled || (parameters.limiter_for_discontinuous_composition_solution[c]
== Parameters<dim>::DGLimiterType::boundary_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 boundary 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."));

Expand Down
Loading
Loading