/***********************************************************************************
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 flux_ct
*
* performs the flux-averaging used to preserve
* the del.B = 0 constraint (see Toth 2000)
*
*/
#include "decs.h"
void
flux_ct ()
{
int i, j;
static double emf[N1 + 1][N2 + 1];
double v1_av, w1, v2_av, w2, w_12;
double rho_av, sqrtrho, b1_av, b2_av, va1, va2;
double emfmm, emfpm, emfmp, emfpp;
/* calculate EMFs */
/* Toth approach: just average */
for (i = 0; i < N1 + 1; i++)
for (j = 0; j < N2 + 1; j++)
{
emf[i][j] = 0.25 * (F1[i][j][B2] + F1[i][j - 1][B2]
- F2[i][j][B1] - F2[i - 1][j][B1]);
}
/* rewrite EMFs as fluxes, after Toth */
for (i = 0; i < N1 + 1; i++)
for (j = 0; j < N2; j++)
{
F1[i][j][B1] = 0.;
F1[i][j][B2] = 0.5 * (emf[i][j] + emf[i][j + 1]);
}
for (i = 0; i < N1; i++)
for (j = 0; j < N2 + 1; j++)
{
F2[i][j][B1] = -0.5 * (emf[i][j] + emf[i + 1][j]);
F2[i][j][B2] = 0.;
}
}