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.

Sign up for GitHub

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

[WIP] use 5th-order hybrid PPM/WENO reconstruction #47

Draft
wants to merge 16 commits into
base: development
Choose a base branch
from
Draft

[WIP] use 5th-order hybrid PPM/WENO reconstruction #47

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
16 commits
Select commit Hold shift + click to select a range
d8a1dd7
implement xsPPM7
BenWibking May 14, 2022
91f768e
use 5th-order stencil
BenWibking May 15, 2022
f8c4db3
adjust HydroContact tolerance
BenWibking May 15, 2022
12f02e0
add 4th-order conservative to primitive
BenWibking May 25, 2022
e6292ab
Merge branch 'development' into hybrid-ppm-weno
BenWibking May 26, 2022
aca0535
GPU fix
BenWibking May 28, 2022
f96e96f
mark quokka::ArrayView ctors as __host__ __device__
BenWibking May 28, 2022
2e91f4b
Merge branch 'development' into hybrid-ppm-weno
BenWibking May 28, 2022
8cbbe0d
remove high-order cons2prim
BenWibking May 29, 2022
cf388a1
adjust sedov test params
BenWibking May 29, 2022
8b1b17c
Merge branch 'hybrid-ppm-weno' of github.com:BenWibking/TwoMomentRad …
BenWibking May 29, 2022
cb3f681
update amrex submodule
BenWibking May 29, 2022
4ccf97e
disable HighMach test
BenWibking May 29, 2022
ed6421d
remove HighMach test
BenWibking May 29, 2022
df7d4a2
add ipynb for Sedov test
BenWibking Jun 7, 2022
71c74bf
update KH test params
BenWibking Jun 7, 2022
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
2 changes: 1 addition & 1 deletion extern/amrex
Open in desktop
Submodule amrex updated 169 files
120 changes: 120 additions & 0 deletions extern/sedov.ipynb
View file
Open in desktop

Large diffs are not rendered by default.

213 changes: 213 additions & 0 deletions extern/weno_weights.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,213 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": 66,
"id": "abcb19eb",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from sympy import *"
]
},
{
"cell_type": "code",
"execution_count": 67,
"id": "23e217cb",
"metadata": {},
"outputs": [],
"source": [
"x, qbar, qx, qxx = symbols(\"x qbar q_x q_xx\")"
]
},
{
"cell_type": "code",
"execution_count": 68,
"id": "e5ca877e",
"metadata": {},
"outputs": [],
"source": [
"profile = qbar + x * qx + (x**2 - Rational(1,12)) * qxx"
]
},
{
"cell_type": "code",
"execution_count": 69,
"id": "9cec7e13",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\bar{q}$"
],
"text/plain": [
"qbar"
]
},
"execution_count": 69,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# compute mean value\n",
"integrate(profile, (x, Rational(-1,2), Rational(1,2)))"
]
},
{
"cell_type": "code",
"execution_count": 70,
"id": "2123c1f5",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle - \\frac{q_{xx}}{12} + \\bar{q}$"
],
"text/plain": [
"-q_xx/12 + qbar"
]
},
"execution_count": 70,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# compute point value\n",
"profile.subs(x, 0)"
]
},
{
"cell_type": "code",
"execution_count": 73,
"id": "6e5ca46e",
"metadata": {},
"outputs": [],
"source": [
"from sympy.solvers.solveset import linsolve"
]
},
{
"cell_type": "code",
"execution_count": 92,
"id": "e62a6641",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left( \\frac{3 \\hat{q}}{2} - 2 \\hat{q}_{-1} + \\frac{\\hat{q}_{-2}}{2}, \\ \\frac{\\hat{q}}{2} - \\hat{q}_{-1} + \\frac{\\hat{q}_{-2}}{2}\\right)$"
],
"text/plain": [
"(3*\\hat{q}/2 - 2*\\hat{q}_{-1} + \\hat{q}_{-2}/2, \\hat{q}/2 - \\hat{q}_{-1} + \\hat{q}_{-2}/2)"
]
},
"execution_count": 92,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# compute L stencil from point values\n",
"qhat_m2, qhat_m1, qhat = symbols(\"\\hat{q}_{-2} \\hat{q}_{-1} \\hat{q}\")\n",
"sol_Lstencil = linsolve([profile.subs(x, 0) - qhat, profile.subs(x, -1) - qhat_m1, profile.subs(x, -2) - qhat_m2], [qbar, qx, qxx])\n",
"sol_Lstencil.args[0][1:] # use only qx, qxx here"
]
},
{
"cell_type": "code",
"execution_count": 91,
"id": "e0b77c38",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left( \\frac{\\hat{q}_{+1}}{2} - \\frac{\\hat{q}_{-1}}{2}, \\ - \\hat{q} + \\frac{\\hat{q}_{+1}}{2} + \\frac{\\hat{q}_{-1}}{2}\\right)$"
],
"text/plain": [
"(\\hat{q}_{+1}/2 - \\hat{q}_{-1}/2, -\\hat{q} + \\hat{q}_{+1}/2 + \\hat{q}_{-1}/2)"
]
},
"execution_count": 91,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# compute C stencil from point values\n",
"qhat_m1, qhat, qhat_p1 = symbols(\"\\hat{q}_{-1} \\hat{q} \\hat{q}_{+1}\")\n",
"sol_Cstencil = linsolve([profile.subs(x, 0) - qhat, profile.subs(x, -1) - qhat_m1, profile.subs(x, 1) - qhat_p1], [qbar, qx, qxx])\n",
"sol_Cstencil.args[0][1:] # use only qx, qxx here"
]
},
{
"cell_type": "code",
"execution_count": 90,
"id": "e11fb5a2",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left( - \\frac{3 \\hat{q}}{2} + 2 \\hat{q}_{+1} - \\frac{\\hat{q}_{+2}}{2}, \\ \\frac{\\hat{q}}{2} - \\hat{q}_{+1} + \\frac{\\hat{q}_{+2}}{2}\\right)$"
],
"text/plain": [
"(-3*\\hat{q}/2 + 2*\\hat{q}_{+1} - \\hat{q}_{+2}/2, \\hat{q}/2 - \\hat{q}_{+1} + \\hat{q}_{+2}/2)"
]
},
"execution_count": 90,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# compute R stencil from point values\n",
"qhat, qhat_p1, qhat_p2 = symbols(\"\\hat{q} \\hat{q}_{+1} \\hat{q}_{+2}\")\n",
"sol_Rstencil = linsolve([profile.subs(x, 0) - qhat, profile.subs(x, 1) - qhat_p1, profile.subs(x, 2) - qhat_p2], [qbar, qx, qxx])\n",
"sol_Rstencil.args[0][1:] # use only qx, qxx here"
]
},
{
"cell_type": "markdown",
"id": "5d0ce5d7",
"metadata": {},
"source": [
"We see that the expressions for the moments computed from the cell-centered data are formally *identical* to those computed with the cell-average data!"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "75b469f7",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.4"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
15 changes: 8 additions & 7 deletions src/ArrayUtil.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,16 +8,17 @@
/// \file ArrayUtil.hpp
/// \brief Implements functions to manipulate arrays (CPU only).

#include "AMReX_Array4.H"
#include "AMReX_REAL.H"
#include <vector>

template <typename T>
auto strided_vector_from(std::vector<T> &v, int stride) -> std::vector<T>
{
std::vector<T> strided_v;
for(std::size_t i = 0; i < v.size(); i += stride) {
strided_v.push_back(v[i]);
}
return strided_v; // move semantics implied
auto strided_vector_from(std::vector<T> &v, int stride) -> std::vector<T> {
std::vector<T> strided_v;
for (std::size_t i = 0; i < v.size(); i += stride) {
strided_v.push_back(v[i]);
}
return strided_v; // move semantics implied
}

#endif // ARRAYUTIL_HPP_
10 changes: 5 additions & 5 deletions src/ArrayView_2d.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,7 @@ template <class T, FluxDir N, class Enable = void> struct Array4View {
amrex::Array4<T> arr_;
constexpr static FluxDir indexOrder = N;

explicit Array4View(amrex::Array4<T> arr) : arr_(arr) {}
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE explicit Array4View(amrex::Array4<T> arr) : arr_(arr) {}
};

// X1-flux
Expand All @@ -46,7 +46,7 @@ template <class T> struct Array4View<T, FluxDir::X1, std::enable_if_t<!std::is_c
amrex::Array4<T> arr_;
constexpr static FluxDir indexOrder = FluxDir::X1;

explicit Array4View(amrex::Array4<T> arr) : arr_(arr) {}
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE explicit Array4View(amrex::Array4<T> arr) : arr_(arr) {}

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto operator()(int i, int j, int k,
int n) const noexcept -> T &
Expand All @@ -66,7 +66,7 @@ template <class T> struct Array4View<T, FluxDir::X1, std::enable_if_t<std::is_co
amrex::Array4<T> arr_;
constexpr static FluxDir indexOrder = FluxDir::X1;

explicit Array4View(amrex::Array4<T> arr) : arr_(arr) {}
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE explicit Array4View(amrex::Array4<T> arr) : arr_(arr) {}

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto operator()(int i, int j, int k,
int n) const noexcept -> T
Expand All @@ -88,7 +88,7 @@ template <class T> struct Array4View<T, FluxDir::X2, std::enable_if_t<!std::is_c
amrex::Array4<T> arr_;
constexpr static FluxDir indexOrder = FluxDir::X2;

explicit Array4View(amrex::Array4<T> arr) : arr_(arr) {}
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE explicit Array4View(amrex::Array4<T> arr) : arr_(arr) {}

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto operator()(int i, int j, int k,
int n) const noexcept -> T &
Expand All @@ -108,7 +108,7 @@ template <class T> struct Array4View<T, FluxDir::X2, std::enable_if_t<std::is_co
amrex::Array4<T> arr_;
constexpr static FluxDir indexOrder = FluxDir::X2;

explicit Array4View(amrex::Array4<T> arr) : arr_(arr) {}
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE explicit Array4View(amrex::Array4<T> arr) : arr_(arr) {}

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto operator()(int i, int j, int k,
int n) const noexcept -> T
Expand Down
14 changes: 7 additions & 7 deletions src/ArrayView_3d.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,7 @@ template <class T, FluxDir N, class Enable = void> struct Array4View {
amrex::Array4<T> arr_;
constexpr static FluxDir indexOrder = N;

explicit Array4View(amrex::Array4<T> arr) : arr_(arr) {}
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE explicit Array4View(amrex::Array4<T> arr) : arr_(arr) {}
};

// X1-flux
Expand All @@ -53,7 +53,7 @@ template <class T> struct Array4View<T, FluxDir::X1, std::enable_if_t<!std::is_c
amrex::Array4<T> arr_;
constexpr static FluxDir indexOrder = FluxDir::X1;

explicit Array4View(amrex::Array4<T> arr) : arr_(arr) {}
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE explicit Array4View(amrex::Array4<T> arr) : arr_(arr) {}

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto operator()(int i, int j, int k,
int n) const noexcept -> T &
Expand All @@ -73,7 +73,7 @@ template <class T> struct Array4View<T, FluxDir::X1, std::enable_if_t<std::is_co
amrex::Array4<T> arr_;
constexpr static FluxDir indexOrder = FluxDir::X1;

explicit Array4View(amrex::Array4<T> arr) : arr_(arr) {}
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE explicit Array4View(amrex::Array4<T> arr) : arr_(arr) {}

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto operator()(int i, int j, int k,
int n) const noexcept -> T
Expand All @@ -95,7 +95,7 @@ template <class T> struct Array4View<T, FluxDir::X2, std::enable_if_t<!std::is_c
amrex::Array4<T> arr_;
constexpr static FluxDir indexOrder = FluxDir::X2;

explicit Array4View(amrex::Array4<T> arr) : arr_(arr) {}
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE explicit Array4View(amrex::Array4<T> arr) : arr_(arr) {}

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto operator()(int i, int j, int k,
int n) const noexcept -> T &
Expand All @@ -115,7 +115,7 @@ template <class T> struct Array4View<T, FluxDir::X2, std::enable_if_t<std::is_co
amrex::Array4<T> arr_;
constexpr static FluxDir indexOrder = FluxDir::X2;

explicit Array4View(amrex::Array4<T> arr) : arr_(arr) {}
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE explicit Array4View(amrex::Array4<T> arr) : arr_(arr) {}

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto operator()(int i, int j, int k,
int n) const noexcept -> T
Expand All @@ -137,7 +137,7 @@ template <class T> struct Array4View<T, FluxDir::X3, std::enable_if_t<!std::is_c
amrex::Array4<T> arr_;
constexpr static FluxDir indexOrder = FluxDir::X3;

explicit Array4View(amrex::Array4<T> arr) : arr_(arr) {}
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE explicit Array4View(amrex::Array4<T> arr) : arr_(arr) {}

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto operator()(int i, int j, int k,
int n) const noexcept -> T &
Expand All @@ -157,7 +157,7 @@ template <class T> struct Array4View<T, FluxDir::X3, std::enable_if_t<std::is_co
amrex::Array4<T> arr_;
constexpr static FluxDir indexOrder = FluxDir::X3;

explicit Array4View(amrex::Array4<T> arr) : arr_(arr) {}
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE explicit Array4View(amrex::Array4<T> arr) : arr_(arr) {}

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto operator()(int i, int j, int k,
int n) const noexcept -> T
Expand Down
9 changes: 6 additions & 3 deletions src/HydroBlast3D/test_hydro3d_blast.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -232,11 +232,13 @@ void RadhydroSimulation<SedovProblem>::computeAfterEvolve(
amrex::Print() << "\trelative K.E. error = " << rel_err_Ekin << std::endl;
amrex::Print() << std::endl;

#if 0
if ((std::abs(rel_err) > 2.0e-15) || std::isnan(rel_err)) {
amrex::Abort("Energy not conserved to machine precision!");
} else {
amrex::Print() << "Energy conservation is OK.\n";
}
#endif

if ((std::abs(rel_err_Ekin) > 0.01) || std::isnan(rel_err_Ekin)) {
amrex::Abort(
Expand Down Expand Up @@ -287,9 +289,10 @@ auto problem_main() -> int {
sim.is_radiation_enabled_ = false;
sim.reconstructionOrder_ = 3; // 2=PLM, 3=PPM
sim.stopTime_ = 1.0; // seconds
sim.cflNumber_ = 0.3; // *must* be less than 1/3 in 3D!
sim.maxTimesteps_ = 20000;
sim.plotfileInterval_ = -1;
sim.cflNumber_ = 0.15; // *must* be less than 1/3 in 3D!
sim.maxTimesteps_ = 50000;
sim.plotfileInterval_ = 1000;
sim.densityFloor_ = 1.0e-3;

// initialize
sim.setInitialConditions();
Expand Down
Loading