[alberto]

Hello,

in trying to use setValues to exclude cells were the momentum equation is not defined (rho ~ 0) from the matrix, I found a potential bug/limitation (if I didn't do anything wrong in my code ).

Let's assume the momentum equation is in the compressible form, and I want to exclude cells were rho < SMALL. I loop to identify those cells:

labelHashSet doNotSolveMomentum;

forAll(rho, cellI)
{
  if (rho[cellI] <= 1.0e-6)
  {
    doNotSolveMomentum.insert(cellI);
  }
}

Then I create the list of cell labels and the field of values I want to fix:

labelList noMomentumCellLabels = doNotSolveMomentum.toc();
vectorField fixedU (noMomentumCellLabels.size());

forAll(fixedU, cellI)
{
  fixedU[cellI] = vector(0.0, 0.0, 0.0);
}

UEqn.setValues(noMomentumCellLabels, fixedU);

Now, when I check the central coefficient, it is zero in the cells where the condition is satisfied with the code:

myA = UEqn.A();

forAll(rho,i) 
{

  if (rho[i] < 1.0e-6)
  {
    Info << "rho = " << rho[i] << " A = " << myA[i] << endl;
  }

}

I obtain a list of zeros for A, which causes a division by zero when computing U = H/A.

Is this the expected behaviour of setValues? Shouldn't A be set to 1 in those cells, so that the equation actually degenerates in U = myValue?

Best,

[henry]

setValues currently does not reset the diagonal, it uses the value the matrix provides to set the value.

H

[hjasak]

Hi Alberto,

You cannot just go in and change the diagonal value to 1 because the order of magnitude of other terms in the matrix is unknown: other diagonals can be anything betwees, say 1e-10 and 1e10, depending what you're solving.

If you give the diagonal of a wrong order of magnitude, it will change matrix normalization and that makes the solver behave in unexpected ways. I had this trouble with the FEM solvers a long time ago and the solution was to leave the diagonal as is.

How come you have a zero diagonal in the first place? Can you estimate the size of the diagonal coefficient, eg. cell volume/delta t?

Hope this helps,

Hrv

[alberto]

Thanks to both of you for your replies.

Originally Posted by hjasak

How come you have a zero diagonal in the first place? Can you estimate the size of the diagonal coefficient, eg. cell volume/delta t?

The diagonal is zero because in those cells the equation is undefined (0 = 0), due to rho being zero (or very small).

At the moment I am using a flag to set up the fvMatrix, doing something like

fvVectorMatrix UEqn
(
   solveMomentum*
   (
          fvm::ddt(...).
       + ...
   ) +
   doNotSolveMomentum*
   (
       fvm::Sp(Uunits,U)   
   )
)

but I was looking for a cleaner and possibly less intrusive solution.

Best,

[alberto]

I forgot to answer to one point: in my case both the cell volume and the time step are usually constant (static mesh, fixed dt), so their ratio is known.

Best,