velocity.c

OpenFVM-v1.1 open source cfd code

	      // Convection           
	      if (parameter.scheme[iu] == UDS)
		{
		  // UDS
		  if (V_GetCmp (&uf, face) > 0.0)
		    xsi = 0.0;
		  else
		    xsi = 1.0;

		}
	      else
		{
		  // CDS
		  xsi = lambda;
		}

	      // Convection
	      app += (1.0 - xsi) * densp * V_GetCmp (&uf, face) * faces[face].Aj / elements[element].Vp;

	      // Diffusion
	      app += viscj * faces[face].Aj / (faces[face].dj + faces[face].kj) / elements[element].Vp;

	      // Convection
	      apn[n] = xsi * densp * V_GetCmp (&uf, face) * faces[face].Aj / elements[element].Vp;

	      // Diffusion
	      apn[n] += -viscj * faces[face].Aj / (faces[face].dj + faces[face].kj) / elements[element].Vp;

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

	      /*
	         // Non-orthogonal correction terms
	         if (parameter.orthof != 0.0)
	         {
	         gradun = Gradient (&xul, &xuf, LOGICAL_TRUE, neighbor);
	         gradvn = Gradient (&xvl, &xvf, LOGICAL_TRUE, neighbor);
	         gradwn = Gradient (&xwl, &xwf, LOGICAL_TRUE, neighbor);

	         bpu += parameter.orthof * 
	         -viscj * faces[face].Aj / (faces[face].dj + faces[face].kj) /
	         elements[element].Vp *
	         (GeoDotVectorVector
	         (gradun,
	         GeoSubVectorVector (faces[face].rnl,
	         elements[neighbor].celement)) -
	         GeoDotVectorVector (gradup,
	         GeoSubVectorVector (faces[face].rpl,
	         elements[element].
	         celement)));

	         bpv += parameter.orthof * 
	         -viscj * faces[face].Aj / (faces[face].dj + faces[face].kj) /
	         elements[element].Vp *
	         (GeoDotVectorVector
	         (gradvn,
	         GeoSubVectorVector (faces[face].rnl,
	         elements[neighbor].celement)) -
	         GeoDotVectorVector (gradvp,
	         GeoSubVectorVector (faces[face].rpl,
	         elements[element].
	         celement)));

	         bpw += parameter.orthof * 
	         -viscj * faces[face].Aj / (faces[face].dj + faces[face].kj) /
	         elements[element].Vp *
	         (GeoDotVectorVector
	         (gradwn,
	         GeoSubVectorVector (faces[face].rnl,
	         elements[neighbor].celement)) -
	         GeoDotVectorVector (gradwp,
	         GeoSubVectorVector (faces[face].rpl,
	         elements[element].
	         celement)));        
	         }
	       */

	    }
	  else
	    {

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

		  // Face between processors
		  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);

		  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;

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

		  if (parameter.scheme[iu] == UDS)
		    {
		      // UDS
		      if (V_GetCmp (&uf, face) > 0.0)
			xsi = 0.0;
		      else
			xsi = 1.0;

		    }
		  else
		    {
		      // CDS
		      xsi = lambda;
		    }

		  // Convection
		  app += (1.0 - xsi) * densp * V_GetCmp (&uf, face) * faces[face].Aj / elements[element].Vp;

		  // Diffusion
		  app += viscj * faces[face].Aj / dj / elements[element].Vp;

		  // Convection
		  apn[n] = xsi * densp * V_GetCmp (&uf, face) * faces[face].Aj / elements[element].Vp;

		  // Diffusion
		  apn[n] += -viscj * faces[face].Aj / dj / elements[element].Vp;

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

		  /*
		     // Non-orthogonal correction terms
		     if (parameter.orthof != 0.0)
		     {
		     gradun = Gradient (&xul, &xuf, LOGICAL_TRUE, faces[face].ghost);    
		     gradvn = Gradient (&xvl, &xvf, LOGICAL_TRUE, faces[face].ghost);
		     gradwn = Gradient (&xwl, &xwf, LOGICAL_TRUE, faces[face].ghost);

		     bpu += parameter.orthof * 
		     -viscj * faces[face].Aj / dj /
		     elements[element].Vp *
		     (GeoDotVectorVector
		     (gradun,
		     GeoSubVectorVector (faces[face].rnl,
		     ghost.celement)) -
		     GeoDotVectorVector (gradup,
		     GeoSubVectorVector (faces[face].rpl,
		     elements
		     [element].
		     celement)));

		     bpv += parameter.orthof * 
		     -viscj * faces[face].Aj / dj /
		     elements[element].Vp *
		     (GeoDotVectorVector
		     (gradvn,
		     GeoSubVectorVector (faces[face].rnl,
		     ghost.celement)) -
		     GeoDotVectorVector (gradvp,
		     GeoSubVectorVector (faces[face].rpl,
		     elements
		     [element].
		     celement)));

		     bpw += parameter.orthof * 
		     -viscj * faces[face].Aj / dj /
		     elements[element].Vp *
		     (GeoDotVectorVector
		     (gradwn,
		     GeoSubVectorVector (faces[face].rnl,
		     ghost.celement)) -
		     GeoDotVectorVector (gradwp,
		     GeoSubVectorVector (faces[face].rpl,
		     elements
		     [element].
		     celement)));
		     }
		   */

		}
	      else
		{

		  if (faces[face].bc != EMPTY)
		    {

		      viscj = V_GetCmp (&viscl, element);

		      if ((faces[face].bc != PERMEABLE || V_GetCmp (&xsl, element) > 0.5) && faces[face].bc != SLIP)
		        {
		           // Diffusion
		           app += viscj * faces[face].Aj / (faces[face].dj + faces[face].kj) / elements[element].Vp;

		           bpu += viscj * faces[face].Aj / (faces[face].dj + faces[face].kj) * V_GetCmp (&xuf, face) / elements[element].Vp;
		           bpv += viscj * faces[face].Aj / (faces[face].dj + faces[face].kj) * V_GetCmp (&xvf, face) / elements[element].Vp;
		           bpw += viscj * faces[face].Aj / (faces[face].dj + faces[face].kj) * V_GetCmp (&xwf, face) / elements[element].Vp;
 		        }

		      // Convection
		      if (parameter.scheme[iu] == UDS)
			{

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

			      // Convection
			      app += densp * V_GetCmp (&uf, face) * faces[face].Aj / elements[element].Vp;

			    }
			  else
			    {

			      // Convection
			      bpu += -densp * V_GetCmp (&uf, face) * faces[face].Aj * V_GetCmp (&xuf, face) / elements[element].Vp;
			      bpv += -densp * V_GetCmp (&uf, face) * faces[face].Aj * V_GetCmp (&xvf, face) / elements[element].Vp;
			      bpw += -densp * V_GetCmp (&uf, face) * faces[face].Aj * V_GetCmp (&xwf, face) / elements[element].Vp;

			    }

			}
		      else
			{

			  // CDS
			  bpu += -densp * V_GetCmp (&uf, face) * faces[face].Aj * V_GetCmp (&xuf, face) / elements[element].Vp;
			  bpv += -densp * V_GetCmp (&uf, face) * faces[face].Aj * V_GetCmp (&xvf, face) / elements[element].Vp;
			  bpw += -densp * V_GetCmp (&uf, face) * faces[face].Aj * V_GetCmp (&xwf, face) / elements[element].Vp;

			}

		      /*
		         // Non-orthogonal correction terms
		         if (parameter.orthof != 0.0)
		         {
		         bpu += parameter.orthof * 
		         viscj * faces[face].Aj / (faces[face].dj + faces[face].kj) /
		         elements[element].Vp * GeoDotVectorVector (gradup,
		         GeoSubVectorVector
		         (faces[face].
		         rpl,
		         elements
		         [element].
		         celement));

		         bpv += parameter.orthof * 
		         viscj * faces[face].Aj / (faces[face].dj + faces[face].kj) /
		         elements[element].Vp * GeoDotVectorVector (gradvp,
		         GeoSubVectorVector
		         (faces[face].
		         rpl,
		         elements
		         [element].
		         celement));

		         bpw += parameter.orthof * 
		         viscj * faces[face].Aj / (faces[face].dj + faces[face].kj) /
		         elements[element].Vp * GeoDotVectorVector (gradwp,
		         GeoSubVectorVector
		         (faces[face].
		         rpl,
		         elements
		         [element].               
		         celement));

		         }
		       */
		    }
		}
	    }
	}

      if (dt > 0)
	{

	  // Unsteady term - Euler

	  app += densp / dt * parameter.st;

	  bpu += densp / dt * V_GetCmp (&xu0l, element);
	  bpv += densp / dt * V_GetCmp (&xv0l, element);
	  bpw += densp / dt * V_GetCmp (&xw0l, element);

	}

      // Source - viscous term
      /*
         gradvisc = Gradient (&viscl, NULL, LOGICAL_FALSE, element);

         bpu += GeoDotVectorVector (gradup, gradvisc);
         bpv += GeoDotVectorVector (gradvp, gradvisc);
         bpw += GeoDotVectorVector (gradwp, gradvisc);
       */

      // Source - gravity
      bpu += densp * g.x;
      bpv += densp * g.y;
      bpw += densp * g.z;

      // Initialize H with source contribution without pressure 
      V_SetCmp (&hu, elements[element].index, bpu);
      V_SetCmp (&hv, elements[element].index, bpv);
      V_SetCmp (&hw, elements[element].index, bpw);

      // Source - pressure 
      /*
         gradp = Gradient (&xpl, &xpf, LOGICAL_TRUE, element);

         bpu += -gradp.x;
         bpv += -gradp.y;
         bpw += -gradp.z;
       */

      if (app == 0.0 || app != app)
	{
	  PetscPrintf (PETSC_COMM_WORLD, "\nError: Problem setting up momentum matrix\n");
	  exit (LOGICAL_ERROR);
	}

      V_SetCmp (&ap, elements[element].index, app);

      /*
      Q_SetEntry (&Am, elements[element].index, elements[element].index, app);

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

      ncols = 0;
      nvals = 0;

      row = elements[element].index;

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

      val[nvals] = app;
      nvals++;

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

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

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

      if (parameter.calc[iu] == LOGICAL_TRUE)
	V_SetCmp (&bu, elements[element].index, bpu);

      if (parameter.calc[iv] == LOGICAL_TRUE)
	V_SetCmp (&bv, elements[element].index, bpv);

      if (parameter.calc[iw] == LOGICAL_TRUE)
	V_SetCmp (&bw, elements[element].index, bpw);

    }

  VecGhostRestoreLocalForm (xu0, &xu0l);
  VecGhostRestoreLocalForm (xv0, &xv0l);
  VecGhostRestoreLocalForm (xw0, &xw0l);

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

  VecGhostRestoreLocalForm (dens, &densl);
  VecGhostRestoreLocalForm (visc, &viscl);

  MatAssemblyBegin (Am, MAT_FINAL_ASSEMBLY);
  MatAssemblyEnd (Am, MAT_FINAL_ASSEMBLY);

  if (parameter.calc[iu] == LOGICAL_TRUE)
    {
      VecAssemblyBegin (bu);
      VecAssemblyEnd (bu);
    }

  if (parameter.calc[iv] == LOGICAL_TRUE)
    {
      VecAssemblyBegin (bv);
      VecAssemblyEnd (bv);
    }

  if (parameter.calc[iw] == LOGICAL_TRUE)
    {
      VecAssemblyBegin (bw);
      VecAssemblyEnd (bw);
    }

  VecAssemblyBegin (ap);
  VecAssemblyEnd (ap);
  VecAssemblyBegin (hu);
  VecAssemblyEnd (hu);
  VecAssemblyBegin (hv);
  VecAssemblyEnd (hv);
  VecAssemblyBegin (hw);
  VecAssemblyEnd (hw);

}

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

  if (parameter.calc[ip] == LOGICAL_TRUE)
    {
  
      // Correct face values      
      CorrectFaceUVW ();

      // Correct cell center   
      CorrectVelocityField ();

  }
  
}

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

  double mres;
  int miter;
  double mtime;
      
  Q_Constr (&Am, nbelements, LOGICAL_FALSE);
  V_Constr (&bu, nbelements, 0);	// Momentum source x-component
  V_Constr (&bv, nbelements, 0);	// Momentum source y-component
  V_Constr (&bw, nbelements, 0);	// Momentum source z-component

  // Store previous time step values
  VecCopy (xu, xu0);
  VecCopy (xv, xv0);
  VecCopy (xw, xw0);

  // Build three momentum matrices for u, v, w velocity components
  if (parameter.calc[ip] == LOGICAL_TRUE)
    {
      BuildMomentumMatrix (dt);
    }

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

  if (parameter.calc[iu] == LOGICAL_TRUE)
    {

      fiter[iu]++;

      // Solve matrix for u velocity component
      SolveMatrix (&Am, &xu, &bu, &miter, &mres, &mtime, parameter.msolver[iu], parameter.mprecond[iu], parameter.miter[iu], parameter.mtol[iu]);

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

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

  if (parameter.calc[iv] == LOGICAL_TRUE)
    {

      fiter[iv]++;

      // Solve matrix for v velocity component
      SolveMatrix (&Am, &xv, &bv, &miter, &mres, &mtime, parameter.msolver[iv], parameter.mprecond[iv], parameter.miter[iv], parameter.mtol[iv]);

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

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

    }

  if (parameter.calc[iw] == LOGICAL_TRUE)
    {

      fiter[iw]++;

      // Solve matrix for w velocity component
      SolveMatrix (&Am, &xw, &bw, &miter, &mres, &mtime, parameter.msolver[iw], parameter.mprecond[iw], parameter.miter[iw], parameter.mtol[iw]);

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

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

  // Calculate correction factors
  if (parameter.calc[ip] == LOGICAL_TRUE)
    {
      CalculateCorrectionFactors ();
    }

  Q_Destr (&Am);

  V_Destr (&bu);
  V_Destr (&bv);
  V_Destr (&bw);
  
}