#include "fiko/dokum/sph_dolum.hpp"

#include <algorithm>
#include <chrono>
#include <cmath>
#include <cstdio>
#include <mutex>
#ifdef _OPENMP
#include <omp.h>
#endif

#include "SPlisHSPlasH/BoundaryModel_Akinci2012.h"
#include "SPlisHSPlasH/DFSPH/TimeStepDFSPH.h"
#include "SPlisHSPlasH/Emitter.h"
#include "SPlisHSPlasH/EmitterSystem.h"
#include "SPlisHSPlasH/FluidModel.h"
#include "SPlisHSPlasH/Simulation.h"
#include "SPlisHSPlasH/StaticRigidBody.h"
#include "SPlisHSPlasH/SurfaceTension/SurfaceTension_Becker2007.h"
#include "SPlisHSPlasH/Vorticity/VorticityConfinement.h"
#include "SPlisHSPlasH/TimeManager.h"
#include "SPlisHSPlasH/TimeStep.h"
#include "SPlisHSPlasH/Viscosity/Viscosity_Standard.h"
#include "SPlisHSPlasH/Viscosity/Viscosity_Weiler2018.h"
#include "SPlisHSPlasH/XSPH.h"
#include "Utilities/Counting.h"
#include "Utilities/Logger.h"
#include "Utilities/Timing.h"

// SPlisHSPlasH tekil durum tanımları (tüm süreçte bir kez).
INIT_LOGGING
INIT_TIMING
INIT_COUNTING

namespace fiko::dokum {

namespace {

using Real = ::Real;  // SPlisHSPlasH Common.h (float)

// Zincir sözleşmesi (dolum.cpp): döküm bölgesi = Ic + Sinir (temas kabuğu).
inline bool metalMi(const VoxelKafes& kafes, u32 idx) {
  return kafes.tip[idx] == HucreTip::Ic || kafes.tip[idx] == HucreTip::Sinir;
}
// SPH akışkan alanı = yalnız gerçek metal içi (Sinir = kalıp kabuğu, katı duvar).
inline bool icMi(const VoxelKafes& kafes, u32 idx) {
  return kafes.tip[idx] == HucreTip::Ic;
}

// yon → ortonormal taban (col0 = yon) — emitter yönelimi.
Matrix3r yonTaban(Vec3f yon) {
  Vector3r x(yon.x, yon.y, yon.z);
  if (x.norm() < 1e-8f) x = Vector3r(0, 0, -1);
  x.normalize();
  Vector3r yardimci = (std::abs(x.z()) < 0.9f) ? Vector3r(0, 0, 1) : Vector3r(1, 0, 0);
  Vector3r y = yardimci.cross(x).normalized();
  Vector3r z = x.cross(y);
  Matrix3r m;
  m.col(0) = x; m.col(1) = y; m.col(2) = z;
  return m;
}

// Metal↔kalıp yüzlerinden Akinci statik sınır partikülleri (dünya=metre).
std::vector<Vector3r> sinirPartikulOrnekle(const VoxelKafes& kafes, f64 birimOlcek,
                                           Vec3d kok, f64 hMesh, f64 aralikM) {
  std::vector<Vector3r> nokta;
  const i32 nx = static_cast<i32>(kafes.nx());
  const i32 ny = static_cast<i32>(kafes.ny());
  const i32 nz = static_cast<i32>(kafes.nz());
  const f64 hFiz = hMesh * birimOlcek;
  const i32 nf = std::max(1, static_cast<i32>(std::ceil(hFiz / aralikM)));
  const f64 adim = hFiz / nf;
  const i32 yonlar[6][3] = {{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,0,-1}};
  for (i32 k = 0; k < nz; ++k)
    for (i32 j = 0; j < ny; ++j)
      for (i32 i = 0; i < nx; ++i) {
        const u32 idx = kafes.indeks(i, j, k);
        if (!icMi(kafes, idx)) continue;
        const Vec3d cm = kafes.merkez(i, j, k);
        const Vector3r merkezM(static_cast<Real>(cm.x * birimOlcek),
                               static_cast<Real>(cm.y * birimOlcek),
                               static_cast<Real>(cm.z * birimOlcek));
        for (const auto& d : yonlar) {
          const i32 ni = i + d[0], nj = j + d[1], nk = k + d[2];
          bool kalip = true;
          if (ni >= 0 && nj >= 0 && nk >= 0 && ni < nx && nj < ny && nk < nz)
            kalip = !icMi(kafes, kafes.indeks(static_cast<u32>(ni),
                                              static_cast<u32>(nj),
                                              static_cast<u32>(nk)));
          if (!kalip) continue;
          // Yüz düzlemi: merkez + 0.5h·d; düzlem içi nf×nf ızgara.
          // İKİ katman (ikincisi dışa ofset + yarım-adım şaşırtmalı): hızlı
          // partikülün büyük CFL adımında tek katmandan tünellemesini keser.
          Vector3r n(static_cast<Real>(d[0]), static_cast<Real>(d[1]), static_cast<Real>(d[2]));
          Vector3r u = (std::abs(n.z()) < 0.5f) ? Vector3r(0, 0, 1) : Vector3r(1, 0, 0);
          Vector3r v = n.cross(u);
          const Vector3r yuzMerkez = merkezM + n * static_cast<Real>(0.5 * hFiz);
          for (i32 katman = 0; katman < 2; ++katman) {
            const Vector3r taban = yuzMerkez + n * static_cast<Real>(katman * aralikM);
            const f64 kaydir = katman ? 0.5 * adim : 0.0;
            for (i32 a = 0; a < nf; ++a)
              for (i32 b = 0; b < nf; ++b) {
                const f64 oa = (a + 0.5) * adim - 0.5 * hFiz + kaydir;
                const f64 ob = (b + 0.5) * adim - 0.5 * hFiz + kaydir;
                nokta.push_back(taban + u * static_cast<Real>(oa) + v * static_cast<Real>(ob));
              }
          }
        }
      }
  (void)kok;
  return nokta;
}

// Tekil SPlisHSPlasH durumunu temizle (yeniden çağrılabilirlik).
void simTemizle() {
  if (SPH::Simulation::hasCurrent()) {
    delete SPH::Simulation::getCurrent();  // dtor: TimeManager dahil siler
    SPH::Simulation::setCurrent(nullptr);
  }
  SPH::TimeManager::setCurrent(nullptr);
}

}  // namespace

// ---- Efektif partikül yarıçapı: CPU + GPU yollarının TEK KAYNAĞI -----------
// İstenen (partikulYaricapMm) ≠ alınan: h/2 tavanı binning sözleşmesinden gelir
// (partikül hücreden büyük olamaz — nokta-binning onu göremez). Ayrıntı ve
// ölçülen maliyet sonuçları için sph_dolum.hpp bildirimi.
f64 sphEtkinYaricapM(const VoxelKafes& kafes, const SphDolumAyar& ayar) {
  if (!kafes.gecerliMi()) return 0.0;
  const f64 hFiz = kafes.h();  // metre
  if (!(hFiz > 0) || !std::isfinite(hFiz)) return 0.0;
  const f64 istenen =
      (ayar.partikulYaricapMm > 0) ? ayar.partikulYaricapMm * 1e-3 : 0.25 * hFiz;
  return std::clamp(istenen, 1e-5, 0.5 * hFiz);
}

bool sphYaricapKirpildiMi(const VoxelKafes& kafes, const SphDolumAyar& ayar) {
  if (!(ayar.partikulYaricapMm > 0) || !kafes.gecerliMi()) return false;
  const f64 istenen = ayar.partikulYaricapMm * 1e-3;
  const f64 etkin = sphEtkinYaricapM(kafes, ayar);
  return etkin > 0 && std::abs(istenen - etkin) > 1e-9 * std::max(1.0, istenen);
}

SphEmitPlan sphEmitPlani(const VoxelKafes& kafes, const DokumKurulum& kurulum,
                         const Alasim& alasim, const SphDolumAyar& ayar) {
  SphEmitPlan p;
  const f64 r = sphEtkinYaricapM(kafes, ayar);
  if (!(r > 0) || !(alasim.yogunluk > 0)) return p;
  for (const auto& g : kurulum.girisler)
    if (g.aktif && g.kutleKg > 0 && g.sureS > 0) {
      p.toplamHacimM3 += g.kutleKg / alasim.yogunluk;
      p.sonPencereS = std::max(p.sonPencereS, g.gecikmeS + g.sureS);
    }
  if (p.toplamHacimM3 <= 0) return p;
  const f64 cap = 2.0 * r;
  const f64 partikulHacim = 0.8 * cap * cap * cap;  // FluidModel m_V kalıbı
  p.gerekliPartikul = 1.3 * p.toplamHacimM3 / partikulHacim;
  p.havuz = static_cast<u32>(
      std::clamp(p.gerekliPartikul, 64.0, static_cast<f64>(ayar.maksPartikul)));
  // Kırpılma = TAVAN kısıyor (u32 kesmesinin <1 partikül farkı kırpılma değil).
  p.kirpildi = p.gerekliPartikul > static_cast<f64>(ayar.maksPartikul) + 0.5;
  return p;
}

SphDolumSonuc sphDolumCoz(const VoxelKafes& kafes, const DokumKurulum& kurulum,
                          const Alasim& alasim, const SphDolumAyar& ayar,
                          IlerlemeCB ilerleme) {
  SphDolumSonuc s;
  const u32 N = kafes.hucreSayi();
  auto& d = s.dolum;
  d.varisZamani.assign(N, -1.0);
  d.varisT.assign(N, 0.0);
  d.coldShut.assign(N, 0);
  d.havaKapani.assign(N, 0);
  if (!kafes.gecerliMi() || !kurulum.aktifGirisVar()) return s;

  u32 toplamMetal = 0, icSayi = 0;
  for (u32 i = 0; i < N; ++i) {
    if (metalMi(kafes, i)) ++toplamMetal;
    if (icMi(kafes, i)) ++icSayi;
  }
  d.toplamMetal = toplamMetal;
  if (icSayi == 0) return s;

  // Geometri ölçekleri: merkez() mesh birimi, h() metre.
  const f64 hFiz = kafes.h();
  const Vec3d m0 = kafes.merkez(0, 0, 0);
  f64 hMesh = hFiz;
  if (kafes.nx() > 1) hMesh = kafes.merkez(1, 0, 0).x - m0.x;
  else if (kafes.ny() > 1) hMesh = kafes.merkez(0, 1, 0).y - m0.y;
  else if (kafes.nz() > 1) hMesh = kafes.merkez(0, 0, 1).z - m0.z;
  const f64 birimOlcek = hFiz / hMesh;
  const Vec3d kok{m0.x - 0.5 * hMesh, m0.y - 0.5 * hMesh, m0.z - 0.5 * hMesh};

  // Partikül yarıçapı [m] — TEK KAYNAK (GPU yolu da aynı fonksiyonu çağırır).
  const f64 r = sphEtkinYaricapM(kafes, ayar);
  if (!(r > 0)) return s;
  const f64 cap = 2.0 * r;

  // Emitter planı + havuz boyutu (TEK KAYNAK — GPU yolu da aynı planı kullanır).
  const SphEmitPlan plan = sphEmitPlani(kafes, kurulum, alasim, ayar);
  if (plan.toplamHacimM3 <= 0) return s;
  const f64 partikulHacim = 0.8 * cap * cap * cap;  // FluidModel m_V kalıbı
  const u32 maksEmit = plan.havuz;
  const f64 maksSure =
      (ayar.maksSure > 0) ? ayar.maksSure : 1.5 * plan.sonPencereS + 0.5;
  s.maksSureS = maksSure;

  // SPlisHSPlasH tekil → serileştir.
  static std::mutex kilit;
  std::scoped_lock kilitle(kilit);
  simTemizle();

#ifdef _OPENMP
  const int ompOnceki = omp_get_max_threads();
  if (ayar.deterministik) omp_set_num_threads(1);
#endif

  const auto baslangic = std::chrono::steady_clock::now();
  bool iptal = false;
  {
    SPH::Simulation* sim = SPH::Simulation::getCurrent();
    sim->init(static_cast<Real>(r), false);
    // Adaptif CFL: kararlı akışta dt tavana (maksAdimDt) çıkar → az adım/hızlı.
    // KRİTİK: kararlılık implicit viskoziteden gelir (bkz SphDolumAyar default'ları);
    // kararsız çözümde v_max fırlar → CFL dt'yi minik tutar → adım patlar. Bu yüzden
    // default'lar implicit + sıkı çözücü (272 adım) — gevşek/explicit = 2765 adım.
    sim->setValue<int>(SPH::Simulation::CFL_METHOD, SPH::Simulation::ENUM_CFL_STANDARD);
    sim->setValue<Real>(SPH::Simulation::CFL_MAX_TIMESTEPSIZE, static_cast<Real>(ayar.maksAdimDt));
    sim->setValue<Real>(SPH::Simulation::CFL_FACTOR,
                        static_cast<Real>(std::clamp(ayar.cflFaktor, 0.05, 1.0)));
    Vec3f gy = kurulum.graviteYon.birim();
    if (gy.uzunlukKare() < 1e-12f) gy = {0, 0, -1};
    Real gvec[3] = {static_cast<Real>(gy.x * ayar.gravite),
                    static_cast<Real>(gy.y * ayar.gravite),
                    static_cast<Real>(gy.z * ayar.gravite)};
    sim->setVecValue(SPH::Simulation::GRAVITATION, gvec);
    sim->setValue<int>(SPH::Simulation::BOUNDARY_HANDLING_METHOD,
                       SPH::Simulation::ENUM_AKINCI2012);

    sim->addFluidModel("dolum", 0u, nullptr, nullptr, nullptr, maksEmit);
    SPH::FluidModel* fm = sim->getFluidModel(0);
    fm->setDensity0(static_cast<Real>(alasim.yogunluk));
    if (ayar.viskoziteImplicit) {
      // Weiler2018 implicit Laplace viskozite [WKBB18, CGF 2018]: matris-serbest
      // CG çözümü → yüksek katsayıda koşulsuz kararlı (Standard >0.02 patlıyordu).
      // Duvar katsayısı = akışkan katsayısı (no-slip sürtünme, SPlisHSPlasH kalıbı).
      fm->setViscosityMethod(SPH::Viscosity_Weiler2018::METHOD_NAME);
      if (auto* visk = fm->getViscosityBase()) {
        visk->setValue<Real>(SPH::Viscosity_Weiler2018::VISCOSITY_COEFFICIENT,
                             static_cast<Real>(ayar.viskozite));
        visk->setValue<Real>(SPH::Viscosity_Weiler2018::VISCOSITY_COEFFICIENT_BOUNDARY,
                             static_cast<Real>(ayar.viskozite));
      }
    } else {
      fm->setValue<int>(SPH::FluidModel::VISCOSITY_METHOD, 1);  // Standard (explicit)
      if (auto* visk = fm->getViscosityBase())
        visk->setValue<Real>(SPH::Viscosity_Standard::VISCOSITY_COEFFICIENT,
                             static_cast<Real>(ayar.viskozite));
    }
    if (ayar.yuzeyGerilim > 0) {
      // Opsiyonel yüzey gerilimi (Becker & Teschner 2007, kohezyon kuvveti).
      fm->setSurfaceTensionMethod(SPH::SurfaceTension_Becker2007::METHOD_NAME);
      if (auto* yg = fm->getSurfaceTensionBase())
        yg->setValue<Real>(SPH::SurfaceTension_Becker2007::SURFACE_TENSION,
                           static_cast<Real>(ayar.yuzeyGerilim));
    }
    // xSPH (Monaghan 1992): hız alanını komşu ortalamasıyla yumuşatır —
    // basınç tabanlı SPH'in sayısal girdap SÜPÜRMESİNİ engeller, akışı
    // "canlı" tutar (sıvı yapışkan/donuk görünmez). computeXSPH() her zaman
    // çağrılır; katsayı 0 iken etkisiz (kütüphanenin kendi varsayılanı).
    if (ayar.xsphKatsayi > 0) {
      if (auto* xs = fm->getXSPH()) {
        xs->setValue<Real>(SPH::XSPH::FLUID_COEFFICIENT,
                           static_cast<Real>(ayar.xsphKatsayi));
        xs->setValue<Real>(SPH::XSPH::BOUNDARY_COEFFICIENT,
                           static_cast<Real>(ayar.xsphKatsayi * 0.5));
      }
    }
    // Girdap sıkıştırma (Spalding 1995): SPH sayısal girdabı kaybettiği için
    // geri beslenir → dar boğaz çıkışında gerçek türbülanslı karışım.
    if (ayar.girdapKatsayi > 0) {
      fm->setVorticityMethod(SPH::VorticityConfinement::METHOD_NAME);
      if (auto* gd = fm->getVorticityBase())
        gd->setValue<Real>(SPH::VorticityConfinement::VORTICITY_COEFFICIENT,
                           static_cast<Real>(ayar.girdapKatsayi));
    }

    // Statik sınır (kalıp duvarları).
    std::vector<Vector3r> sinir =
        sinirPartikulOrnekle(kafes, birimOlcek, kok, hMesh, 1.5 * r);
    s.sinirPartikulSayi = static_cast<u32>(sinir.size());
    auto* rb = new SPH::StaticRigidBody();
    rb->setWorldSpacePosition(Vector3r::Zero());
    rb->setWorldSpaceRotation(Matrix3r::Identity());
    rb->setPosition0(Vector3r::Zero());
    rb->setRotation0(Quaternionr::Identity());
    auto* bm = new SPH::BoundaryModel_Akinci2012();
    bm->initModel(rb, static_cast<unsigned int>(sinir.size()), sinir.data());
    sim->addBoundaryModel(bm);

    // Girişler → dairesel emitter (kendi zaman penceresi + debiden hız).
    // Hız kalibrasyonu: emitter d/v aralıkla nIn partikül dilimi yayar →
    // hacimsel debi = nIn·Vp·v/d. Hedef debiyi TAM tutturmak için v buradan.
    std::vector<Vector3r> emitPoz;  // kaynak kimliği ataması (Material Trace)
    for (const auto& g : kurulum.girisler) {
      if (!g.aktif || g.kutleKg <= 0 || g.sureS <= 0) continue;
      const f64 debiM3s = (g.kutleKg / alasim.yogunluk) / g.sureS;
      const u32 genislik =
          std::max(1u, static_cast<u32>(std::round(g.capMm * 1e-3 / cap)));
      u32 nIn = 0;  // Emitter::emitParticlesCircle dilim örüntüsünün aynısı
      {
        const f64 yariCap2 = 0.25 * genislik * genislik * cap * cap;
        const f64 bas = -0.5 * (genislik - 1) * cap;
        for (u32 a = 0; a < genislik; ++a)
          for (u32 b = 0; b < genislik; ++b) {
            const f64 x = bas + a * cap, y = bas + b * cap;
            if (x * x + y * y <= yariCap2) ++nIn;
          }
        nIn = std::max(nIn, 1u);
      }
      const f64 hiz = debiM3s * cap / (nIn * partikulHacim);
      const Matrix3r rot = yonTaban(g.yon);
      Vector3r poz(static_cast<Real>(g.konum.x * birimOlcek),
                   static_cast<Real>(g.konum.y * birimOlcek),
                   static_cast<Real>(g.konum.z * birimOlcek));
      poz += rot.col(0) * static_cast<Real>(sim->getSupportRadius());  // kaviteye it
      emitPoz.push_back(poz);
      fm->getEmitterSystem()->addEmitter(genislik, 1u, poz, rot,
                                         static_cast<Real>(std::max(hiz, 1e-6)), 1u);
      SPH::Emitter* em = fm->getEmitterSystem()->getEmitters().back();
      em->setObjectId(0u);
      em->setEmitStartTime(static_cast<Real>(g.gecikmeS));
      em->setEmitEndTime(static_cast<Real>(g.gecikmeS + g.sureS));
    }

    sim->setSimulationMethod(static_cast<int>(SPH::SimulationMethods::DFSPH));
    // DFSPH basınç çözücü ayarı: görsel dolum için tolerans gevşet + iterasyon
    // tavanı düşür + divergence solver'ı kapat (adım başına ikinci pahalı blok).
    // Not: setValue MAX_ERROR yüzdedir (SPlisHSPlasH iç ölçek).
    if (auto* ts = dynamic_cast<SPH::TimeStepDFSPH*>(sim->getTimeStep())) {
      ts->setValue<Real>(SPH::TimeStepDFSPH::MAX_ERROR,
                         static_cast<Real>(std::clamp(ayar.basincHata, 1e-3, 1.0)));
      ts->setValue<unsigned int>(SPH::TimeStepDFSPH::MAX_ITERATIONS,
                                 std::max(1u, ayar.basincMaksIter));
      ts->setValue<bool>(SPH::TimeStepDFSPH::USE_DIVERGENCE_SOLVER,
                         ayar.divergenceCozucu);
      if (ayar.divergenceCozucu) {
        ts->setValue<Real>(SPH::TimeStepDFSPH::MAX_ERROR_V,
                           static_cast<Real>(std::clamp(ayar.divergenceHata, 1e-3, 2.0)));
        ts->setValue<unsigned int>(SPH::TimeStepDFSPH::MAX_ITERATIONS_V,
                                   std::max(1u, ayar.divergenceMaksIter));
      }
    }
    sim->performNeighborhoodSearchSort();
    sim->updateBoundaryVolume();
    SPH::TimeManager::getCurrent()->setTimeStepSize(
        static_cast<Real>(std::min(0.001, ayar.maksAdimDt)));

    // ---- Zaman döngüsü ------------------------------------------------------
    const u32 kareSayi = std::max(1u, ayar.kareSayi);
    const f64 kareAralik = maksSure / kareSayi;
    f64 sonrakiKare = kareAralik;
    const f64 hucreHacim = kafes.hucreHacim();
    std::vector<u16> sayac(N, 0);
    std::vector<u8> dolu(N, 0);  // anlık kalıcı-dolum durumu (histerezis)
    std::vector<f64> hizToplam(N, 0.0);
    std::vector<u32> hizAdet(N, 0);
    s.hizMaks.assign(N, 0.0f);
    u32 dolan = 0;
    f64 t = 0, sonDegisimT = 0;

    // ---- Isı + oksit partikül durumu (indeks = FluidModel havuzu; döngü içi
    //      sort YOK ve emitter reuse kapalı → indeksler kararlı) --------------
    const bool isi = ayar.isiKuplaj;
    const f64 tLo = std::min(ayar.kalipSicaklik, ayar.dokumSicaklik);
    const f64 tHi = std::max(ayar.kalipSicaklik, ayar.dokumSicaklik);
    std::vector<f32> pT(maksEmit, static_cast<f32>(ayar.dokumSicaklik));
    std::vector<f32> pTyeni(isi ? maksEmit : 0, 0.f);
    std::vector<u8> pKaynak(maksEmit, 0);  // doğduğu emitter (Material Trace)
    std::vector<u8> pOksit(maksEmit, 0);   // 1 = oksitlenmiş (kalıcı, Campbell bifilm)
    const f64 kritikHiz =
        (ayar.kritikHiz > 0) ? ayar.kritikHiz : kritikGirisHiz(alasim);
    // Adım-içi hücre istatistiği: cephe sıcaklığı + iki-kaynak birleşme takibi.
    std::vector<f32> tAdim(N, 0.f);           // hücredeki partikül T toplamı
    std::vector<f32> cepheT(N, -1.f);         // varış anı ort. T [°C] (-1 = yok)
    const bool birlesmeIzle = isi && emitPoz.size() >= 2;
    std::vector<u8> kynA(N, 255), kynB(N, 255);   // hücredeki ilk iki farklı kaynak
    std::vector<f32> tKynA(N, 0.f), tKynB(N, 0.f);
    std::vector<u16> sKynA(N, 0), sKynB(N, 0);
    std::vector<u8> ciftOnce(N, 0);  // önceki adımda hücre iki-kaynaklı mıydı

#ifdef FIKO_SPH_PERF
    auto* tsPerf = dynamic_cast<SPH::TimeStepDFSPH*>(sim->getTimeStep());
    f64 perfTopMs = 0;
#endif
    while (t < maksSure) {
      // Komşu listeleri adım BAŞINDA aranır → adım sonrası sorgular yalnız
      // adım öncesi aktif partiküller (nOnce) için geçerli (emisyon adım
      // SONUNDA olur, yeni doğanların komşu listesi bir sonraki adımda).
      const unsigned int nOnce = fm->numActiveParticles();
      const f64 tOnce = t;
#ifdef FIKO_SPH_PERF
      const auto adimBas = std::chrono::steady_clock::now();
#endif
      sim->getTimeStep()->step();
      ++s.adimSayi;
      t = SPH::TimeManager::getCurrent()->getTime();
      const f64 dtAdim = t - tOnce;
      const unsigned int nAktif = fm->numActiveParticles();
#ifdef FIKO_SPH_PERF
      {
        const f64 adimMs =
            std::chrono::duration<f64, std::milli>(
                std::chrono::steady_clock::now() - adimBas).count();
        perfTopMs += adimMs;
        if (s.adimSayi % 50 == 0 || s.adimSayi == 1) {
          const int it = tsPerf ? tsPerf->getNumIterations() : -1;
          std::fprintf(stderr,
              "[perf] adim=%llu t=%.4f dt=%.5f n=%u iter=%d adimMs=%.2f ort=%.2f\n",
              (unsigned long long)s.adimSayi, t, dtAdim, nAktif, it, adimMs,
              perfTopMs / s.adimSayi);
        }
      }
#endif

      // Yeni doğan partiküller: T = dokumSicaklik + en yakın emitter kaynağı.
      for (unsigned int p = nOnce; p < nAktif; ++p) {
        pT[p] = static_cast<f32>(ayar.dokumSicaklik);
        if (emitPoz.size() > 1) {
          const Vector3r& x = fm->getPosition(p);
          u8 enYakin = 0;
          f64 enK = 1e30;
          for (size_t e = 0; e < emitPoz.size(); ++e) {
            const f64 dd = (x - emitPoz[e]).squaredNorm();
            if (dd < enK) { enK = dd; enYakin = static_cast<u8>(e); }
          }
          pKaynak[p] = enYakin;
        }
      }

      // ---- Isı adımı (explicit Euler; dt << h²/α → kararlı) -----------------
      if (isi && nOnce > 0 && dtAdim > 0) {
        const Real hDestek = sim->getSupportRadius();
        const f64 eta2 = 0.01 * static_cast<f64>(hDestek) * hDestek;
        const f64 Lc = std::cbrt(static_cast<f64>(fm->getVolume(0)));
        // (a) Duvar Robin: q = U·(T−T_kalıp); partikül temas alanı A≈V^(2/3)
        //     → dT/dt = −U·(T−T_kalıp)/(ρ·c·V^(1/3)).
        const f64 duvarKat = ayar.duvarU / (alasim.yogunluk * alasim.c * Lc);
        // (b) Partiküller arası iletim — Cleary & Monaghan 1999 (J. Comput.
        //     Phys. 148, "Conduction modelling using smoothed particle
        //     hydrodynamics"): c·dT_i/dt = Σ_j m_j/(ρ_iρ_j)·(4k_ik_j/(k_i+k_j))
        //     ·(T_i−T_j)·(x_ij·∇W_ij)/(|x_ij|²+0.01h²). Tek malzeme → 2k.
        const f64 ikiK = 2.0 * alasim.k;
        const unsigned int nPS = sim->numberOfPointSets();
#ifdef _OPENMP
#pragma omp parallel for schedule(static)
#endif
        for (i64 pi = 0; pi < static_cast<i64>(nOnce); ++pi) {
          const auto p = static_cast<unsigned int>(pi);
          const Vector3r& xi = fm->getPosition(p);
          const f64 rhoI = fm->getDensity(p);
          const f64 Ti = pT[p];
          f64 dT = 0;
          const unsigned int nf = sim->numberOfNeighbors(0, 0, p);
          for (unsigned int jj = 0; jj < nf; ++jj) {
            const unsigned int q = sim->getNeighbor(0, 0, p, jj);
            const Vector3r xij = xi - fm->getPosition(q);
            const f64 r2 = xij.squaredNorm();
            dT += static_cast<f64>(fm->getMass(q)) * ikiK * (Ti - pT[q]) *
                  static_cast<f64>(xij.dot(sim->gradW(xij))) /
                  (rhoI * fm->getDensity(q) * (r2 + eta2) * alasim.c);
          }
          bool temas = false;  // temas = herhangi bir sınır partikülü komşuluğu
          for (unsigned int pid = 1; pid < nPS && !temas; ++pid)
            temas = sim->numberOfNeighbors(0, pid, p) > 0;
          if (temas) dT -= duvarKat * (Ti - ayar.kalipSicaklik);
          // Enerji sınırı: kalıp sıcaklığı altına inemez, döküm T üstüne çıkamaz.
          pTyeni[p] = static_cast<f32>(std::clamp(Ti + dtAdim * dT, tLo, tHi));
        }
        std::copy(pTyeni.begin(), pTyeni.begin() + nOnce, pT.begin());
      }

      // Partikülleri hücrelere bağla (doluluk + hız + T + kaynak istatistiği).
      std::fill(sayac.begin(), sayac.end(), u16{0});
      // Hız toplamı/sayacı da ADIMLIK olmalı — aksi halde tum sim boyunca
      // birikir (s.hizOrt ve kare hiz alani yalnizca SON ADIMDA dolmus hucrelerin
      // ort. hizini verir → akis alani zaman icinde yanlis/donuk kalir).
      std::fill(hizToplam.begin(), hizToplam.end(), 0.0);
      std::fill(hizAdet.begin(), hizAdet.end(), 0u);
      if (isi) std::fill(tAdim.begin(), tAdim.end(), 0.f);
      if (birlesmeIzle) {
        std::fill(kynA.begin(), kynA.end(), u8{255});
        std::fill(kynB.begin(), kynB.end(), u8{255});
        std::fill(tKynA.begin(), tKynA.end(), 0.f);
        std::fill(tKynB.begin(), tKynB.end(), 0.f);
        std::fill(sKynA.begin(), sKynA.end(), u16{0});
        std::fill(sKynB.begin(), sKynB.end(), u16{0});
      }
      for (unsigned int p = 0; p < nAktif; ++p) {
        const Vector3r& x = fm->getPosition(p);
        const f64 v = fm->getVelocity(p).norm();
        // Oksit/bifilm (Campbell "Castings": serbest yüzey + hız > kritik →
        // yüzey oksit filmi katlanıp sıvıya karışır; işaret KALICI). Serbest
        // yüzey vekili: akışkan komşu sayısı eşik altı (FLOW-3D defect
        // tracking'in partikül-tabanlı karşılığı).
        if (!pOksit[p] && p < nOnce && v > kritikHiz &&
            sim->numberOfNeighbors(0, 0, p) < ayar.yuzeyKomsuEsik)
          pOksit[p] = 1;
        const f64 mi = (x.x() / birimOlcek - kok.x) / hMesh;
        const f64 mj = (x.y() / birimOlcek - kok.y) / hMesh;
        const f64 mk = (x.z() / birimOlcek - kok.z) / hMesh;
        if (mi < 0 || mj < 0 || mk < 0) continue;
        const u32 i = static_cast<u32>(mi), j = static_cast<u32>(mj), k = static_cast<u32>(mk);
        if (i >= kafes.nx() || j >= kafes.ny() || k >= kafes.nz()) continue;
        const u32 idx = kafes.indeks(i, j, k);
        if (sayac[idx] < 65535) ++sayac[idx];
        hizToplam[idx] += v;
        ++hizAdet[idx];
        if (v > s.hizMaks[idx]) s.hizMaks[idx] = static_cast<f32>(v);
        if (v > s.maksHiz) s.maksHiz = v;
        if (isi) tAdim[idx] += pT[p];
        if (birlesmeIzle) {
          const u8 ky = pKaynak[p];
          if (kynA[idx] == 255 || kynA[idx] == ky) {
            kynA[idx] = ky; tKynA[idx] += pT[p]; ++sKynA[idx];
          } else if (kynB[idx] == 255 || kynB[idx] == ky) {
            kynB[idx] = ky; tKynB[idx] += pT[p]; ++sKynB[idx];
          }
        }
      }

      // Varış zamanı = İLK eşik aşımı (seviye gelişi). Histerezisli dolu[]
      // durumu yalnız NİHAİ geçerlilik için: sonda boş kalan (sıçrayıp
      // boşalmış) hücrenin varışı çözüm sonunda geçersiz kılınır.
      const f64 doluPay = static_cast<f64>(fm->getVolume(0)) / hucreHacim;
      const f64 coldShutEsikT = alasim.tLikidus - ayar.coldShutDT;
      for (u32 c = 0; c < N; ++c) {
        if (!icMi(kafes, c)) continue;
        // GERÇEK cold-shut: hücre İLK KEZ iki FARKLI kaynağın malzemesini
        // birlikte görüyor (cephe birleşme anı) ve İKİ tarafın da T'si <
        // likidus−ΔT → oksitli cepheler kaynaşamaz (Campbell cold shut /
        // bifilm dikişi). Tek-kaynak→çift-kaynak GEÇİŞİNDE denetlenir; dolu
        // eşiğine bağlı değil (sıçramayla erken dolan hücrede birleşme
        // dolu-geçişiyle çakışmaz).
        if (birlesmeIzle) {
          const bool cift = kynB[c] != 255 && sKynA[c] > 0 && sKynB[c] > 0;
          if (cift && !ciftOnce[c] &&
              tKynA[c] / sKynA[c] < coldShutEsikT &&
              tKynB[c] / sKynB[c] < coldShutEsikT)
            d.coldShut[c] = 1;
          ciftOnce[c] = cift ? u8{1} : u8{0};
        }
        const f64 oran = sayac[c] * doluPay;
        if (!dolu[c] && oran >= ayar.doluEsik) {
          dolu[c] = 1;
          if (d.varisZamani[c] < 0.0) {
            d.varisZamani[c] = t;  // ilk aşım kalır
            // Cephe sıcaklığı: varış anındaki partikül T ortalaması →
            // katılaşma baslangicT köprüsü (dolum.varisT'e sonda yazılır).
            if (isi && sayac[c] > 0)
              cepheT[c] = tAdim[c] / static_cast<f32>(sayac[c]);
          }
          ++dolan;
          sonDegisimT = t;
        } else if (dolu[c] && oran < 0.6 * ayar.doluEsik) {
          dolu[c] = 0;
          --dolan;
          sonDegisimT = t;
        }
      }

      // Animasyon karesi.
      // BUGDUUZTMESI: `while (t >= sonrakiKare)` kareyi ANLIK GÖNDERIYORDU — aradaki
      // zaman farkinda cozucu IKI ADIM daha ilerlemis oluyordu. dt buyuk (sakin
      // faz, dt≈maksAdimDt=4 ms) iken kare araligi (maksSure/kareSayi) cok kucuk
      // kalir → bir ciftte BIRDEN FAZLA kare ayni (guncel) dolulukla yazilir,
      // sonraki kareler atlanti. Sonuc: animasyon "donuk/zıplayan" — akis kare
      // basina buyuk adim atar, aradaki akis gorunmez. Ayrica sayac[] bu adimda
      // hesaplandigi icin kare t anini DEGIL, t+dt anini temsil ediyordu.
      // DOGRUSU: kare anini t'ye YAKLASTIR (kare araliginin geriden gelmesine izin
      // ver) ve kacirilan kareleri ATLAMA — hepsini ayni (dogru) dolulukla doldur.
      // Bu, zaman karelerini ADIM sayisina degil SIM'a yaslar (deterministik).
      while (sonrakiKare <= t && s.kareler.size() < kareSayi) {
        SphDolumKare kare;
        kare.sure = sonrakiKare;
        kare.doluluk.resize(N);
        for (u32 c = 0; c < N; ++c)
          kare.doluluk[c] = static_cast<u8>(
              std::min(255.0, sayac[c] * doluPay * 255.0));
        // Hucre-bası ORT. hız alanı (GÖRSEL — akis "yasiyor" izlenimi icin).
        // hizAdet bu adimda SIFIRLANDI ve partikul dongusunde dolduruldu → su an
        // ki deger t+dt anini temsil eder (dolulukla ayni kayan pencere).
        // Bos hucreler 0 (yuzey=0 hiz, mandala hizli → kontrast kirilir).
        kare.hizOrt.resize(N);
        for (u32 c = 0; c < N; ++c)
          kare.hizOrt[c] = hizAdet[c] ? static_cast<f32>(hizToplam[c] / hizAdet[c]) : 0.0f;
        if (nAktif > 0) {
          f64 toplamT = 0;
          for (unsigned int p = 0; p < nAktif; ++p) toplamT += pT[p];
          kare.ortT = static_cast<f32>(toplamT / nAktif);
        }
        // GÖRSEL partikül kaydı: hücre binning'i (yukarıda) fizik için aynen
        // kaldı; burada konumun KENDİSİ saklanır (jet/sıçrama/damla hücreye
        // ezilmeden görünsün). Kapalıysa vektör boş → eski davranış.
        if (ayar.partikulKaydet && nAktif > 0) {
          const unsigned int adim = std::max(1u, ayar.partikulAltOrnek);
          kare.partikul.reserve((nAktif + adim - 1) / adim);
          for (unsigned int p = 0; p < nAktif; p += adim) {
            const Vector3r& x = fm->getPosition(p);
            SphPartikul sp;
            sp.konum = {static_cast<f32>(x.x() / birimOlcek),
                        static_cast<f32>(x.y() / birimOlcek),
                        static_cast<f32>(x.z() / birimOlcek)};
            sp.hiz = static_cast<f32>(fm->getVelocity(p).norm());
            sp.sicaklik = pT[p];
            kare.partikul.push_back(sp);
          }
        }
        s.kareler.push_back(std::move(kare));
        sonrakiKare += kareAralik;
      }

      if (ilerleme && (s.adimSayi % 16 == 0)) {
        const f64 oran = std::min(1.0, std::max(t / maksSure,
            static_cast<f64>(dolan) / icSayi));
        if (!ilerleme(oran, "SPH dolum")) { iptal = true; break; }
      }
      // Erken çıkış: akış pencereleri bitti + dolum durumu duruldu.
      if (t > plan.sonPencereS && (dolan == icSayi || t - sonDegisimT > 0.2)) break;
    }

    // Sonda boş kalan (sıçrayıp boşalmış) hücrelerin varışı geçersiz.
    for (u32 c = 0; c < N; ++c)
      if (icMi(kafes, c) && !dolu[c] && d.varisZamani[c] >= 0.0)
        d.varisZamani[c] = -1.0;
    dolan = 0;
    for (u32 c = 0; c < N; ++c)
      if (icMi(kafes, c) && d.varisZamani[c] >= 0.0) ++dolan;

    // Sinir (temas kabuğu) hücrelere varışı komşu iç hücreden yay — zincir
    // sözleşmesi (Ic+Sinir = döküm bölgesi) korunur, misrun doğru kalır.
    u32 dolanKabuk = 0;
    for (u32 k = 0; k < kafes.nz(); ++k)
      for (u32 j = 0; j < kafes.ny(); ++j)
        for (u32 i = 0; i < kafes.nx(); ++i) {
          const u32 idx = kafes.indeks(i, j, k);
          if (kafes.tip[idx] != HucreTip::Sinir || d.varisZamani[idx] >= 0.0)
            continue;
          f64 enErken = -1.0;
          u32 enErkenIdx = idx;
          const i64 nx = kafes.nx(), ny = kafes.ny(), nz = kafes.nz();
          const i64 komsu[6][3] = {{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,0,-1}};
          for (const auto& dk : komsu) {
            const i64 ni = i + dk[0], nj = j + dk[1], nk = k + dk[2];
            if (ni < 0 || nj < 0 || nk < 0 || ni >= nx || nj >= ny || nk >= nz)
              continue;
            const u32 nidx = kafes.indeks(static_cast<u32>(ni), static_cast<u32>(nj),
                                          static_cast<u32>(nk));
            if (icMi(kafes, nidx) && d.varisZamani[nidx] >= 0.0 &&
                (enErken < 0 || d.varisZamani[nidx] < enErken)) {
              enErken = d.varisZamani[nidx];
              enErkenIdx = nidx;
            }
          }
          if (enErken >= 0) {
            d.varisZamani[idx] = enErken;
            cepheT[idx] = cepheT[enErkenIdx];  // kabuk cephe T = kaynak iç hücre
            ++dolanKabuk;
          }
        }

    // ---- Sonuçları haritala -------------------------------------------------
    d.dolanHucre = dolan + dolanKabuk;
    d.dolanHacim = d.dolanHucre * hucreHacim;
    f64 sonVaris = 0;
    f64 tMin = ayar.dokumSicaklik;
    d.sira.clear();
    d.sira.reserve(dolan);
    for (u32 c = 0; c < N; ++c) {
      if (d.varisZamani[c] >= 0.0) {
        d.sira.push_back(c);
        sonVaris = std::max(sonVaris, d.varisZamani[c]);
        // Cephe sıcaklığı → varisT (katılaşma baslangicT köprüsü; ısı kuplajı
        // kapalıysa/kaydedilmediyse döküm sıcaklığı sabiti).
        d.varisT[c] = (isi && cepheT[c] >= 0.f) ? static_cast<f64>(cepheT[c])
                                                : ayar.dokumSicaklik;
        tMin = std::min(tMin, d.varisT[c]);
      } else if (d.coldShut[c]) {
        d.coldShut[c] = 0;  // sonda boşalan hücre: cold-shut işareti düşer
      }
    }
    s.cepheTMin = tMin;
    d.coldShutSayi = 0;
    for (u32 c = 0; c < N; ++c)
      if (d.coldShut[c]) ++d.coldShutSayi;
    std::sort(d.sira.begin(), d.sira.end(), [&](u32 a, u32 b) {
      return d.varisZamani[a] < d.varisZamani[b];
    });
    d.toplamSure = sonVaris;
    d.misrun = !iptal && (d.dolanHucre < toplamMetal);
    d.hazir = !iptal;

    s.hizOrt.assign(N, 0.0f);
    for (u32 c = 0; c < N; ++c)
      if (hizAdet[c] > 0)
        s.hizOrt[c] = static_cast<f32>(hizToplam[c] / hizAdet[c]);
    s.partikulSayi = fm->numActiveParticles();
    s.partikulKutleKg = static_cast<f64>(s.partikulSayi) *
                        static_cast<f64>(fm->getVolume(0)) * alasim.yogunluk;

    // ---- Oksit yoğunluğu alanı: oksitli partikülün SON konumu → hücre ------
    // (oksitli hacim / hücre hacmi; bifilmler nihai yerlerinde kusur adayı).
    s.oksit.assign(N, 0.f);
    const f64 oksitPay = static_cast<f64>(fm->getVolume(0)) / hucreHacim;
    for (unsigned int p = 0; p < fm->numActiveParticles(); ++p) {
      if (!pOksit[p]) continue;
      ++s.oksitliPartikul;
      const Vector3r& x = fm->getPosition(p);
      const f64 mi = (x.x() / birimOlcek - kok.x) / hMesh;
      const f64 mj = (x.y() / birimOlcek - kok.y) / hMesh;
      const f64 mk = (x.z() / birimOlcek - kok.z) / hMesh;
      if (mi < 0 || mj < 0 || mk < 0) continue;
      const u32 i = static_cast<u32>(mi), j = static_cast<u32>(mj), k = static_cast<u32>(mk);
      if (i >= kafes.nx() || j >= kafes.ny() || k >= kafes.nz()) continue;
      s.oksit[kafes.indeks(i, j, k)] += static_cast<f32>(oksitPay);
    }
  }
  simTemizle();
#ifdef _OPENMP
  if (ayar.deterministik) omp_set_num_threads(ompOnceki);
#endif

  s.cozumSuresiSn =
      std::chrono::duration<f64>(std::chrono::steady_clock::now() - baslangic).count();
  s.iptal = iptal;
  s.hazir = !iptal;
  return s;
}

}  // namespace fiko::dokum