#include "fiko/dokum/katilasma_fdm.hpp"

#include <algorithm>
#include <cmath>

namespace fiko::dokum {

// ---- Entalpi ↔ sıcaklık (H = ρ[cT + L(1−fs(T))]) ---------------------------
f64 entalpi(const Alasim& a, f64 T) {
  return a.yogunluk * (a.c * T + a.gizliIsi * (1.0 - a.katiOran(T)));
}

f64 sicaklik(const Alasim& a, f64 H) {
  // H(T) parçalı: katı (fs=1) H=ρcT ve sıvı (fs=0) H=ρ(cT+L) LİNEER → analitik
  // ters (çözücünün en sıcak yolu; eski geniş-aralık 60-iter bisection darboğazdı).
  const f64 rc = a.yogunluk * a.c;
  const f64 hSol = rc * a.tSolidus;                       // fs=1 → H=ρcT
  if (H <= hSol) return H / rc;                           // tam katı
  const f64 hLik = a.yogunluk * (a.c * a.tLikidus + a.gizliIsi);  // fs=0
  if (H >= hLik) return (H / a.yogunluk - a.gizliIsi) / a.c;      // tam sıvı
  const f64 dT = a.tLikidus - a.tSolidus;
  if (dT <= 1e-9) return a.tSolidus;                      // izotermal mushy
  if (a.fsModel == FsModel::Lineer || a.fsModel == FsModel::Lever) {
    // fs=(Tl−T)/dT lineer → H=ρ[cT+L(T−Ts)/dT] → analitik lineer ters.
    return (H / a.yogunluk + a.gizliIsi * a.tSolidus / dT) / (a.c + a.gizliIsi / dT);
  }
  // Scheil (nonlineer): DAR aralık [Ts,Tl] bisection (32 iter → <1e-8 °C).
  f64 lo = a.tSolidus, hi = a.tLikidus;
  for (int it = 0; it < 32; ++it) {
    const f64 mid = 0.5 * (lo + hi);
    if (entalpi(a, mid) < H) lo = mid; else hi = mid;
  }
  return 0.5 * (lo + hi);
}

KatilasmaSonuc katilasmaCoz(VoxelKafes& kafes, const Alasim& alasim,
                            const Kalip& kalip, const FdmAyar& ayar,
                            f64 niyamaEsik, IlerlemeCB ilerleme,
                            const std::vector<f64>* baslangicT,
                            const std::vector<f64>* hucreH) {
  KatilasmaSonuc r;
  if (!kafes.gecerliMi()) return r;

  const u32 nx = kafes.nx(), ny = kafes.ny(), nz = kafes.nz();
  const u32 N = kafes.hucreSayi();
  const f64 h = kafes.h();          // hücre boyutu [m]
  const f64 h2 = h * h;
  const f64 rho = alasim.yogunluk, c = alasim.c, k = alasim.k;
  const f64 alpha = alasim.yayinim();

  // Metal hücre kümesi (Ic + Sinir = döküm bölgesi). Dış = kalıp emici sınır.
  auto metalMi = [&](u32 idx) {
    return kafes.tip[idx] == HucreTip::Ic || kafes.tip[idx] == HucreTip::Sinir;
  };

  // ---- Başlangıç koşulu: metal = dökümSıcaklık, kalıp = tBaslangic ----
  // Döküm köprüsü: baslangicT verilirse per-hücre (dolum sonrası soğumuş cephe)
  // başlangıç sıcaklığı kullanılır; en az likidus (metal döküm anında hep sıvı).
  const bool kopru = baslangicT && baslangicT->size() == N;
  u32 metalSayi = 0;
  for (u32 i = 0; i < N; ++i) {
    if (metalMi(i)) {
      const f64 T0 = kopru ? std::max((*baslangicT)[i], alasim.tLikidus) : ayar.dokumSicaklik;
      kafes.T[i] = T0;
      kafes.H[i] = entalpi(alasim, T0);
      kafes.fs[i] = alasim.katiOran(T0);
      ++metalSayi;
    } else {
      kafes.T[i] = kalip.tBaslangic;
      kafes.H[i] = 0.0;
      kafes.fs[i] = 0.0;
    }
  }
  if (metalSayi == 0) return r;

  // ---- Kararlı zaman adımı (3B açık difüzyon: dt ≤ cfl·h²/(6α)) ----
  // 6 komşu → kararlılık sınırı h²/(6α). Robin sınır kaybı ek kararsızlık:
  //   T güncelleme katsayısı 1 − dt·(6α/h² + 6·U/(ρc·h)) ≥ 0 olmalı. En kötü
  //   U (t→0 arayüz) ile birleşik sınır → güvenli dt (cfl<1 pay).
  const f64 rc = rho * c;
  // Chill maskesi (per-hücre arayüz HTC; 0 = varsayılan). Kararlılık için en
  // kötü (en büyük) HTC ile dt sınırlanır.
  const f64* hH = (hucreH && hucreH->size() == N) ? hucreH->data() : nullptr;
  f64 Umaks = kalip.arayuzH;  // t→0 en büyük arayüz (en kötü kararlılık)
  if (hH)
    for (u32 i = 0; i < N; ++i) Umaks = std::max(Umaks, hH[i]);
  const f64 difOran = 6.0 * alpha / h2;                 // iç difüzyon oranı [1/s]
  const f64 sinirOran = 6.0 * Umaks / (rc * h);         // en kötü Robin oranı [1/s]
  f64 dt = ayar.cflGuvenlik / std::max(difOran + sinirOran, 1e-12);
  dt = std::min(dt, ayar.maksZaman);
  if (!(dt > 0) || !std::isfinite(dt)) return r;

  // Post-process için gerekli izleme alanları.
  r.katilasmaSuresi.assign(N, -1.0);   // -1 = henüz katılaşmadı
  r.sonSicaklik.assign(N, 0.0);
  std::vector<f64> oncekiT(N, 0.0);
  // Niyama örnekleme: fs solidus civarına (fs≈0.7) ilk ulaşıldığı anda G,Ṙ yakala.
  std::vector<f64> niyamaG(N, 0.0), niyamaR(N, 0.0);
  std::vector<u8> niyamaAlindi(N, 0);
  const f64 niyamaFs = 0.7;

  const f64 hMm = h * 1000.0;  // hücre boyutu mm (Niyama birimi mm)

  std::vector<f64> yeniH(N);
  f64 t = 0.0;
  u32 adim = 0;
  u32 katilasan = 0;  // fs≥hedefFs olan metal hücre sayısı (hepsi → dur)
  const f64 hedefFs = std::min(ayar.durmaKatiOran, 0.99);

  auto komsuT = [&](u32 i, u32 j, u32 kk, f64 kendiT) -> f64 {
    // Sınır dışı → kendi (Neumann sıfır-akı, simetri).
    if (i >= nx || j >= ny || kk >= nz) return kendiT;
    return kafes.T[kafes.indeks(i, j, kk)];
  };

  // OpenMP: küçük gridlerde (testler) thread açma maliyeti işten büyük → if-eşiği.
  const bool paralel = N > 20000;
  u32 niyamaKalan = metalSayi;  // tüm hücreler örneklenince Niyama taraması atlanır

  while (adim < ayar.maksAdim && t < ayar.maksZaman) {
    // Efektif arayüz U(t) (yarı-sonsuz kalıp emici, zamanla azalır).
    const f64 U = ayar.yariSonsuzKalip ? kalip.etkinH(t) : kalip.arayuzH;

#pragma omp parallel for schedule(static) if (paralel)
    for (int kk = 0; kk < (int)nz; ++kk)
      for (u32 j = 0; j < ny; ++j)
        for (u32 i = 0; i < nx; ++i) {
          const u32 idx = kafes.indeks(i, j, (u32)kk);
          if (!metalMi(idx)) continue;
          const f64 Tc = kafes.T[idx];
          // Chill etki bölgesi: bu hücrenin kalıp arayüz HTC'si maskeden gelir
          // (sabit — metal chill zamanla doymaz varsayımı, D1).
          const f64 Uc = (hH && hH[idx] > 0.0) ? hH[idx] : U;

          // Laplacian: 6 komşu. Metal komşu → iletim; kalıp/dış komşu → Robin kaybı.
          f64 laplasyonAkis = 0.0;  // ∇·(k∇T) toplamı (birim hacim başına W/m³·h²/k biçimi)
          f64 sinirKayip = 0.0;     // Robin arayüzden ısı kaybı [W/m³ · ... ]
          const int di[6] = {-1, 1, 0, 0, 0, 0};
          const int dj[6] = {0, 0, -1, 1, 0, 0};
          const int dk[6] = {0, 0, 0, 0, -1, 1};
          for (int f = 0; f < 6; ++f) {
            const int ni = (int)i + di[f], nj = (int)j + dj[f], nk = kk + dk[f];
            const bool disarda =
                ni < 0 || nj < 0 || nk < 0 || ni >= (int)nx || nj >= (int)ny || nk >= (int)nz;
            const u32 nidx = disarda ? idx : kafes.indeks(ni, nj, nk);
            if (!disarda && metalMi(nidx)) {
              // İletim: k·(Tn−Tc)/h²  (açık difüzyon).
              laplasyonAkis += (kafes.T[nidx] - Tc);
            } else {
              // Kalıp/dış yüz: Robin q = U·(Tc − T_kalip) [W/m²]. Hacimsel: q·A/V = q/h.
              sinirKayip += Uc * (Tc - kalip.tBaslangic);
            }
          }
          // İç difüzyon ısı akısı [W/m³] = k·Σ(Tn−Tc)/h².
          const f64 difuzyon = k * laplasyonAkis / h2;
          const f64 sinir = sinirKayip / h;  // A/V = 6·(h²)/(h³)= 6/h; her yüz q/h
          const f64 dH = (difuzyon - sinir) * dt;  // ΔH = (∇·k∇T − sınır) dt [J/m³]
          yeniH[idx] = kafes.H[idx] + dH;
        }

    // H → T, fs güncelle + katılaşma anı (per-idx bağımsız → paralel + reduction).
    int yeniKatilasan = 0;
#pragma omp parallel for schedule(static) reduction(+ : yeniKatilasan) if (paralel)
    for (int sidx = 0; sidx < (int)N; ++sidx) {
      const u32 idx = (u32)sidx;
      if (!metalMi(idx)) continue;
      oncekiT[idx] = kafes.T[idx];
      kafes.H[idx] = yeniH[idx];
      const f64 Tn = sicaklik(alasim, kafes.H[idx]);
      kafes.T[idx] = Tn;
      const f64 fsn = alasim.katiOran(Tn);
      kafes.fs[idx] = fsn;
      if (fsn >= hedefFs && r.katilasmaSuresi[idx] < 0) {
        r.katilasmaSuresi[idx] = t + dt;
        ++yeniKatilasan;
      }
    }
    katilasan += (u32)yeniKatilasan;

    // Niyama: fs solidus civarı (~0.7) geçilirken G=|∇T|, Ṙ=|dT/dt| örnekle.
    // Tüm metal hücreler örneklendiyse tarama tümüyle atlanır (geç adımlar bedava).
    if (niyamaKalan > 0) {
      int alinan = 0;
#pragma omp parallel for schedule(static) reduction(+ : alinan) if (paralel)
      for (int kk = 0; kk < (int)nz; ++kk)
        for (u32 j = 0; j < ny; ++j)
          for (u32 i = 0; i < nx; ++i) {
            const u32 idx = kafes.indeks(i, j, (u32)kk);
            if (!metalMi(idx) || niyamaAlindi[idx]) continue;
            if (kafes.fs[idx] < niyamaFs) continue;
            // Sıcaklık gradyanı büyüklüğü (merkezi fark, mm cinsinden).
            const f64 tx = komsuT(i + 1, j, kk, kafes.T[idx]) - komsuT(i - 1, j, kk, kafes.T[idx]);
            const f64 ty = komsuT(i, j + 1, kk, kafes.T[idx]) - komsuT(i, j - 1, kk, kafes.T[idx]);
            const f64 tz = komsuT(i, j, kk + 1, kafes.T[idx]) - komsuT(i, j, kk - 1, kafes.T[idx]);
            const f64 G = std::sqrt(tx * tx + ty * ty + tz * tz) / (2.0 * hMm);  // °C/mm
            const f64 Rdot = std::abs(kafes.T[idx] - oncekiT[idx]) / dt;         // °C/s
            niyamaG[idx] = G;
            niyamaR[idx] = std::max(Rdot, 1e-9);
            niyamaAlindi[idx] = 1;
            ++alinan;
          }
      niyamaKalan -= std::min((u32)alinan, niyamaKalan);
    }

    t += dt;
    ++adim;
    // Chvorinov: TÜM metal hücreler katılaşınca dur (son katılaşan = hot-spot).
    // Ortalama-fs değil (merkez hücre erken durmayı bozar → V/A² ilişkisi kaybolur).
    if (ilerleme && (adim % 32 == 0)) {
      if (!ilerleme(std::clamp((f64)katilasan / (f64)metalSayi, 0.0, 1.0),
                    "Katılaşma çözülüyor"))
        break;  // iptal
    }
    if (katilasan >= metalSayi) break;
  }

  // ---- Sonuç doldur: Niyama, SDAS, hot-spot ----
  r.niyama.assign(N, 0.0);
  r.sdas.assign(N, 0.0);
  r.niyamaMin = 1e30; r.niyamaMax = 0.0;
  r.toplamSure = t;
  r.hotSpotSure = -1.0;
  bool niyamaVar = false;

  for (u32 idx = 0; idx < N; ++idx) {
    r.sonSicaklik[idx] = kafes.T[idx];
    if (!metalMi(idx)) continue;
    // Hot-spot = en geç katılaşan hücre (besleyici ipucu).
    const f64 ts = r.katilasmaSuresi[idx];
    if (ts > r.hotSpotSure) { r.hotSpotSure = ts; r.hotSpotIndeks = idx; }
    if (ts > r.toplamSure) r.toplamSure = ts;
    // Niyama Ny = G/√Ṙ  [(°C·s)^½/mm].
    if (niyamaAlindi[idx]) {
      const f64 Ny = niyamaG[idx] / std::sqrt(niyamaR[idx]);
      r.niyama[idx] = Ny;
      r.niyamaMin = std::min(r.niyamaMin, Ny);
      r.niyamaMax = std::max(r.niyamaMax, Ny);
      if (Ny < niyamaEsik) ++r.niyamaKritik;
      niyamaVar = true;
      // SDAS: λ₂ = K·Ṙ^(−1/3) [m] → µm. (soğuma hızı katılaşma anındaki Ṙ).
      const f64 lambda = alasim.sdasK * std::pow(niyamaR[idx], -1.0 / 3.0);
      r.sdas[idx] = lambda * 1e6;
    }
  }
  if (!niyamaVar) { r.niyamaMin = 0; r.niyamaMax = 0; }
  r.hazir = true;
  return r;
}

}  // namespace fiko::dokum