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 "variables.h"
#include "vector.h"
#include "matrix.h"
#include "itersolv.h"

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

#include "parallel.h"
#include "gradient.h"
#include "geocalc.h"
#include "globals.h"
#include "setup.h"
#include "msolver.h"

#include "gamma.h"

double
CalculateMaxCourantNumber (double dt, int interface)
{

  unsigned int i, j;

  unsigned int element, face;

  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 (&xsl, element), 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) * faces[face].Aj * dt / elements[element].Vp, 0.0);
	  Cpp += cs * Cj;
	      
	}

      maxCp = LMAX (maxCp, Cpp);

      V_SetCmp (&Co, elements[element].index, Cpp);

    }
    
  VecAssemblyBegin (Co);
  VecAssemblyEnd (Co);

  VecGhostUpdateBegin (Co, INSERT_VALUES, SCATTER_FORWARD);
  VecGhostUpdateEnd (Co, INSERT_VALUES, SCATTER_FORWARD);

  return maxCp;

}

void
CorrectFaceS ()
{

  int i;

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

  msh_element ghost;

  double betaj;

  VecGhostGetLocalForm (xs, &xsl);
  
  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) > 0.0)
	    betaj = V_GetCmp (&betaf, face);
	  else
	    betaj = 1.0 - V_GetCmp (&betaf, face);

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

	}
      else
	{

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

	      ghost.index = faces[face].physreg;

	      ghost.celement.x = V_GetCmp (&cexl, faces[face].ghost);
	      ghost.celement.y = V_GetCmp (&ceyl, faces[face].ghost);
	      ghost.celement.z = V_GetCmp (&cezl, faces[face].ghost);

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

	      V_SetCmp (&xsf, face, LMAX (LMIN((1.0 - betaj) * V_GetCmp (&xsl, element) + betaj * V_GetCmp (&xsl, faces[face].ghost), 1.0), 0.0));

	    }
	  else
	    {

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

		  // zero gradient

		  V_SetCmp (&xsf, face, V_GetCmp (&xsl, element));

		}

	    }

	}

    }

  VecGhostRestoreLocalForm (xs, &xsl);
    
  VecAssemblyBegin (xsf);
  VecAssemblyEnd (xsf);

}

void
PredictBeta ()
{

  int i;

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

  double dot, l1, l2;

  double su, sdn, sjnCBC, sjnUQ, sjn;

  double qj, tetaj;

  double betaj;

  double Cod;

  double ang;

  msh_vector grads;

  VecGhostGetLocalForm (xs, &xsl);
  VecGhostGetLocalForm (Co, &Col);
  
  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) != 0.0)
		{

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

		      acceptor = neighbor;
		      donor = element;

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

		      dot = GeoDotVectorVector (grads, faces[face].d);

		      l1 = GeoMagVector (grads);

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

		    }
		  else
		    {

		      acceptor = element;
		      donor = neighbor;

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

		      dot = GeoDotVectorVector (grads, faces[pair].d);

		      l1 = GeoMagVector (grads);

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

		    }

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

		  Cod = LMIN (V_GetCmp (&Col, donor), 1.0);

		  if (LABS (V_GetCmp (&xsl, acceptor) - su) > SMALL)
		    {

		      sdn =
			(V_GetCmp (&xsl, donor) -
			 su) / (V_GetCmp (&xsl, acceptor) - su);

		      if (sdn >= 0.0 && sdn <= 1.0 && LABS (Cod) > SMALL)
			sjnCBC = LMIN (1, 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);
			}

		    }

		}

	    }
	  else
	    {

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

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

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

			  acceptor = faces[face].ghost;
			  donor = element;

			  grads.x = 0.0;
			  grads.y = 0.0;
			  grads.z = 0.0;

			  dot = GeoDotVectorVector (grads, faces[face].d);

			  l1 = GeoMagVector (grads);

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

			}
		      else
			{

			  acceptor = element;
			  donor = faces[face].ghost;

			  grads.x = 0.0;
			  grads.y = 0.0;
			  grads.z = 0.0;

			  dot = GeoDotVectorVector (grads, faces[face].d);

			  l1 = GeoMagVector (grads);

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

			}

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

		      Cod = LMIN (V_GetCmp (&Col, donor), 1.0);

		      if (LABS (V_GetCmp (&xsl, acceptor) - su) > SMALL)
			{

			  sdn = (V_GetCmp (&xsl, donor) - su) / (V_GetCmp (&xsl, acceptor) - su);

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

			  if (sdn >= 0.0 && sdn <= 1.0)
			    sjnUQ = LMIN ((8.0 * V_GetCmp (&Col, donor) * sdn + (1.0 - V_GetCmp (&Col, donor)) * (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, betaj);

    }

  VecGhostRestoreLocalForm (xs, &xsl);
  VecGhostRestoreLocalForm (Co, &Col);
    
  VecAssemblyBegin (betaf);
  VecAssemblyEnd (betaf);

}

void
CorrectBeta (double dt)
{

  int i;

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

  double Cj;

  double cbetaj, betaj;

  double ds, Ep, Em;

  VecGhostGetLocalForm (xs0, &xs0l);
  VecGhostGetLocalForm (xs, &xsl);
  
  for (i = 0; i < nbfaces; i++)
    {

      face = i;

      element = faces[face].element;

      pair = faces[face].pair;

      cbetaj = 0.0;

      betaj = V_GetCmp (&betaf, face);

      if (betaj < 1E-2) continue;

      if (pair != -1)
	{

	  neighbor = faces[pair].element;

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

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

		  acceptor = neighbor;
		  donor = element;

		}
	      else
		{

		  acceptor = element;
		  donor = neighbor;

		}

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

	      ds =
		0.5 * (V_GetCmp (&xs0l, acceptor) +
		       V_GetCmp (&xsl, acceptor)) - 0.5 *
		(V_GetCmp (&xs0l, donor) + V_GetCmp (&xsl, donor));

	      if (V_GetCmp (&xsl, donor) < 0.0)
		{
		  Em = LMAX (-V_GetCmp (&xsl, donor), 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 (&xsl, donor) > 1.0)
		{

		  Ep = LMAX (V_GetCmp (&xsl, donor) - 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);
			}
		    }
		}

	    }

	}
      else
	{

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

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

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

		      acceptor = faces[face].ghost;
		      donor = element;

		    }
		  else
		    {
		      acceptor = element;
		      donor = faces[face].ghost;

		    }

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

		  ds = 0.5 * (V_GetCmp (&xs0l, acceptor) + V_GetCmp (&xsl, acceptor)) - 0.5 * (V_GetCmp (&xs0l, donor) + V_GetCmp (&xsl, donor));

		  if (V_GetCmp (&xsl, donor) < 0.0)
		    {
		      Em = LMAX (-V_GetCmp (&xsl, donor), 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 (&xsl, donor) > 1.0)
		    {

		      Ep = LMAX (V_GetCmp (&xsl, donor) - 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, betaj);

    }

  VecGhostRestoreLocalForm (xs0, &xs0l);
  VecGhostRestoreLocalForm (xs, &xsl);

  VecAssemblyBegin (betaf);
  VecAssemblyEnd (betaf);

}

void
BoundScalar (Vec * x, Vec * xl, double min, double max)
{

  int i;

  int element;

  VecGhostGetLocalForm (*x, xl);

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

      element = i;

      V_SetCmp (x, elements[element].index, LMAX (LMIN (V_GetCmp (xl, element), max), min));

    }