Saya mencoba menerapkan versi 2D makalah Foster dan Fedkiw, "Animasi Praktis Cairan" di sini: http://physbam.stanford.edu/~fedkiw/papers/stanford2001-02.pdf
Sebagian besar semuanya berfungsi, kecuali bagian 8: "Konservasi Misa." Di sana, kami membuat matriks persamaan untuk menghitung tekanan yang diperlukan untuk membuat cairan divergen bebas.
Saya percaya kode saya cocok dengan kertas, namun saya mendapatkan matriks yang tidak dapat diselesaikan selama konservasi langkah massal.
Berikut langkah-langkah saya untuk menghasilkan matriks A:
- Atur entri diagonal ke negatif dari jumlah sel cair yang berdekatan ke sel i.
- Atur entri dan menjadi 1 jika kedua sel i dan j memiliki cairan.
Perhatikan bahwa, dalam implementasi saya, sel , dalam kotak cair sesuai dengan baris gridWidth dalam matriks.
Makalah ini menyebutkan, "Objek statis dan sel kosong tidak mengganggu struktur ini. Dalam hal ini istilah tekanan dan kecepatan dapat menghilang dari kedua sisi", jadi saya menghapus kolom dan baris untuk sel yang tidak memiliki cairan.
Jadi pertanyaan saya adalah: Mengapa matriks saya tunggal? Apakah saya melewatkan semacam kondisi batas di tempat lain di koran? Apakah ini fakta bahwa implementasi saya adalah 2D?
Berikut ini adalah contoh matriks dari implementasi saya untuk kisi 2x2 di mana sel pada 0,0 tidak memiliki cairan:
-1 0 1
0 -1 1
1 1 -2
Edit
Penelitian saya membuat saya percaya bahwa saya tidak benar menangani kondisi batas.
Pertama-tama, pada titik ini saya dapat mengatakan bahwa matriks saya mewakili persamaan Poisson tekanan diskrit. Ini adalah analog diskrit dari penerapan operator Laplacian yang menggabungkan perubahan tekanan lokal ke divergensi sel.
Sejauh yang saya bisa mengerti, karena kita berurusan dengan perbedaan tekanan, kondisi batas diperlukan untuk "menjangkar" tekanan ke nilai referensi absolut. Kalau tidak, mungkin ada sejumlah solusi tak terbatas untuk set persamaan.
Dalam catatan ini , 3 cara berbeda diberikan untuk menerapkan kondisi batas, sejauh yang saya pahami:
Dirichlet - menentukan nilai absolut di batas.
Neummann - menentukan turunan pada batas.
Robin - menentukan beberapa jenis kombinasi linear dari nilai absolut dan turunannya pada batas.
Makalah Foster dan Fedki tidak menyebutkan semua ini, tetapi saya percaya bahwa mereka menegakkan kondisi batas Dirichlet, penting karena pernyataan ini pada akhir 7.1.2, "Tekanan dalam sel permukaan diatur ke tekanan atmosfer."
Saya telah membaca catatan yang saya tautkan beberapa kali dan masih tidak mengerti perhitungan matematika yang terjadi. Bagaimana tepatnya kita menegakkan syarat batas ini? Melihat implementasi lain, tampaknya ada semacam gagasan tentang sel "Hantu" yang terletak pada batas.
Di sini saya telah menautkan ke beberapa sumber yang mungkin bermanfaat bagi orang lain yang membaca ini.
Catatan tentang kondisi batas untuk Matriks Poisson
Postingan Ilmu Komputasi StackExchange pada kondisi batas Neumann
Pos Ilmu Komputasi StackExchange di Poisson Solver
Berikut adalah kode yang saya gunakan untuk menghasilkan matriks. Perhatikan bahwa, alih-alih menghapus kolom dan baris secara eksplisit, saya membuat dan menggunakan peta dari indeks sel cair ke kolom / baris matriks terakhir.
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
}
}
}
sumber
Jawaban:
Dari cuplikan kode Anda dan hasil Anda untuk contoh 2x2, saya dapat melihat bahwa Anda benar-benar mensimulasikan domain dengan hanya syarat batas Neumann (dinding selip). Dalam hal ini, sistem berisi ruang kosong dan matriks Anda tunggal.
Jika ini adalah konfigurasi simulasi yang Anda inginkan (yaitu tanpa Dirichlet (tekanan) BC), Anda perlu memproyeksikan ruang nol dari solusi Anda. Ini mudah jika Anda menggunakan gradien konjugat (CG) seperti yang disarankan dalam makalah itu. Di setiap iterasi dari iterasi CG Anda, cukup ambil vektor solusi saat ini , dan lakukan mana adalah ruang null dinormalisasi dari operator gradien: , .x x′=(I−u^u^T)x=x−(u^⋅x)u^ u^ u=(1,1,…,1) u^=u∥u∥
Kalau tidak, jika Anda ingin mensimulasikan udara (batas bebas atau Dirichlet BC), Anda perlu membedakan dinding dan sel udara (yaitu memiliki boolean
hasLiquid
tidak cukup), dan menerapkan diskritisasi yang benar untuk mereka (lihat di bawah).Sebagai catatan akhir, entri diagonal Anda negatif. Anda mungkin ingin membalik tanda sehingga metode CG bekerja.
Di bawah ini saya ingin menunjukkan lebih banyak detail. Pertimbangkan proses proyeksi tekanan. Nyatakan kecepatan sebelum proyeksi tekanan sebagai . Itu bisa saja berbeda, jadi kami menghitung tekanan untuk memperbaikinya dan mendapatkan kecepatan bebas divergensi . Yaitu, Ambil divergence darinya, dan karena bebas dari perbedaan, Misalkan tidak ada tekanan Dirichlet BC presesnt dan kami memiliki satu solusi , maka untuk konstantav∗ vn+1 vn+1=v∗−Δtρ∇P vn+1 ∇⋅∇P=∇⋅v∗ P0 P0+c c juga merupakan solusi karena .
adalah ruang nol yang ingin kita proyeksikan.∇⋅∇(P0+c)=∇⋅∇P0=∇⋅v∗ c
Untuk menangani Dirichlet BC, mari kita pertimbangkan kasus 1D sebagai contoh. Asumsikan Anda menggunakan grid terhuyung-huyung, di mana tekanan berada di pusat-pusat grid dan kecepatan terletak di wajah antara node dan . Maka diskritisasi umum untuk satu sel adalah Misalkan adalah sel udara, yaitu tekanannya telah ditentukan, maka istilah dipindahkan ke sisi kanan dan menghilang dari matriks. Perhatikan bahwa jumlah istilah diagonal masih dua. Itu sebabnya saya mengatakan contoh 2x2 Anda tidak mengandung Dirichlet BC.pi vi+1/2 i i+1 pi+1−pi−(pi−pi−1)Δx2=rhs pi+1 pi+1 pi
Dengan Dirichlet atau Neumann BC, matriks selalu pasti positif simetris. Itu sebabnya penulis berkata
sumber