pressure.c

OpenFVM-v1.1 open source cfd code

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

		  dj = GeoMagVector (GeoSubVectorVector
				     (ghost.celement,
				      elements[element].celement));

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

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

		  lambda = 0.5;

		  apj =
		    V_GetCmp (&apl, faces[face].ghost) * lambda +
		    V_GetCmp (&apl, element) * (1.0 - lambda);

		  Huj =
		    V_GetCmp (&hul, faces[face].ghost) * lambda +
		    V_GetCmp (&hul, element) * (1.0 - lambda);

		  Hvj =
		    V_GetCmp (&hvl, faces[face].ghost) * lambda +
		    V_GetCmp (&hvl, element) * (1.0 - lambda);

		  Hwj =
		    V_GetCmp (&hwl, faces[face].ghost) * lambda +
		    V_GetCmp (&hwl, element) * (1.0 - lambda);

		  Hf = Huj * faces[face].n.x +
		    Hvj * faces[face].n.y + Hwj * faces[face].n.z;

		  acp += -1.0 / (apj * dj) * faces[face].Aj;

		  acn[n] = 1.0 / (apj * dj) * faces[face].Aj;

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

		  V_SetCmp (&uf, face, 1.0 / apj * Hf);

		  bcp += V_GetCmp (&uf, face) * faces[face].Aj;

		  /*
		     if (parameter.orthof != 0.0)                  
		     gradpn = Gradient (&xpl, &xpf, LOGICAL_TRUE,  faces[face].ghost);

		     if (parameter.orthof != 0.0)
		     {
		     // Non-orthogonal correction term
		     bcp += -1.0 * parameter.orthof / (apj * dj) * faces[face].Aj *
		     (GeoDotVectorVector
		     (gradpn,
		     GeoSubVectorVector (faces[face].rnl,
		     ghost.celement)) -
		     GeoDotVectorVector (gradpp,
		     GeoSubVectorVector (faces[face].rpl,
		     elements[element].
		     celement))); 
		     }
		   */

		}
	      else
		{

		  apj = V_GetCmp (&apl, element);

		  Huj = V_GetCmp (&hul, element);
		  Hvj = V_GetCmp (&hvl, element);
		  Hwj = V_GetCmp (&hwl, element);

		  Hf =
		    Huj / apj * faces[face].n.x +
		    Hvj / apj * faces[face].n.y + Hwj / apj * faces[face].n.z;

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

		      acp +=
			-1.0 / (apj * (faces[face].dj + faces[face].kj)) *
			faces[face].Aj * (1.0 - V_GetCmp (&xsl, element));

		      bcp +=
			-1.0 / (apj * (faces[face].dj + faces[face].kj)) *
			faces[face].Aj * V_GetCmp (&xpf,
						   face) * (1.0 -
							    V_GetCmp (&xsl,
								      element));

		      V_SetCmp (&uf, face,
				1.0 / apj * Hf * (1.0 -
						  V_GetCmp (&xsl, element)));

		      bcp += V_GetCmp (&uf, face) * faces[face].Aj;

		      if (parameter.orthof != 0.0)
			{
			  // Non-orthogonal correction term               
			  bcp +=
			    1.0 * parameter.orthof / (apj *
						      (faces[face].dj +
						       faces[face].kj)) *
			    (GeoDotVectorVector
			     (gradpp,
			      GeoSubVectorVector (faces[face].rpl,
						  elements[element].
						  celement))) *
			    faces[face].Aj * (1.0 - V_GetCmp (&xsl, element));
			}

		    }

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

		      // velocity gradient = 0
		      // specified pressure

		      acp +=
			-1.0 / (apj * (faces[face].dj + faces[face].kj)) *
			faces[face].Aj;

		      bcp +=
			-1.0 / (apj * (faces[face].dj + faces[face].kj)) *
			faces[face].Aj * V_GetCmp (&xpf, face);

		      V_SetCmp (&uf, face, Hf / apj);

		      bcp += V_GetCmp (&uf, face) * faces[face].Aj;

		      if (parameter.orthof != 0.0)
			{
			  // Non-orthogonal correction term
			  bcp +=
			    1.0 * parameter.orthof / (apj *
						      (faces[face].dj +
						       faces[face].kj)) *
			    (GeoDotVectorVector
			     (gradpp,
			      GeoSubVectorVector (faces[face].rpl,
						  elements[element].
						  celement))) *
			    faces[face].Aj;
			}

		    }

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

		      // specified pressure 
		      // velocity gradient = 0

		      acp +=
			-1.0 / (apj * (faces[face].dj + faces[face].kj)) *
			faces[face].Aj;

		      bcp +=
			-1.0 / (apj * (faces[face].dj + faces[face].kj)) *
			faces[face].Aj * V_GetCmp (&xpf, face);

		      V_SetCmp (&uf, face, 1.0 / apj * Hf);

		      bcp += V_GetCmp (&uf, face) * faces[face].Aj;

		      if (parameter.orthof != 0.0)
			{
			  // Non-orthogonal correction term
			  bcp +=
			    1.0 * parameter.orthof / (apj *
						      (faces[face].dj +
						       faces[face].kj)) *
			    (GeoDotVectorVector
			     (gradpp,
			      GeoSubVectorVector (faces[face].rpl,
						  elements[element].
						  celement))) *
			    faces[face].Aj;
			}

		    }

		  if (faces[face].bc == INLET ||
		      faces[face].bc == MOVINGWALL ||
		      faces[face].bc == WALL ||
		      faces[face].bc == ADIABATICWALL ||
		      faces[face].bc == SURFACE)
		    {

		      // pressure gradient = 0
		      // specified velocity
		      V_SetCmp (&uf, face,
				V_GetCmp (&xuf, face) * faces[face].n.x +
				V_GetCmp (&xvf, face) * faces[face].n.y +
				V_GetCmp (&xwf, face) * faces[face].n.z);

		      bcp += V_GetCmp (&uf, face) * faces[face].Aj;

		    }

		}

	    }

	}

      if (acp == 0.0 || acp != acp)
	{
	  PetscPrintf (PETSC_COMM_WORLD, "\nError: Problem setting up continuity matrix\n");
	  exit (LOGICAL_ERROR);
	}
	
      /*
      Q_SetEntry (&Ac, elements[element].index, elements[element].index, acp);

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

      ncols = 0;
      nvals = 0;

      row = elements[element].index;

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

      val[nvals] = acp;
      nvals++;

      for (j = 0; j < n; j++)
	{
	  if (ani[j] > elements[element].index)
	    {
      		col[ncols] = ani[j];
	        ncols++;

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

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

      V_SetCmp (&bp, elements[element].index, bcp);

    }

  VecGhostRestoreLocalForm (ap, &apl);
  VecGhostRestoreLocalForm (hu, &hul);
  VecGhostRestoreLocalForm (hv, &hvl);
  VecGhostRestoreLocalForm (hw, &hwl);

  VecGhostRestoreLocalForm (xu, &xul);
  VecGhostRestoreLocalForm (xv, &xvl);
  VecGhostRestoreLocalForm (xw, &xwl);
  VecGhostRestoreLocalForm (xp, &xpl);
  VecGhostRestoreLocalForm (xs, &xsl);

  VecGhostRestoreLocalForm (dens, &densl);

  MatAssemblyBegin (Ac, MAT_FINAL_ASSEMBLY);
  MatAssemblyEnd (Ac, MAT_FINAL_ASSEMBLY);
  VecAssemblyBegin (bp);
  VecAssemblyEnd (bp);

}

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

  int i;

  double mres;
  int miter;
  double mtime;

  double presc;

  if (parameter.calc[ip] == LOGICAL_FALSE)
    return;
    
  V_Constr (&xpp, nbelements, 0);	// Pressure at cell center - previous iteration

  // Store previous time step values
  VecCopy (xp, xp0);

  fiter[ip]++;
   
  for (i = 0; i <= parameter.northocor; i++)
    {

      Q_Constr (&Ac, nbelements, LOGICAL_TRUE);
      V_Constr (&bp, nbelements, 0);	// Continuity source

      // Store previous iteration values
      VecCopy (xp, xpp);

      // Build the continuity matrix (mass conservation)          
      BuildContinuityMatrix (dt);

      if (pchecks == LOGICAL_TRUE)
	{
	  if (!CheckIfDiagonalMatrix (&Ac))
	    {
	      PetscPrintf (PETSC_COMM_WORLD,
			   "\nWarning: Continuity matrix is not diagonal dominant\n");
	      //MatView(Ac, PETSC_VIEWER_STDOUT_WORLD);
	      WriteMatrix (&Ac);
	      WriteVector (&bp);
	      //exit (LOGICAL_ERROR);
	    }
	}

      // Solve matrix to get pressure p
      SolveMatrix (&Ac, &xp, &bp, &miter, &mres, &mtime,
		   parameter.msolver[ip], parameter.mprecond[ip],
		   parameter.miter[ip], parameter.mtol[ip]);

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

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

      // Calculate pressure convergence
      VecWAXPY (temp1, -1.0, xp, xpp);
      VecNorm (temp1, NORM_2, &presc);

      if (verbose == LOGICAL_TRUE)
	PetscPrintf (PETSC_COMM_WORLD, "\nNon-orthogonality error (continuity): %+E\n", presc);

      if (presc < parameter.mtol[ip])
	break;

      CorrectFaceP ();

      Q_Destr (&Ac);
      V_Destr (&bp);

    }

  V_Destr (&xpp);
    
}