/***********************************************************************************
Copyright 2008 Charles F. Gammie and Xiaoyue Guan
HAM version 1.0 (released Jan 15, 2008)
This file is part of HAM. HAM is a program that solves non-relativistic
hyperbolic partial differential equation in conservative form using
high-resolution shock-capturing techniques. This version of HAM has been
configured to solve the magnetohydrodynamic equations of motion in
axisymmetry to evolve a shearing box model.
You are morally obligated to cite the following paper in his/her
scientific literature that results from use of any part of HAM:
[1] Guan, X. \& Gammie, C.F. \ 2008, Astrophysical Journal, Supplement,
174, 145
The latest version of the HAM is available on our code distribution
website:
http://rainman.astro.uiuc.edu/codelib/
HAM 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 of the License, or
(at your option) any later version.
HAM 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 HAM; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
***********************************************************************************/
/*
* Subroutine slope-lim
*
* performs the slope-limiting for the numerical flux calculation
*
*/
#include "decs.h"
double
slope_lim (double y1, double y2, double y3)
{
double Dqm, Dqp, Dqc, s;
/* woodward, or monotonized central, slope limiter */
if (LIM == MC)
{
Dqm = 2. * (y2 - y1);
Dqp = 2. * (y3 - y2);
Dqc = 0.5 * (y3 - y1);
s = Dqm * Dqp;
if (s <= 0.)
return 0.;
else
{
if (fabs (Dqm) < fabs (Dqp) && fabs (Dqm) < fabs (Dqc))
return (Dqm);
else if (fabs (Dqp) < fabs (Dqc))
return (Dqp);
else
return (Dqc);
}
}
/* van leer slope limiter */
else if (LIM == VANL)
{
Dqm = (y2 - y1);
Dqp = (y3 - y2);
s = Dqm * Dqp;
if (s <= 0.)
return 0.;
else
return (2. * s / (Dqm + Dqp));
}
/* minmod slope limiter (crude but robust) */
else if (LIM == MINM)
{
Dqm = (y2 - y1);
Dqp = (y3 - y2);
s = Dqm * Dqp;
if (s <= 0.)
return 0.;
else if (fabs (Dqm) < fabs (Dqp))
return Dqm;
else
return Dqp;
}
fprintf (stderr, "unknown slope LIMiter\n");
exit (10);
return (0.);
}