Conservando la masa en la simulación de líquidos

8

Estoy tratando de implementar una versión 2D del documento de Foster y Fedkiw, "Animación práctica de líquidos" aquí: http://physbam.stanford.edu/~fedkiw/papers/stanford2001-02.pdf

Casi todo funciona, excepto la sección 8: "Conservación de la misa". Allí, configuramos una matriz de ecuaciones para calcular las presiones necesarias para liberar el líquido divergente.

Creo que mi código coincide con el documento, sin embargo, obtengo una matriz insoluble durante la conservación del paso de masa.

Aquí están mis pasos para generar la matriz A:

  1. Establezca las entradas diagonales en el negativo del número de celdas líquidas adyacentes a la celda i.UNAyo,yo
  2. Establecer las entradas UNAyo,j y UNAj,yo a 1 si ambas celdas i y j tienen líquido.

Tenga en cuenta que, en mi implementación, la celda yo,j en la rejilla líquida corresponde a la fila yo+gridWidthj en la matriz

El documento menciona: "El objeto estático y las celdas vacías no interrumpen esta estructura. En ese caso, los términos de presión y velocidad pueden desaparecer de ambos lados", por lo que elimino las columnas y filas de las celdas que no tienen líquido.

Entonces mi pregunta es: ¿Por qué mi matriz es singular? ¿Me estoy perdiendo algún tipo de condición límite en algún otro lugar del periódico? ¿Es el hecho de que mi implementación es 2D?

Aquí hay una matriz de ejemplo de mi implementación para una cuadrícula de 2x2 donde la celda a 0,0 no tiene líquido:

-1   0   1

 0  -1   1

 1   1  -2

Editar

Mi investigación me ha llevado a creer que no estoy manejando adecuadamente las condiciones de contorno.

En primer lugar, en este punto puedo decir que mi matriz representa la ecuación de Poisson de presión discreta. Es el análogo discreto de aplicar el operador laplaciano acoplando los cambios de presión locales a la divergencia celular.

Hasta donde puedo entender, dado que estamos lidiando con diferencias de presión, se necesitan condiciones de frontera para "anclar" las presiones a un valor de referencia absoluto. De lo contrario, puede haber un número infinito de soluciones para el conjunto de ecuaciones.

En estas notas , se dan 3 formas diferentes de aplicar condiciones de contorno, a mi entender:

  1. Dirichlet: especifica valores absolutos en los límites.

  2. Neummann: especifica la derivada en los límites.

  3. Robin: especifica algún tipo de combinación lineal del valor absoluto y la derivada en los límites.

El artículo de Foster y Fedki no menciona ninguno de estos, pero creo que imponen las condiciones de frontera de Dirichlet, notables debido a esta declaración al final de 7.1.2, "La presión en una celda de superficie se ajusta a la presión atmosférica".

He leído las notas que vinculé varias veces y todavía no entiendo las matemáticas. ¿Cómo aplicamos exactamente estas condiciones límite? Mirando otras implementaciones, parece haber algún tipo de noción de celdas "fantasmas" que se encuentran en el límite.

Aquí he vinculado a algunas fuentes que pueden ser útiles para otros que leen esto.

Notas sobre condiciones de contorno para matrices de Poisson

Computational Science StackExchange post en condiciones de contorno de Neumann

Computational Science StackExchange post en Poisson Solver

Implementación de agua Physbam


Aquí está el código que uso para generar la matriz. Tenga en cuenta que, en lugar de eliminar explícitamente columnas y filas, genero y uso un mapa de índices de celdas líquidas a las columnas / filas de la matriz final.

for (int i = 0; i < cells.length; i++) {
  for (int j = 0; j < cells[i].length; j++) {
    FluidGridCell cell = cells[i][j];

    if (!cell.hasLiquid)
      continue;

    // get indices for the grid and matrix
    int gridIndex = i + cells.length * j;
    int matrixIndex = gridIndexToMatrixIndex.get((Integer)gridIndex);

    // count the number of adjacent liquid cells
    int adjacentLiquidCellCount = 0;
    if (i != 0) {
      if (cells[i-1][j].hasLiquid)
        adjacentLiquidCellCount++;
    }
    if (i != cells.length-1) {
      if (cells[i+1][j].hasLiquid)
        adjacentLiquidCellCount++;
    }
    if (j != 0) {
      if (cells[i][j-1].hasLiquid)
      adjacentLiquidCellCount++;
    }
    if (j != cells[0].length-1) {
      if (cells[i][j+1].hasLiquid)
        adjacentLiquidCellCount++;
    }

    // the diagonal entries are the negative count of liquid cells
    liquidMatrix.setEntry(matrixIndex, // column
                          matrixIndex, // row
                          -adjacentLiquidCellCount); // value

    // set off-diagonal values of the pressure matrix
    if (cell.hasLiquid) {
      if (i != 0) {
        if (cells[i-1][j].hasLiquid) {
          int adjacentGridIndex = (i-1) + j * cells.length;
          int adjacentMatrixIndex = gridIndexToMatrixIndex.get((Integer)adjacentGridIndex);
          liquidMatrix.setEntry(matrixIndex, // column
                                adjacentMatrixIndex, // row
                                1.0); // value
          liquidMatrix.setEntry(adjacentMatrixIndex, // column
                                matrixIndex, // row
                                1.0); // value
        }
      }
      if (i != cells.length-1) {
        if (cells[i+1][j].hasLiquid) {
          int adjacentGridIndex = (i+1) + j * cells.length;
          int adjacentMatrixIndex = gridIndexToMatrixIndex.get((Integer)adjacentGridIndex);
          liquidMatrix.setEntry(matrixIndex, // column
                                adjacentMatrixIndex, // row
                                1.0); // value
          liquidMatrix.setEntry(adjacentMatrixIndex, // column
                                matrixIndex, // row
                                1.0); // value
        }
      }
      if (j != 0) {
        if (cells[i][j-1].hasLiquid) {
          int adjacentGridIndex = i + (j-1) * cells.length;
          int adjacentMatrixIndex = gridIndexToMatrixIndex.get((Integer)adjacentGridIndex);
          liquidMatrix.setEntry(matrixIndex, // column
                                adjacentMatrixIndex, // row
                                1.0); // value
          liquidMatrix.setEntry(adjacentMatrixIndex, // column
                                matrixIndex, // row
                                1.0); // value
        }
      }
      if (j != cells[0].length-1) {
        if (cells[i][j+1].hasLiquid) {
          int adjacentGridIndex = i + (j+1) * cells.length;
          int adjacentMatrixIndex = gridIndexToMatrixIndex.get((Integer)adjacentGridIndex);
          liquidMatrix.setEntry(matrixIndex, // column
                                adjacentMatrixIndex, // row
                                1.0); // value
          liquidMatrix.setEntry(adjacentMatrixIndex, // column
                                matrixIndex, // row
                                1.0); // value
        }
      }
    }
Jared cuenta
fuente
No está claro por qué los objetos estáticos y las celdas vacías deberían permitir la eliminación de filas y columnas. ¿Está configurando estas filas y columnas a cero o eliminándolas por completo para obtener una matriz más pequeña?
trichoplax
En caso de que el problema esté en otro lugar que no sea donde adivina, sería útil ver el código, si es algo que está feliz de compartir. Idealmente un MCVE
trichoplax
Hola trichoplax. Una matriz con una fila o columna completamente cero sería singular, hasta donde yo sé, por lo que, en cambio, los elimino de la matriz para hacer una matriz más pequeña (así como sus entradas correspondientes en el vector b).
Jared cuenta el
Editaré un MCVE esta noche cuando esté cerca de mi computadora con la fuente.
Jared cuenta el
También sospeché que tal vez estaba haciendo una suposición incorrecta en otro lugar del código, sin embargo, esto solo se refiere a la estructura de la matriz (y si es singular o no). Lo único que se me ocurre es lo que califica como una "celda de superficie" frente a una celda de aire o una celda de líquido. Si esta es una celda de líquido adyacente a una celda de aire, ¿hay algo diferente que debería hacer con sus columnas / filas correspondientes?
Jared cuenta el

Respuestas:

2

A partir de su fragmento de código y su resultado para el ejemplo 2x2, puedo ver que en realidad está simulando un dominio con solo condiciones de contorno de Neumann (pared de deslizamiento). En este caso, el sistema contiene un espacio nulo y su matriz es singular.

Si esta es la configuración de simulación que desea (es decir, sin Dirichlet (presión) BC), deberá proyectar el espacio nulo de su solución. Esto es sencillo si está utilizando el gradiente conjugado (CG) como se sugiere en ese documento. En cada iteración de su iteración CG, simplemente tome el vector de solución actualX, y hacer

X=(yo-tu^tu^T)X=X-(tu^X)tu^
dónde tu^ es el espacio nulo normalizado del operador de gradiente: tu=(1,1,...,1), tu^=tutu.

De lo contrario, si desea simular el aire (límite libre o Dirichlet BC), deberá distinguir un muro y una celda de aire (es decir, tener un booleano hasLiquidno es suficiente) y aplicar la discretización correcta para ellos (ver más abajo).

Como nota final, sus entradas diagonales son negativas. Es posible que desee voltear los signos para que funcione el método CG.


A continuación me gustaría mostrar más detalles. Considere el proceso de proyección de presión. Denote la velocidad antes de la proyección de presión comov. Podría ser divergente, por lo que calculamos la presión para corregirlo y obtener la velocidad libre de divergenciavnorte+1. Es decir,

vnorte+1=v-ΔtρPAGS
Tome la divergencia de eso, y desde vnorte+1 es divergencia libre,
PAGS=v
Supongamos que no hay presión Dirichlet BC presente y tenemos una solución PAGS0 0, entonces PAGS0 0+C para cualquier constante C también es una solución porque (PAGS0 0+C)=PAGS0 0=v. C es el espacio nulo que queremos proyectar.

Para manejar el Dirichlet BC, consideremos el caso 1D como ejemplo. Suponga que está utilizando una cuadrícula escalonada, donde las presionespagsyo están ubicados en centros de rejilla y velocidades vyo+1/ /2 están ubicados en las caras entre nodos yo y yo+1. Entonces la discretización general para una celda es

pagsyo+1-pagsyo-(pagsyo-pagsyo-1)ΔX2=rhs
Suponer que pagsyo+1 es una celda de aire, es decir, se ha especificado su presión, entonces el pagsyo+1El término se mueve hacia el lado derecho y desaparece de la matriz. Tenga en cuenta que el recuento del término diagonalpagsyoTodavía son dos. Es por eso que dije que su ejemplo de 2x2 no contenía un Dirichlet BC.

Con Dirichlet o Neumann BC, la matriz siempre es simétrica positiva definida. Por eso los autores dijeron

Static object and empty cells don’t disrupt this structure.
In that case pressure and velocity terms can disappear from both sides
TheBusyTypist
fuente