gamma.c

OpenFVM-v1.1 open source cfd code

/***************************************************************************
 *   Copyright (C) 2004-2008 by OpenFVM team                               *
 *   http://sourceforge.net/projects/openfvm/                              *
 *                                                                         *
 *   This program 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.                                   *
 *                                                                         *
 *   This program 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 this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

#include "globals.h"

#include "mesh.h"
#include "material.h"
#include "bcond.h"
#include "param.h"

#include "variables.h"
#include "gradient.h"
#include "geocalc.h"
#include "setup.h"
#include "msolver.h"

#include "gamma.h"

double
CalculateMaxCourantNumber (double dt, int interface)
{

  unsigned int i, j;

  register unsigned int face;
  register unsigned int element;

  double Cpp;
  double Cj;

  double s, cs;
  
  double maxCp;

  maxCp = 0.0;

  for (i = 0; i < nbelements; i++)
    {

      element = i;

      if (interface == 1)
      {     
      	s = LMIN(LMAX(V_GetCmp (&xs, element + 1), 0.0), 1.0);
      	cs = (1.0 - s) * (1.0 - s) * s * s * 16.0;    
      }
      else
      {
        cs = 1.0;
      }
      
      Cpp = 0.0;

      for (j = 0; j < elements[element].nbfaces; j++)
	{

	  face = elements[element].face[j];

          Cj = LMAX (-V_GetCmp (&uf, face + 1) * faces[face].Aj * dt / elements[element].Vp, 0.0);
          Cpp += cs * Cj;

	}

      maxCp = LMAX (maxCp, Cpp);

      V_SetCmp (&Co, element + 1, Cpp);

    }

  return maxCp;

}

void
CorrectFaceGamma ()
{

  int i;

  register unsigned int face, pair;
  register unsigned int element, neighbor;

  double betaj;

  for (i = 0; i < nbfaces; i++)
    {

      face = i;

      element = faces[face].element;

      pair = faces[face].pair;

      if (pair != -1)
	{

	  neighbor = faces[pair].element;

	  if (V_GetCmp (&uf, face + 1) > 0.0)
	    betaj = V_GetCmp (&betaf, face + 1);
	  else
	    betaj = 1.0 - V_GetCmp (&betaf, face + 1);

	  V_SetCmp (&xsf, face + 1,
		    LMAX (LMIN
			  ((1.0 - betaj) * V_GetCmp (&xs,
						     element + 1) +
			   betaj * V_GetCmp (&xs, neighbor + 1), 1.0), 0.0));

	}
      else
	{

	  if (faces[face].bc == OUTLET)
	    {

	      // zero gradient

	      V_SetCmp (&xsf, face + 1, V_GetCmp (&xs, element + 1));

	    }

	}

    }

}

void
PredictBeta ()
{

  unsigned int i;

  register unsigned int face, pair;
  register unsigned int element, neighbor;
  register unsigned int donor, acceptor;

  double dot, l1, l2;

  double su, sdn;
  double sjnCBC, sjnUQ;
  double sjn;

  double qj, tetaj;

  double betaj;

  double Cod;

  double ang;

  msh_vector grads;

  for (i = 0; i < nbfaces; i++)
    {

      face = i;

      element = faces[face].element;

      pair = faces[face].pair;

      betaj = 0.0;

      if (parameter.ncicsamcor != 0)
	{

	  if (pair != -1)
	    {

	      neighbor = faces[pair].element;

	      if (V_GetCmp (&uf, face + 1) != 0.0)
		{

		  if (V_GetCmp (&uf, face + 1) > 0.0)
		    {

		      acceptor = neighbor;
		      donor = element;

		      grads = Gradient (&xs, &xsf, LOGICAL_FALSE, donor);

		      dot = grads.x * faces[face].d.x + grads.y * faces[face].d.y + grads.z * faces[face].d.z;

		      l1 = GeoMagVector (grads);

		      l2 = GeoMagVector (faces[face].d);

		    }
		  else
		    {

		      acceptor = element;
		      donor = neighbor;

		      grads = Gradient (&xs, &xsf, LOGICAL_FALSE, donor);

		      dot = grads.x * faces[pair].d.x + grads.y * faces[pair].d.y + grads.z * faces[pair].d.z;

		      l1 = GeoMagVector (grads);

		      l2 = GeoMagVector (faces[pair].d);

		    }

		  su = LMIN (LMAX (V_GetCmp (&xs, acceptor + 1) - 2 * dot, 0.0), 1.0);

		  Cod = LMIN(V_GetCmp (&Co, donor + 1), 1.0);

		  if (LABS (V_GetCmp (&xs, acceptor + 1) - su) > SMALL)
		    {

		      sdn = (V_GetCmp (&xs, donor + 1) - su) / (V_GetCmp (&xs, acceptor + 1) - su);

		      if (sdn >= 0.0 && sdn <= 1.0 && LABS (Cod) > SMALL)
			  sjnCBC = LMIN (1.0, sdn / Cod);
		      else
			sjnCBC = sdn;

		      if (sdn >= 0.0 && sdn <= 1.0)
			  sjnUQ =  LMIN ((8.0 * Cod * sdn + (1.0 - Cod) * (6.0 * sdn + 3.0)) / 8.0, sjnCBC);
		      else
			  sjnUQ = sdn;

		      if (LABS (l1 * l2) > SMALL)
			ang = LABS (dot / (l1 * l2));
		      else
			ang = LABS (dot / SMALL);

		      if (ang > 1.0)
			ang = 1.0;

		      tetaj = acos (ang);

       	              qj = LMIN (parameter.kq * 0.5 * (cos (2 * tetaj) + 1.0), 1.0);

		      sjn = qj * sjnCBC + (1 - qj) * sjnUQ;

		      if (LABS (1.0 - sdn) > SMALL)
			{
			  betaj = LMIN (LMAX ((sjn - sdn) / (1.0 - sdn), 0.0), 1.0);
			}
		    }

		}

	    }

	}

      V_SetCmp (&betaf, face + 1, betaj);

    }

}

void
CorrectBeta (double dt)
{

  unsigned int i;

  unsigned int face, pair;
  register unsigned int element;
  unsigned int neighbor;
  unsigned int donor, acceptor;

  double Cj;

  double cbetaj, betaj;

  double ds, Ep, Em;

  for (i = 0; i < nbfaces; i++)
    {

      face = i;

      element = faces[face].element;

      pair = faces[face].pair;

      cbetaj = 0.0;

      betaj = V_GetCmp (&betaf, face + 1);

      if (betaj < 1E-2) continue;

      if (pair != -1)
	{

	  neighbor = faces[pair].element;

	  if (V_GetCmp (&uf, face + 1) != 0.0)
	    {

	      if (V_GetCmp (&uf, face + 1) > 0.0)
		{

		  acceptor = neighbor;
		  donor = element;

		}
	      else
		{
		  acceptor = element;
		  donor = neighbor;

		}

	      Cj =
		LMIN (LMAX
		      (-V_GetCmp (&uf, face + 1) * faces[face].Aj * dt /
		       elements[element].Vp, 0.0), 1.0);

	      ds =
		0.5 * (V_GetCmp (&xs0, acceptor + 1) +
		       V_GetCmp (&xs, acceptor + 1)) -
		0.5 * (V_GetCmp (&xs0, donor + 1) +
		       V_GetCmp (&xs, donor + 1));

	      if (V_GetCmp (&xs, donor + 1) < 0.0)
		{
		  Em = LMAX (-V_GetCmp (&xs, donor + 1), 0.0);

		  // Donor value < 0.0 Ex: sd = -0.1 -> Em = +0.1 
		  if (Em > SMALL && Cj > SMALL)
		    {
		      if (ds > Em)
			{
			  cbetaj =
			    Em * (2 + Cj -
				  2 * Cj * betaj) / (2 * Cj * (ds - Em));

			  cbetaj = LMIN (cbetaj, betaj);
			}
		    }

		}


	      if (V_GetCmp (&xs, donor + 1) > 1.0)
		{

		  Ep = LMAX (V_GetCmp (&xs, donor + 1) - 1.0, 0.0);

		  // Donor value > 1.0 Ex: sd = 1.1 -> Ep = +0.1 
		  if (Ep > SMALL && Cj > SMALL)
		    {
		      if (ds < -Ep)
			{
			  cbetaj =
			    Ep * (2 + Cj -
				  2 * Cj * betaj) / (2 * Cj * (-ds - Ep));

			  cbetaj = LMIN (cbetaj, betaj);
			}
		    }
		}

	    }

	}

      betaj -= cbetaj;

      betaj = LMAX (betaj, 0.0);

      V_SetCmp (&betaf, face + 1, betaj);

    }

}

void
BoundScalar (Vector * x, double min, double max)
{

  unsigned int i;

  register unsigned int element;

  for (i = 0; i < nbelements; i++)
    {

      element = i;

      V_SetCmp (x, element + 1,
		LMAX (LMIN (V_GetCmp (x, element + 1), max), min));
    }

}

void
SmoothScalar (Vector * xm, Vector * x, int n)
{

  unsigned int i, j, k;

  unsigned int face, pair;

  register unsigned int element;
  unsigned int neighbor;

  double sj;

  //double dNf, dPf;
  double lambda;

  double sum1, sum2;

  double *sa;
  double *sm;

  sa = calloc (nbelements, sizeof (double));
  sm = calloc (nbelements, sizeof (double));

  for (i = 0; i < nbelements; i++)
    {

      element = i;

      sa[i] = V_GetCmp (x, element + 1);
      sm[i] = 0.0;

    }

  for (k = 0; k < n; k++)
    {

      for (i = 0; i < nbelements; i++)
	{

	  element = i;

	  sum1 = 0.0;

	  sum2 = 0.0;

	  for (j = 0; j < elements[element].nbfaces; j++)