gamma.c

OpenFVM-v1.1 open source cfd code

	    {

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

	      pair = faces[face].pair;

	      if (pair != -1)
		{
		  neighbor = faces[pair].element;

		  /*
		  dNf =
		    GeoMagVector (GeoSubVectorVector
				  (elements[neighbor].celement,
				   faces[pair].cface));
		  dPf =
		    GeoMagVector (GeoSubVectorVector
				  (elements[element].celement,
				   faces[face].cface));

		  lambda = dPf / (dPf + dNf);
		  */
		  
	          lambda = 0.5;

		  sj = sa[neighbor] * lambda + sa[element] * (1.0 - lambda);

		}
	      else
		{

		  sj = sa[element];

		}

	      sum1 += sj * faces[face].Aj;

	      sum2 += faces[face].Aj;

	    }

	  sm[i] = sum1 / sum2;

	}

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

	  element = i;

	  sa[i] = sm[i];
	}

    }

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

      element = i;

      V_SetCmp (xm, element + 1, sm[i]);

    }

  free (sa);
  free (sm);

}

void
CalculateMassFraction ()
{

  unsigned int i;

  register unsigned int element;

  double f[2];
  double vol[2];

  vol[0] = 0.0;
  vol[1] = 0.0;

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

      f[0] = (1.0 - V_GetCmp (&xs, element + 1));
      f[1] = V_GetCmp (&xs, element + 1);

      vol[0] += f[0] * V_GetCmp (&dens, element + 1) * elements[element].Vp;
      vol[1] += f[1] * V_GetCmp (&dens, element + 1) * elements[element].Vp;
    }

  printf ("\nMass of fluid %d: %+E kg\n", 0, vol[0]);
  printf ("\nMass of fluid %d: %+E kg\n", 1, vol[1]);

}

void
BuildVolumeOfFluidMatrix (double dt)
{

  unsigned int i, j, n;

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

  double aip;
  double ain[MAXFACES];
  unsigned int ani[MAXFACES];
  double bip;

  double betaj;

  // Equation: ds/dt + div(s*U) = 0

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

      element = i;

      aip = 0;

      bip = 0;

      n = 0;

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

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

	  pair = faces[face].pair;

	  if (parameter.scheme[is] == UDS)
	    {

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

	    }
	  else
	    {

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

	    }

	  if (pair != -1)
	    {

	      neighbor = faces[pair].element;

	      aip += 0.5 * (1.0 - betaj) * V_GetCmp (&uf, face + 1) * faces[face].Aj;

	      ain[n] = 0.5 * betaj * V_GetCmp (&uf, face + 1) * faces[face].Aj;

	      ani[n] = neighbor;
	      n++;

	      bip += -0.5 * (1.0 - betaj) * V_GetCmp (&uf, face + 1) * faces[face].Aj * V_GetCmp (&xs, element + 1);

	      bip += -0.5 * betaj * V_GetCmp (&uf, face + 1) * faces[face].Aj * V_GetCmp (&xs, neighbor + 1);

	    }
	  else
	    {

	      bip += -1.0 * V_GetCmp (&uf, face + 1) * faces[face].Aj * V_GetCmp (&xsf, face + 1);

	    }

	}

      if (dt > 0.0)
	{

	  aip += elements[element].Vp / dt;

	  bip += elements[element].Vp / dt * V_GetCmp (&xs0, element + 1);

	}

      if (aip == 0.0 || aip != aip)
	{
	  printf ("\nError: Problem setting up volume-of-fluid matrix\n");
	  exit (LOGICAL_ERROR);
	}

      Q_SetLen (&As, element + 1, n + 1);

      Q_SetEntry (&As, element + 1, 0, element + 1, aip);

      for (j = 0; j < n; j++)
	{
	  Q_SetEntry (&As, element + 1, j + 1, ani[j] + 1, ain[j]);
	}

      V_SetCmp (&bs, element + 1, bip);

    }

}

void
SolveVolumeOfFluidExplicit (double dt)
{

  unsigned int i, j, n;

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

  double aip;
  double ain[MAXFACES];
  unsigned int ani[MAXFACES];
  double bip;

  double betaj;

  double sums;

  // Equation: ds/dt + div(s*U) = 0

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

      element = i;

      aip = 0;

      bip = 0;

      n = 0;

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

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

	  pair = faces[face].pair;

	  if (parameter.scheme[is] == UDS)
	    {

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

	    }
	  else
	    {

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

	    }

	  if (pair != -1)
	    {

	      neighbor = faces[pair].element;

	      aip += 0.5 * (1.0 - betaj) * V_GetCmp (&uf, face + 1) * faces[face].Aj;

	      ain[n] = 0.5 * betaj * V_GetCmp (&uf, face + 1) * faces[face].Aj;

	      ani[n] = neighbor;
	      n++;

	      bip += -0.5 * (1.0 - betaj) * V_GetCmp (&uf, face + 1) * faces[face].Aj * V_GetCmp (&xs, element + 1);

	      bip += -0.5 * betaj * V_GetCmp (&uf, face + 1) * faces[face].Aj * V_GetCmp (&xs, neighbor + 1);

	    }
	  else
	    {

	      bip += -1.0 * V_GetCmp (&uf, face + 1) * faces[face].Aj * V_GetCmp (&xsf, face + 1);

	    }

	}

      if (dt > 0.0)
	{

	  aip += elements[element].Vp / dt;

	  bip += elements[element].Vp / dt * V_GetCmp (&xs0, element + 1);

	}

      sums = 0.0;

      for (j = 0; j < n; j++)
	{
	  sums += ain[j] * V_GetCmp (&xs0, ani[j] + 1);
	}

      V_SetCmp (&xs, element + 1, (bip - sums) / aip);

    }

}

void
CalculateGamma (char *var, int *fiter, double dt, double *maxCp, int verbose,
		int pchecks)
{

  unsigned int i, j;

  double mres;
  int miter;
  double mtime;

  if (parameter.calc[is] == LOGICAL_FALSE)
    return;

  if (parameter.vofastemp == LOGICAL_TRUE)
  {
     // The value of T is assigned to s

    Asgn_VV (&xs, &xT);
    return;

  } 

  *maxCp = CalculateMaxCourantNumber (dt, 1);

  parameter.ncicsamsteps = LMIN ((int) (*maxCp * 2) + 1, 100);

  if (verbose == LOGICAL_TRUE)
    printf ("\nMaximum Courant number at interface: %.3f\n", *maxCp);

  V_Constr (&betaf, "CICSAM interpolation factor", nbfaces, Normal, True);

  for (i = 0; i < parameter.ncicsamsteps; i++)
    {

      *maxCp = CalculateMaxCourantNumber (dt / parameter.ncicsamsteps, 0);

      // Store previous time step values 
      Asgn_VV (&xs0, &xs);

      fiter[is]++;

      // Predict beta - CICSAM       
      PredictBeta (&betaf, &xsf, &xs, &Co, &uf);

      for (j = 0; j <= parameter.ncicsamcor; j++)
	{

	  if (parameter.timemethod[is] == IMPLICITEULER)
	    {
	      Q_Constr (&As, "Indicator function matrix", nbelements, False,
			Rowws, Normal, True);
	      V_Constr (&bs, "Gamma source", nbelements, Normal, True);
	    }

	  if (parameter.timemethod[is] == EXPLICITEULER)
	    {

	      // Matrix free VOF 
	      SolveVolumeOfFluidExplicit (dt / parameter.ncicsamsteps);

	    }

	  if (parameter.timemethod[is] == IMPLICITEULER)
	    {

	      // Build VOF matrix
	      BuildVolumeOfFluidMatrix (dt / parameter.ncicsamsteps);

	      if (pchecks == LOGICAL_TRUE)
		{
		  if (!CheckIfDiagonalMatrix (&As))
		    {
		      printf
			("\nWarning: Volume-of-fluid matrix is not diagonal dominant\n");
		      WriteMatrix (&As, LOGICAL_FALSE);
		      WriteVector (&bs);
		      //exit (LOGICAL_ERROR);
		    }
		}

	      // Set matrix solution accuracy
	      SetRTCAccuracy (parameter.mtol[is]);

	      // Solve matrix to get indicator function s                   
	      SolveMatrix (&As, &xs, &bs, &miter, &mres, &mtime,
			   parameter.msolver[is], parameter.mprecond[is],
			   parameter.miter[is]);

	      if (verbose == LOGICAL_TRUE)
		printf
		  ("\nMatrix %c Number of iterations: %d Residual: %+E Time: %+E\n",
		   var[is], miter, mres, mtime);

	      if ((mres > parameter.mtol[is] && miter == parameter.miter[is])
		  || LASResult () != LASOK)
		{
		  printf ("\nError: Problem solving matrix %c\n", var[is]);
		  exit (LOGICAL_ERROR);
		}

	    }

	  // Correct beta 
	  CorrectBeta (dt / parameter.ncicsamsteps);

	  // Correct face and boundary 
	  CorrectFaceGamma ();

	  if (parameter.timemethod[is] == IMPLICITEULER)
	    {
	      Q_Destr (&As);
	      V_Destr (&bs);
	    }

	}

      if (parameter.fill == LOGICAL_TRUE)
	{
	  // Bound volume fraction
	  BoundScalar (&xs, 0.0, 1.0);
	}

      // Smooth volume fraction
      //SmoothScalar (&xsm, &xs, 2);

    }

  if (pchecks == LOGICAL_TRUE)
    {
      // Calculate mass fractions 
      CalculateMassFraction ();
    }

  V_Destr (&betaf);

}