gamma.c

OpenFVM-v1.1 open source cfd code

  VecGhostRestoreLocalForm (*x, xl);

  VecAssemblyBegin (*x);
  VecAssemblyEnd (*x);

  VecGhostUpdateBegin (*x, INSERT_VALUES, SCATTER_FORWARD);
  VecGhostUpdateEnd (*x, INSERT_VALUES, SCATTER_FORWARD);
  
}

void
SmoothScalar (Vec * x, Vec * xl, int n)
{

  unsigned int i, j, k;

  register unsigned int face, pair;
  register unsigned int element, 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));

  VecGhostGetLocalForm (*x, xl);

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

      element = i;

      sa[i] = V_GetCmp (xl, element);

    }

  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++)
	    {

	      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++)
	{

	  sa[i] = sm[i];

	}

    }

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

      element = i;

      V_SetCmp (xl, element, sm[i]);

    }

  VecGhostRestoreLocalForm (*x, xl);

  free (sa);
  free (sm);

  VecAssemblyBegin (*x);
  VecAssemblyEnd (*x);

  VecGhostUpdateBegin (*x, INSERT_VALUES, SCATTER_FORWARD);
  VecGhostUpdateEnd (*x, INSERT_VALUES, SCATTER_FORWARD);

}

void
CalculateMassFraction ()
{

  unsigned int i;

  register unsigned int element;

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

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

  VecGhostGetLocalForm (xs, &xsl);
  VecGhostGetLocalForm (dens, &densl);
  
  for (i = 0; i < nbelements; i++)
    {
      element = i;

      f[1] = V_GetCmp (&xsl, element);
      f[0] = 1 - f[1];

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

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

  VecGhostRestoreLocalForm (xs, &xsl);
  VecGhostRestoreLocalForm (dens, &densl);
  
}

void
BuildVolumeOfFluidMatrix (double dt)
{

  unsigned int i, j, n;

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

  msh_element ghost;

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

  double betaj;

  // MatSetValues
  int row;
  int ncols;
  int col[MAXFACES+1];
  int nvals;
  double val[MAXFACES+1];

  VecGhostGetLocalForm (xs0, &xs0l);
  VecGhostGetLocalForm (xs, &xsl);
  
  // 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) > 0.0)
		betaj = 0.0;
	      else
		betaj = 1.0;

	    }
	  else
	    {

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

	    }

	  if (pair != -1)
	    {

	      neighbor = faces[pair].element;

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

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

	      ani[n] = elements[neighbor].index;
	      n++;

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

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

	    }
	  else
	    {

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

		  ghost.index = faces[face].physreg;

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

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

		  ani[n] = ghost.index;
		  n++;

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

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

		}
	      else
		{

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

	}

      if (dt > 0.0)
	{

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

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

	}

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

      /*
      Q_SetEntry (&As, elements[element].index, elements[element].index, aip);

      for (j = 0; j < n; j++)
	{
	  Q_SetEntry (&As, elements[element].index, ani[j], ain[j]);
	}
       */

      ncols = 0;
      nvals = 0;

      row = elements[element].index;

      col[ncols] = elements[element].index;
      ncols++;

      val[nvals] = aip;
      nvals++;

      for (j = 0; j < n; j++)
	{
		col[ncols] = ani[j];
	        ncols++;

		val[nvals] = ain[j];
	        nvals++;
	}

      Q_SetEntries (&As, 1, &row, ncols, col, val);

      V_SetCmp (&bs, elements[element].index, bip);

    }

  VecGhostRestoreLocalForm (xs0, &xs0l);
  VecGhostRestoreLocalForm (xs, &xsl);
    
  MatAssemblyBegin (As, MAT_FINAL_ASSEMBLY);
  MatAssemblyEnd (As, MAT_FINAL_ASSEMBLY);
  VecAssemblyBegin (bs);
  VecAssemblyEnd (bs);

}

void
SolveVolumeOfFluidExplicit (double dt)
{

  unsigned int i, j, n;

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

  msh_element ghost;

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

  double betaj;

  double sums;

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

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

  VecGhostGetLocalForm (xs0, &xs0l);
  VecGhostGetLocalForm (xs, &xsl);
  
  // 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) > 0.0)
		betaj = 0.0;
	      else
		betaj = 1.0;

	    }
	  else
	    {

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

	    }

	  if (pair != -1)
	    {

	      neighbor = faces[pair].element;

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

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

	      ani[n] = neighbor;
	      n++;

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

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

	    }
	  else
	    {

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

		  ghost.index = faces[face].physreg;

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

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

		  ani[n] = faces[face].ghost;
		  n++;

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

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

		}
	      else
		{

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

	}

      if (dt > 0.0)
	{

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

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

	}

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

      sums = 0.0;

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

      V_SetCmp (&xs, elements[element].index, (bip - sums) / aip);

    }

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

}

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

  int i, j;

  int loc;
  
  double mres;
  int miter;
  double mtime;

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

  CalculateMaxCourantNumber (dt, 1);

  VecMax (Co, &loc, maxCp);
  
  parameter.ncicsamsteps = LMIN ((int) (*maxCp * 2) + 1, 100);
 
  if (verbose == LOGICAL_TRUE)
    PetscPrintf (PETSC_COMM_WORLD, "\nMaximum Courant number at interface: %.3f\n", *maxCp);
  
  V_Constr (&betaf, nbfaces, 1);	// CICSAM interpolation factor

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

      CalculateMaxCourantNumber (dt / parameter.ncicsamsteps, 0);

      VecMax (Co, &loc, maxCp);

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

      fiter[is]++;

      // Predict beta - CICSAM       
      PredictBeta ();

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

	  if (parameter.timemethod[is] == IMPLICITEULER)
	    {
	      Q_Constr (&As, nbelements, LOGICAL_FALSE);
	      V_Constr (&bs, nbelements, 0);	// Gamma source
	    }

	  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))
		    {
		      PetscPrintf (PETSC_COMM_WORLD, "\nWarning: Volume-of-fluid matrix is not diagonal dominant\n");
		      //MatView(As, PETSC_VIEWER_STDOUT_WORLD);
		      WriteMatrix (&As);
		      WriteVector (&bs);
		      //exit (LOGICAL_ERROR);
		    }
		}

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

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

	      if (pchecks == LOGICAL_TRUE)
		{
		  if (mres > parameter.mtol[is]
		      && miter == parameter.miter[is])
		    {
		      PetscPrintf (PETSC_COMM_WORLD,
				   "\nError: Problem solving matrix %c\n",
				   var[is]);
		      exit (LOGICAL_ERROR);
		    }
		}

	    }

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

	  CorrectFaceS ();

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

	}

      if (parameter.fill == LOGICAL_TRUE)
	{
	  // Bound indicator function      
	  BoundScalar (&xs, &xsl, 0.0, 1.0);
	}

      //SmoothScalar (&xs, &xsl, 2);

    }

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

  V_Destr (&betaf);

}