GPU-Accelerated Poisson Solver for Complex Geometries Karmaşık Geometrilerde GPU-Hızlandırılmış Poisson Çözücü

The Poisson Equation in Physics and Engineering Fizik ve Mühendislikte Poisson Denklemi

The Poisson equation $\nabla^2\phi = f$ is one of the most fundamental partial differential equations in physics. It links a field's curvature (how the potential bends in space) to the source that generates it. Whenever we ask "given some source distribution, what does the resulting field look like everywhere?", we are solving a Poisson problem. Poisson denklemi $\nabla^2\phi = f$, fiziğin en temel kısmi diferansiyel denklemlerinden biridir. Bir alanın eğriliğini (potansiyelin uzayda nasıl büküldüğünü) onu üreten kaynağa bağlar. "Verilen bir kaynak dağılımı için ortaya çıkan alan her yerde nasıl görünür?" sorusunu her sorduğumuzda, bir Poisson problemi çözüyoruzdur.

Where Does It Appear? Nerelerde Karşımıza Çıkar?

Why Fast Solving Matters Hızlı Çözüm Neden Kritik

In CFD applications, the Poisson solve is performed at every single time step, often consuming 60–80% of total computation time. For a transient 3D simulation with $10^5$ time steps on a $1024^3$ grid, the Poisson solver is called hundreds of thousands of times. A 2× improvement in Poisson solve speed can translate to weeks of saved wall-clock time on cluster resources. This is precisely why GPU-accelerated solvers are not a luxury — they are a practical necessity for production-scale simulations. CFD uygulamalarında, Poisson çözümü her bir zaman adımında yapılır ve toplam hesaplama süresinin genellikle %60–80'ini tüketir. $1024^3$ grid üzerinde $10^5$ zaman adımlı bir geçişli 3B simülasyon için Poisson çözücü yüz binlerce kez çağrılır. Poisson çözüm hızındaki 2× iyileşme, küme kaynaklarında haftalarca duvar saati süresinden tasarruf anlamına gelebilir. GPU-hızlandırılmış çözücülerin bir lüks değil, üretim ölçeğinde simülasyonlar için pratik bir gereklilik olmasının sebebi budur.

Challenges in the Literature Literatürdeki Problemler

Building a high-performance Poisson solver for realistic geometries is far from straightforward. Four key challenges recur across the literature: Gerçekçi geometriler için yüksek performanslı bir Poisson çözücü oluşturmak hiç de kolay değildir. Literatürde dört temel problem sürekli karşımıza çıkar:

Our project addresses all four: we use binary masking to handle arbitrary geometries (§3), asynchronous halo exchange to hide communication latency (§5), multi-grid and neural-network warm starts for fast convergence (§6), and careful boundary treatment at domain contours. Projemiz dördünü de ele alır: keyfi geometrileri işlemek için binary maskeleme (§3), iletişim gecikmesini gizlemek için asenkron halo değişimi (§5), hızlı yakınsama için multi-grid ve sinir ağı sıcak başlangıçları (§6) ve domain konturlarında dikkatli sınır işleme kullanırız.

From Maxwell to Poisson Maxwell'den Poisson'a

The Poisson equation is not an arbitrary mathematical construction — it emerges directly from the fundamental laws of electromagnetism. In the electrostatic limit, two of Maxwell's equations combine to produce it. Poisson denklemi keyfi bir matematiksel yapı değildir — doğrudan elektromanyetizmanın temel yasalarından çıkar. Elektrostatik limitte, Maxwell denklemlerinden ikisi birleşerek onu üretir.

Here $\phi(\mathbf{r})$ is the electrostatic potential and $\rho(\mathbf{r})$ is the charge density. When the charge distribution is known and boundary conditions are specified, finding $\phi$ everywhere in the domain is the core computational problem. For grounded conductors with complex shapes — like aircraft fuselage cross-sections — we apply Dirichlet boundary conditions ($\phi = 0$ on the conductor surface). Burada $\phi(\mathbf{r})$ elektrostatik potansiyel, $\rho(\mathbf{r})$ ise yük yoğunluğudur. Yük dağılımı biliniyor ve sınır koşulları belirlendiğinde, $\phi$'yi domainin her yerinde bulmak temel hesaplama problemidir. Karmaşık şekillere sahip topraklanmış iletkenler için — uçak gövde kesitleri gibi — Dirichlet sınır koşulları uygulanır (iletken yüzeyinde $\phi = 0$).

Finite Difference Method Sonlu Farklar Yöntemi

To solve the Poisson equation numerically, we discretize the continuous domain onto a uniform Cartesian grid with spacing $h$. The Laplacian operator $\nabla^2$ is approximated using the finite difference method (FDM), which we derive from the Taylor expansion. Poisson denklemini sayısal olarak çözmek için, sürekli domaini $h$ adım büyüklüğüne sahip düzenli bir Kartezyen ızgara üzerine diskritize ederiz. Laplasyen operatörü $\nabla^2$, Taylor açılımından türettiğimiz sonlu farklar yöntemi (FDM) ile yaklaşıklanır.

Taylor Series Derivation of the Second Derivative İkinci Türevin Taylor Serisi Türetmesi

The 5-Point Stencil 5 Noktalı Stencil

Applying the central difference formula (Eq. 2) to both $x$ and $y$ directions on a 2D grid where $(x_i, y_j) = (ih, jh)$: Merkezi fark formülünü (Denk. 2) $(x_i, y_j) = (ih, jh)$ olan 2B ızgara üzerinde hem $x$ hem $y$ yönlerinde uygularsak:

The stencil pattern shows how each grid point depends on its four nearest neighbors. This local dependency is what makes the method highly parallelizable on GPU architectures — but it also introduces the halo problem when the grid is split across multiple processors. Stencil paterni, her ızgara noktasının en yakın dört komşusuna nasıl bağlı olduğunu gösterir. Bu yerel bağımlılık, yöntemi GPU mimarilerinde yüksek oranda paralelleştirilebilir kılar — ancak grid birden fazla işlemciye bölündüğünde halo problemini de ortaya çıkarır.

This discretization transforms the continuous PDE into a large sparse linear system $A\mathbf{x} = \mathbf{b}$, where $A$ is a pentadiagonal matrix, $\mathbf{x}$ contains the unknown potential values at each grid point, and $\mathbf{b}$ encodes the source distribution and boundary conditions. For a $500 \times 500$ grid, this means a system of $250{,}000$ unknowns. For $2048 \times 2048$, it grows to over 4 million — well beyond direct solvers, motivating the iterative approaches in §4. Bu diskritizasyon, sürekli PDE'yi büyük bir seyrek lineer sisteme $A\mathbf{x} = \mathbf{b}$ dönüştürür; burada $A$ pentadiyagonal bir matristir, $\mathbf{x}$ her ızgara noktasındaki bilinmeyen potansiyel değerlerini içerir ve $\mathbf{b}$ kaynak dağılımını ve sınır koşullarını kodlar. $500 \times 500$'lük bir ızgara için bu, $250{,}000$ bilinmeyenli bir sistem demektir. $2048 \times 2048$ için 4 milyonun üzerine çıkar — doğrudan çözücülerin çok ötesinde, §4'teki iteratif yaklaşımları motive eder.

Domain Representation with Binary Masking Binary Maskeleme ile Domain Temsili

Real-world physical systems rarely possess simple rectangular shapes. Our approach uses a binary mask matrix $M(i,j)$ on the Cartesian grid: Gerçek dünya fiziksel sistemleri nadiren basit dikdörtgen şekillere sahiptir. Yaklaşımımız, Kartezyen ızgara üzerinde bir ikili maske matrisi $M(i,j)$ kullanır:

The boundary of the domain — where $M$ transitions from 1 to 0 — is where we enforce Dirichlet conditions ($\phi = 0$ on the contour). Any geometry representable as a black-and-white image can serve as input. We demonstrated this by using a 2D cross-section of the TUSAŞ KAAN fighter aircraft as the computational domain. Domainin sınırı — $M$'nin 1'den 0'a geçtiği yer — Dirichlet koşullarını uyguladığımız yerdir (kontür üzerinde $\phi = 0$). Siyah-beyaz görsel olarak temsil edilebilen herhangi bir geometri girdi olarak kullanılabilir. Bunu TUSAŞ KAAN savaş uçağının 2B kesitini hesaplama domaini olarak kullanarak gösterdik.

Image-to-Mask Pipeline Görselden Maskeye Dönüşüm Süreci

The key innovation is that any geometry can be imported as a standard image file. The conversion from a photograph or CAD rendering to a computable binary mask follows a four-step pipeline: Temel yenilik, herhangi bir geometrinin standart bir görüntü dosyası olarak içe aktarılabilmesidir. Bir fotoğraftan veya CAD görselinden hesaplanabilir bir binary maskeye dönüşüm dört adımlı bir pipeline izler:

Test Geometry: Heart-Shaped Domain Test Geometrisi: Kalp Şekilli Domain

For systematic benchmarking, we used an algebraic heart curve: Sistematik karşılaştırma testleri için cebirsel bir kalp eğrisi kullandık:

Mask ImplementationMask Uygulaması

int isInsideHeart(double x, double y) {
    double expr = (x * x + y * y - 1.0);
    double val = (expr * expr * expr) - (x * x * y * y * y);
    return (val <= 0.0) ? 1 : 0;
}

// build the mask matrix over the Cartesian grid
for (int i = 0; i < GRID_SIZE; i++) {
    for (int j = 0; j < GRID_SIZE; j++) {
        int idx = i * GRID_SIZE + j;
        mask[idx] = isInsideHeart(xVals[i], yVals[j]);
    }
}

Why Iterative Methods? Neden İteratif Yöntemler?

For a $2048\times2048$ grid, the linear system has over 4 million unknowns. Direct methods (LU decomposition, Cholesky) would require $\mathcal{O}(N^3)$ operations and $\mathcal{O}(N^2)$ memory — completely infeasible. Iterative methods start from an initial guess and progressively refine it, exploiting the sparse structure of the matrix. The question is: which iterative method gives the best balance of convergence speed, parallelizability, and implementation complexity? $2048\times2048$ grid için lineer sistemde 4 milyonun üzerinde bilinmeyen vardır. Doğrudan yöntemler (LU ayrıştırma, Cholesky) $\mathcal{O}(N^3)$ işlem ve $\mathcal{O}(N^2)$ bellek gerektirir — tamamen uygulanabilir değildir. İteratif yöntemler bir başlangıç tahmininden başlayıp matrisin seyrek yapısından faydalanarak onu aşamalı olarak iyileştirir. Soru şudur: hangi iteratif yöntem yakınsama hızı, paralelleştirilebilirlik ve uygulama karmaşıklığı arasında en iyi dengeyi verir?

Ahmet Ali Akkurt implemented all six methods from scratch in C/C++, ran systematic benchmarks on the heart domain, and characterized their convergence properties. This section represents his primary research contribution to the project. Ahmet Ali Akkurt altı yöntemin tamamını C/C++ ile sıfırdan kodladı, kalp domaininde sistematik benchmark'lar çalıştırdı ve yakınsama özelliklerini karakterize etti. Bu bölüm onun projeye temel araştırma katkısını temsil eder.

1. Jacobi Iteration1. Jacobi İterasyonu

2. Gauss-Seidel2. Gauss-Seidel

3. Red-Black Gauss-Seidel3. Red-Black Gauss-Seidel

4. Averaged Red-Black Gauss-Seidel (ARBGS)4. Ortalamalı Red-Black Gauss-Seidel (ARBGS)

5. Conjugate Gradient (CG)5. Eşlenik Gradyan (CG)

6. Preconditioned Conjugate Gradient (PCG)6. Ön Koşullu Eşlenik Gradyan (PCG)

Benchmark ComparisonBenchmark Karşılaştırması

Ahmet ran all six methods on a 500×500 heart domain with $\rho(x,y) = \sin(x)\cos(y)$ and tolerance $10^{-10}$. The table summarizes both convergence characteristics and GPU suitability: Ahmet altı yöntemin tamamını 500×500 kalp domaininde $\rho(x,y) = \sin(x)\cos(y)$ ve $10^{-10}$ tolerans ile çalıştırdı. Tablo hem yakınsama özelliklerini hem de GPU uygunluğunu özetler:

Why We Chose Jacobi Despite Slower Convergence Daha Yavaş Yakınsamaya Rağmen Neden Jacobi Seçtik

The benchmark results above make the case for CG/PCG look overwhelming — 180 iterations versus 180,000. So why would anyone pick Jacobi? The answer lies in the project's primary goal: multi-GPU scalability as a first priority. Yukarıdaki benchmark sonuçları CG/PCG lehine baskın görünür — 180 iterasyona karşı 180,000. Peki neden Jacobi seçilir? Cevap projenin birincil hedefinde yatar: ilk öncelik olarak çoklu GPU ölçeklenebilirliği.

Jacobi's update rule has zero data dependency within each iteration. Every thread independently reads from the previous iteration's buffer and writes to a new buffer. There are no race conditions, no sequential sweeps, no global reductions during the stencil compute. This translates directly to: Jacobi'nin güncelleme kuralının her iterasyon içinde sıfır veri bağımlılığı vardır. Her thread önceki iterasyonun tamponundan bağımsız olarak okur ve yeni bir tampona yazar. Yarış koşulu yok, sıralı tarama yok, stencil hesaplaması sırasında global indirgeme yok. Bu doğrudan şuna dönüşür:

CG and PCG converge orders of magnitude faster, but they require global MPI_Allreduce calls for dot products at every iteration and sequential preconditioner solves. On a 2-GPU system connected via Ethernet, these synchronization points become the dominant cost. We compensate for Jacobi's slow convergence with multi-grid warm starts (§7) and ML-assisted initial guesses (§7), bringing the effective iteration count down by 79%. CG ve PCG büyüklük mertebeleri daha hızlı yakınsar, ancak her iterasyonda iç çarpımlar için global MPI_Allreduce çağrıları ve sıralı preconditioner çözümleri gerektirir. Ethernet üzerinden bağlanan 2 GPU'lu bir sistemde bu senkronizasyon noktaları baskın maliyet olur. Jacobi'nin yavaş yakınsamasını multi-grid sıcak başlangıçları (§7) ve ML destekli başlangıç tahminleri (§7) ile telafi ederek efektif iterasyon sayısını %79 düşürürüz.

Kernel Development Path Kernel Geliştirme Yolu

The CUDA implementation evolved through three stages, each addressing a specific performance bottleneck: CUDA uygulaması her biri belirli bir performans darboğazını ele alan üç aşamadan geçerek evrildi:

GPU architecture, system design, and CUDA implementation represent Berkay Yılmaz's primary contribution.

GPU mimarisi, sistem tasarımı ve CUDA uygulaması Berkay Yılmaz'ın projeye temel katkısını temsil eder.

Halo Exchange Between GPUs GPU'lar Arası Halo Değişimi

In the distributed setting, each GPU owns a horizontal strip of the global grid. Jacobi's stencil needs neighbors — including rows owned by adjacent GPUs. We solve this with ghost layers: each local domain is padded by one extra row on top and bottom. The halo exchange protocol ensures these ghost rows contain up-to-date values from the neighboring GPU: Dağıtık sistemde her GPU, global gridin yatay bir dilimini sahiplenir. Jacobi stencil'i komşulara ihtiyaç duyar — komşu GPU'ların sahip olduğu satırlar dahil. Bunu ghost layer'larla çözeriz: her lokal domain üst ve altta birer ekstra satırla genişletilir. Halo değişim protokolü, bu ghost satırların komşu GPU'dan güncel değerler içermesini sağlar:

Cascadic Multi-Grid Strategy Kademeli Multi-Grid Stratejisi

Starting Jacobi from $\phi^{(0)} = 0$ on a $2048^2$ grid wastes thousands of iterations converging gross features that could be captured cheaply on a coarser grid. Our cascadic approach solves hierarchically: $2048^2$ grid üzerinde Jacobi'yi $\phi^{(0)} = 0$'dan başlatmak, kaba bir gridde ucuza yakalanabilecek genel özellikleri yakınsamak için binlerce iterasyon harcar. Kademeli yaklaşımımız hiyerarşik olarak çözer:

At each interpolation, domain boundaries and Dirichlet conditions ($\phi = 0$ on the contour) are preserved with boundary-aware algorithms. The CUDA implementation uses texture memory for hardware-accelerated bilinear filtering. Result: 79% iteration reduction on the target layer. Her interpolasyonda, domain sınırları ve Dirichlet koşulları ($\phi = 0$ kontür üzerinde) sınır-farkında algoritmalarla korunur. CUDA uygulaması donanım hızlandırmalı bilineer filtreleme için texture memory kullanır. Sonuç: hedef katmanda %79 iterasyon azalması.

ML-Assisted Initial Guess: The Warm-Start Idea ML Destekli Başlangıç Tahmini: Sıcak Başlangıç Fikri

Multi-grid provides a smooth initial guess, but it still starts each level from scratch. The fundamental insight behind our ML approach is this: the mapping from source distribution $\rho$ to solution $\phi$ is a well-defined operator, and a neural network can learn to approximate it. Multi-grid yumuşak bir başlangıç tahmini sağlar, ama her seviyeye yine sıfırdan başlar. ML yaklaşımımızın temel kavrayışı şudur: kaynak dağılımı $\rho$'dan çözüm $\phi$'ye eşleme iyi tanımlı bir operatördür ve bir sinir ağı bunu yaklaşıklamayı öğrenebilir.

Mathematical Motivation Matematiksel Motivasyon

Consider the iterative solver as minimizing the residual $r^{(k)} = b - A\phi^{(k)}$. Starting from $\phi^{(0)} = 0$, the initial residual norm $\|r^{(0)}\| = \|b\|$ is maximal. If we can instead start from $\phi^{(0)} = \hat{\phi}_\text{NN}$ where $\hat{\phi}_\text{NN} \approx \phi^*$ (the true solution), the initial residual drops dramatically: İteratif çözücüyü artık $r^{(k)} = b - A\phi^{(k)}$'yı minimize etme olarak düşünelim. $\phi^{(0)} = 0$'dan başlayarak ilk artık normu $\|r^{(0)}\| = \|b\|$ maksimumdur. Bunun yerine $\hat{\phi}_\text{NN} \approx \phi^*$ (gerçek çözüm) olmak üzere $\phi^{(0)} = \hat{\phi}_\text{NN}$'den başlayabilirsek, ilk artık dramatik olarak düşer:

The network doesn't need to produce an exact solution — it only needs to capture the low-frequency structure of $\phi$. The iterative solver handles the remaining high-frequency corrections efficiently, since Jacobi is actually quite good at smoothing high-frequency errors (spectral radius near zero for short-wavelength modes). Ağın kesin bir çözüm üretmesi gerekmez — yalnızca $\phi$'nin düşük frekanslı yapısını yakalaması yeter. İteratif çözücü kalan yüksek frekanslı düzeltmeleri verimli bir şekilde ele alır, çünkü Jacobi aslında yüksek frekanslı hataları yumuşatmakta oldukça iyidir (kısa dalga boylu modlar için spektral yarıçap sıfıra yakın).

Network Architecture Ağ Mimarisi

We use a fully-connected feedforward network that takes the flattened source distribution $\rho$ (plus the binary mask $M$) as input and outputs the flattened potential field $\hat{\phi}$. The architecture is intentionally kept shallow and lightweight — inference must be faster than the iterations it saves. Düzleştirilmiş kaynak dağılımı $\rho$'yu (artı binary mask $M$'yi) girdi olarak alan ve düzleştirilmiş potansiyel alanı $\hat{\phi}$'yi çıkaran tam bağlantılı ileri beslemeli bir ağ kullanırız. Mimari kasıtlı olarak sığ ve hafif tutulur — çıkarım, tasarruf ettiği iterasyonlardan daha hızlı olmalıdır.

Training & Design Rationale Eğitim ve Tasarım Gerekçesi

The network optimizes the ratio of iteration savings to inference cost. Its ~50K parameters capture the first 10–15 Fourier modes (~85% of solution energy); Jacobi corrects the high-frequency remainder. U-Net, PINN, and FNO alternatives offer higher accuracy but 10–40× higher latency — planned for future work. Ağ, iterasyon tasarrufunun çıkarım maliyetine oranını optimize eder. ~50K parametresi ilk 10–15 Fourier modunu (~%85 çözüm enerjisi) yakalar; Jacobi yüksek frekanslı kalanı düzeltir. U-Net, PINN ve FNO alternatifleri daha yüksek doğruluk sunar ancak 10–40× daha yüksek gecikmeyle — gelecek çalışmalar için planlanmıştır.

Dynamic Adaptive Convergence Control Dinamik Adaptif Yakınsama Kontrolü

Standard distributed solvers check global convergence every iteration via MPI_Allreduce, creating a communication bottleneck. Our algorithm adjusts the check interval based on problem size, reducing communication overhead by 90%: Standart dağıtık çözücüler her iterasyonda MPI_Allreduce ile global yakınsamayı kontrol eder ve iletişim darboğazı yaratır. Algoritmamız kontrol aralığını problem boyutuna göre ayarlayarak iletişim yükünü %90 azaltır:

Interactive Solver Interface İnteraktif Çözücü Arayüzü

Beyond the command-line solver, we developed a desktop GUI (Tkinter/Python frontend + CUDA backend) that enables interactive exploration of the Poisson solver. The interface lets users draw arbitrary domains with a brush tool, import images as geometry masks, select source functions from a library, and visualize results in real-time with 2D heatmaps and 3D surface plots. Komut satırı çözücünün ötesinde, Poisson çözücünün interaktif keşfini sağlayan bir masaüstü GUI (Tkinter/Python arayüzü + CUDA backend) geliştirdik. Arayüz kullanıcıların fırça aracıyla keyfi domainler çizmesine, görüntüleri geometri maskesi olarak içe aktarmasına, bir kütüphaneden kaynak fonksiyonları seçmesine ve sonuçları 2B ısı haritaları ve 3B yüzey grafikleri ile gerçek zamanlı görselleştirmesine olanak tanır.

The GUI serves as both a research tool — enabling rapid experimentation with different geometries and charge distributions — and a demonstration platform for the TUSAŞ LiftUp program presentations. Users can interactively adjust grid resolution (128 to 2048), toggle multi-grid acceleration, and watch the solution converge live via an embedded log console. GUI hem bir araştırma aracı — farklı geometriler ve yük dağılımları ile hızlı deneme yapılmasını sağlar — hem de TUSAŞ LiftUp programı sunumları için bir gösteri platformu olarak hizmet eder. Kullanıcılar grid çözünürlüğünü (128'den 2048'e) interaktif olarak ayarlayabilir, multi-grid hızlandırmayı açıp kapatabilir ve çözümün gömülü bir log konsolu aracılığıyla canlı olarak yakınsama sürecini izleyebilir.

Benchmark ResultsBenchmark Sonuçları

The solver was tested on two NVIDIA RTX A4000 GPUs connected via Gigabit Ethernet (1 Gbps) — deliberately a bandwidth-limited setup to stress-test the async optimizations. Grid sizes ranged from $128^2$ to $2048^2$, with the TUSAŞ KAAN aircraft cross-section as the test geometry under Dirichlet conditions. Çözücü, Gigabit Ethernet (1 Gbps) üzerinden bağlanan iki NVIDIA RTX A4000 GPU'da test edildi — asenkron optimizasyonları stres-test etmek için kasıtlı olarak bant genişliği sınırlı bir kurulum. Grid boyutları $128^2$'den $2048^2$'ye değişti; test geometrisi olarak TUSAŞ KAAN uçağı kesiti kullanıldı.

2048×2048 Configuration Comparison2048×2048 Konfigürasyon Karşılaştırması

Speedup vs Grid SizeGrid Boyutuna Göre Hızlanma

At small grids ($128^2$, $256^2$), communication overhead dominates and multi-GPU is slower. From $512^2$ onward, the compute-to-communication ratio favors distribution: 2.17× for point charge and 2.47× for Gaussian at $2048^2$. With multi-grid warmup, the total speedup reaches 3.2× (solution time dropping from 290.8s to 88.7s for Gaussian charge). Küçük gridlerde ($128^2$, $256^2$) iletişim yükü baskındır. $512^2$'den itibaren hesaplama/iletişim oranı dağıtımı destekler: $2048^2$'de nokta yük için 2.17×, Gaussian için 2.47×. Multi-grid ısınmasıyla toplam hızlanma 3.2×'e ulaşır (Gaussian yük için çözüm süresi 290.8s'den 88.7s'ye düşer).

Solution VisualizationÇözüm Görselleştirmesi

The potential field $\phi(x,y)$ computed on the heart domain with $\rho(x,y) = \sin(x)\cos(y)$ — Dirichlet boundary ($\phi=0$ on the contour): Kalp domaininde $\rho(x,y) = \sin(x)\cos(y)$ ile hesaplanan potansiyel alanı $\phi(x,y)$ — Dirichlet sınır koşulu ($\phi=0$ kontür üzerinde):