Commit 054eb4c4 authored by Philippe Helluy's avatar Philippe Helluy

end review 1

parent bb5dac00
......@@ -1237,7 +1237,10 @@ the transport operator then reduces to a simple shift.
Before the collision step, we have thus
\begin_inset Formula
\begin{equation}
\v f_{1,i,j}^{n+1,-}=\v f_{1,i+1,j}^{n},\quad\v f_{2,i,j}^{n+1,-}=\v f_{2,i-1,j}^{n},\quad\v f_{3,i,j}^{n+1,-}=\v f_{3,i,j+1}^{n},\quad\v f_{4,i,j}^{n+1,-}=\v f_{4,i,j-1}^{n}.\label{eq:shift-algorithm}
\begin{array}{c}
\v f_{1,i,j}^{n+1,-}=\v f_{1,i+1,j}^{n},\v f_{2,i,j}^{n+1,-}=\v f_{2,i-1,j}^{n},\\
\v f_{3,i,j}^{n+1,-}=\v f_{3,i,j+1}^{n},\quad\v f_{4,i,j}^{n+1,-}=\v f_{4,i,j-1}^{n}.
\end{array}\label{eq:shift-algorithm}
\end{equation}
\end_inset
......
......@@ -75,8 +75,9 @@
\begin_body
\begin_layout Standard
Answers to reviewer 1 (most of the changes are highlighted in green)
\begin_layout Subsection*
Answers to reviewer 1 (most of the corresponding changes are highlighted
in green)
\end_layout
\begin_layout Enumerate
......@@ -196,29 +197,36 @@ One of the authors deeply studied Riemann problems for MHD with Godunov
\begin_inset Newline newline
\end_inset
We have written an additional remark at the end of section 3.2
\end_layout
\begin_layout Enumerate
In the numerical tests you choose epsilon=10^{-3}.
What happens for different scales of this value?
\begin_inset Newline newline
\end_inset
We have not tried other values of this parameter.
epsilon has been fixed to this value for comparison with physics papers.
The objective here is really to catch the linear behavior of the instability.
\end_layout
\begin_layout Standard
Other comments:
\end_layout
\begin_layout Enumerate
\begin_inset Quotes eld
\end_inset
Please write MagnetoHydroDynamic extensively at the beginning of the paper
and then replace it by the acronym.
The same applies to GPU (in the abstract is used the acronym, but never
introduced).
\begin_inset Quotes erd
\end_inset
\begin_layout Standard
We have corrected all the signaled typos.
Thanks for pointing then to us.
\end_layout
\begin_layout Subsection*
Answers to reviewer 2 (most of the corresponding changes are highlighted
in red)
\end_layout
\begin_inset Newline newline
\end_inset
\begin_layout Standard
We have corrected the writing at several places.
\end_layout
\end_body
......
......@@ -10,6 +10,17 @@
year={2020}
}
@article{bouchut2010multiwave,
title={A multiwave approximate Riemann solver for ideal MHD based on relaxation II: numerical implementation with 3 and 5 waves},
author={Bouchut, Fran{\c{c}}ois and Klingenberg, Christian and Waagan, Knut},
journal={Numerische Mathematik},
volume={115},
number={4},
pages={647--679},
year={2010},
publisher={Springer}
}
@article{klockner2012pycuda,
title={PyCUDA and PyOpenCL: A scripting-based approach to GPU run-time code generation},
author={Kl{\"o}ckner, Andreas and Pinto, Nicolas and Lee, Yunsup and Catanzaro, Bryan and Ivanov, Paul and Fasih, Ahmed},
......
......@@ -106,7 +106,7 @@
\abstract{This paper is devoted to the simulation of \revA{Magnetohydrohynamics} flows with complex
structures. This kind flows present instabilities that generate shock
waves. We propose a robust and precise numerical method based on the
waves. We propose a robust and accurate numerical method based on the
Lattice Boltzmann methodology. We explain how to adjust the numerical
viscosity in order to obtain stable, \revA{accurate results in smooth or discontinuous parts of the flow} and reduced divergence
errors. This method can handle shock wave and is almost second order.
......@@ -148,7 +148,7 @@
on a regular grid. An important aspect of the LBM is the choice of
the relaxation parameter in the collision step. The choice of the
parameter allows adjusting the numerical viscosity of the LBM scheme.
We provide a analysis in a simplified one-dimensional framework which
We provide an analysis in a simplified one-dimensional framework which
shows that it is possible to adjust more precisely the numerical viscosity
with a generalized matrix relaxation parameter. We also show that
the divergence cleaning effect is improved if the relaxation parameter
......@@ -237,9 +237,9 @@
\[
\v u=(u_{1},u_{2},u_{3})^{T},\quad\v B=(B_{1},B_{2},B_{3})^{T},
\]
the pressure is defined by
the pressure is given \revA{by a perfect-gas law with a constant polytropic exponent $\gamma>1$ }
\[
p=(\gamma-1)(Q-\rho\frac{\v u\cdot\v u}{2}-\frac{\v B\cdot\v B}{2}).
p=(\gamma-1)(Q-\rho\frac{\v u\cdot\v u}{2}-\frac{\v B\cdot\v B}{2}),
\]
The other variables are the density $\rho$, the total energy $Q$,
......@@ -424,7 +424,7 @@
of the original equations. See for instance \cite{dellar2013interpretation}).
See also below for a mathematical analysis in the one-dimensional
case.
\section{vide}
%\section{vide}
\section{Numerical method}
......@@ -474,8 +474,11 @@ $
$
the transport operator then reduces to a simple shift. Before the
collision step, we have thus
\begin{equation}
\v f_{1,i,j}^{n+1,-}=\v f_{1,i+1,j}^{n},\quad\v f_{2,i,j}^{n+1,-}=\v f_{2,i-1,j}^{n},\quad\v f_{3,i,j}^{n+1,-}=\v f_{3,i,j+1}^{n},\quad\v f_{4,i,j}^{n+1,-}=\v f_{4,i,j-1}^{n}.\label{eq:shift-algorithm}
\begin{equation}\label{eq:shift-algorithm}\left\{
\begin{array}{l}
\v f_{1,i,j}^{n+1,-}=\v f_{1,i+1,j}^{n},\v f_{2,i,j}^{n+1,-}=\v f_{2,i-1,j}^{n},\\
\v f_{3,i,j}^{n+1,-}=\v f_{3,i,j+1}^{n},\quad\v f_{4,i,j}^{n+1,-}=\v f_{4,i,j-1}^{n}.
\end{array}\right.
\end{equation}
......@@ -490,7 +493,7 @@ $
\begin{equation}
\v f_{k,i,j}^{n+1}=\v f_{k}^{eq}(\v w_{i,j}^{n+1}).\label{eq:equilibrium}
\end{equation}
As stated above, it is more precise to consider an over-relaxation
As stated above, it is more accurate to consider an over-relaxation
\[
\v f_{k,i,j}^{n+1}=2\v f_{k}^{eq}(\v w_{i,j}^{n+1})-\v f_{k,i,j}^{n+1,-}.
\]
......@@ -505,6 +508,7 @@ $
\end{equation}
where $\omega=\frac{2\Delta t}{2\tau+\Delta t}$. Typically, one will
choose $\omega=1.9$ in our simulations.
\revA{We see that the whole algorithm is extremely simple. It is a succession of "shifts" (\ref{eq:shift-algorithm}) and "collisions" (\ref{eq:collision}. For $\omega=2$ the scheme is second order but unstable in shocks. For $\omega=1$ the scheme is very robust, entropy dissipative, but quite diffusive. It would be interesting to construct a rigorous strategy for choosing locally the optimal value of $\omega$. Let us mention that many finite volume schemes have been designed for solving MHD equations, with specific Riemann solvers. We can mention for instance the solver of \cite{bouchut2010multiwave}, based on a relaxation approach, with proven stability and accuracy features. However, the solver presented in \cite{bouchut2010multiwave} is more complicated to program and less computationally efficient. We can thus compensate the slightly lower accuracy of the kinetic scheme by finer meshes, without increasing too much the computational load.}
\subsection{Boundary conditions}
......@@ -531,7 +535,7 @@ $
For the analysis, it is possible to replace the scalar relaxation
parameter $\omega$ by a matrix $\vom$, for more generality. In the
following, we understand the comparison of matrices in the usual way,
following, we \revA{establish the comparison between matrices} in the usual way,
by the comparison of the associated quadratic form (the resulting
order is thus not total).
......@@ -581,7 +585,7 @@ $
In order to check practically the accuracy of the approximation, we
apply the above analysis for a simplified system of two conservation
laws
laws.
We consider the one-dimensional isothermal Euler equations with a
diagonal diffusion of $\epsilon\partial_{xx}\left(\rho,\rho u\right)^{T}$.
......@@ -597,7 +601,7 @@ $
\end{align}
This system is a very simplified version of the MHD equation, where
the magnetic field is assumed to vanish and the gas is supposed to
the magnetic field is assumed to vanish and the gas is supposed to be
isothermal. We have chosen to use this specific non-physical diffusion
for our first test since it is exactly the type of diffusion that
is apparent in a standard finite volume code using a Lax-Friedrichs
......@@ -635,15 +639,15 @@ $
For our first test of the matrix relaxation, we solve (\ref{eq: Isothermal Navier-Stokes})
comparing the Lattice Boltzmann sheme using (\ref{eq:relax-with-diffusion})
as relaxation matrix and a standard explicit centered finite volume
scheme for approximating (\eqref{eq: Isothermal Navier-Stokes}).
scheme for approximating (\ref{eq: Isothermal Navier-Stokes}).
In this centered scheme, the time step is taken very small in such
way that the stability condition is satisfied and that the time integration
error can be neglected. In other words, the equivalent PDE of both
schemes is (\ref{eq:equiv-eq-2order}). Hence, given the parameters
in both schemes are set to represent the same diffusion $\epsilon$,
one should get the same type of diffusion for both schemes.\\
To test this we take the simple test case of a stationary viscous
shock:\\
To test this we take the simple case of a stationary viscous
shock.\\
The initial data for this shock tube problem are
\begin{align}
\wvec_{l}=\left(\rho_{l},\rho_{l}u_{l}\right)^{T}\ \ \ \wvec_{r}=\left(\frac{\rho_{l}u_{l}^{2}}{c},\frac{c}{u_{l}}\right)^{T}
......@@ -657,7 +661,7 @@ $
scheme we set $\lambda=40$. The results on a $128$ cell grid at
time $T=0.1$ are shown in Figure \ref{Fg: ViscouseShock}. The test
shows that the matrix relaxation lattice Boltzmann scheme provides
the correct diffusion and therefore results in the right viscose profile.
the correct diffusion and therefore results in the right viscous profile.
\begin{figure}[h!]
\centering{}\includegraphics[width=0.4\textwidth]{png/FVvsMatRel}\includegraphics[width=0.4\textwidth]{png/FVvsMatRelZoom}
\caption{{Viscous Shock Test with $\epsilon=0.1$: Comparison of the viscous
......@@ -1068,7 +1072,7 @@ $
0
\end{matrix}\right),
\]
with $\epsilon=1.0\cdot10^{-3}$.
with $\epsilon=10^{-3}$.
The asymptotic magnetic field strength for large $r$ is unity, and
thus defines our normalization. We also point out that we consider
......@@ -1212,7 +1216,6 @@ $
of the kinetic energy growth rate according to the grid and the order
of the numerical method. The converged value for the growth rate is
near 1.3 for an initial perturbation $\epsilon=1.0\cdot10^{-3}$.
In
The simulations of the tilt instability that we perform are post-processed
at regular time intervals. In Figures \ref{fig:kinetic-1} and \ref{fig:kinetic-2},
......@@ -1231,9 +1234,9 @@ $
The converged growth rate is close to 1.45 for our simulations. It
is higher than the results of \cite{richard1990tilt} and \cite{lankalapalli2007adaptive}.
In \cite{keppens2014interacting} for a similar configuration the
growth rate is evaluated to 1.498. This probably very close to the
exact rate because the simulation is conducted with a very precise
adaptive scheme. Our results seem thus to be quite precise on the
growth rate is evaluated to 1.498. This is probably very close to the
exact rate because the simulation is conducted with a very accurate
adaptive scheme. Our results seem thus to be quite correct on the
finest mesh.
%\tikzsetnextfilename{kinetic_512-dirichlet}
......@@ -1406,6 +1409,6 @@ $
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\bibliographystyle{plain}
\bibliography{mhd_lbm.bib}
\bibliography{spring_mhd_lbm.bib}
\end{document}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment