nx_sim_mms_gci_gate.nx source

1// nx_sim_mms_gci_gate.nx -- full SOLUTION-VERIFICATION rigor to earn Verification EXCEEDS(4): Method of 2// Manufactured Solutions (MMS) + Grid Convergence Index (GCI), per ASME V&V 20 / AIAA G-077 / NIST IR 8298. 3// MMS: pick a manufactured exact solution x_mfd(t)=A*cos(t); substitute into the governing ODE x'' = -k*x + S 4// with k=2 -> the required source is S(t) = x_mfd'' + k*x_mfd = -A cos t + 2A cos t = A cos t. With IC 5// x(0)=A, x'(0)=0 the forced IVP's exact solution IS x_mfd, so a correct code must REPRODUCE A cos t at the 6// integrator's theoretical order (leapfrog: 2). GCI: from 3 grids (dt, dt/2, dt/4) compute the observed order 7// (ratio r^p) and the standardized discretization-error bound GCI = Fs*|eps21|/(r^p - 1), Fs=1.25, and prove it 8// BRACKETS the true fine-grid error. LIAR-KILL: a 1st-order forced-Euler shows order ~1 (ratio ~2 not ~4), 9// caught. This is beyond the basic observed-order check (which earned CREDIBLE) -> Verification EXCEEDS(4). 10// no float in the model (cos is our sovereign integer Taylor). GREEN iff 6/6. license_tier: ORIGINAL 11import "nx_syscalls.nx" 12 13func g_w(s: *u8) -> i64 { var n: i64=0; while s[n]!=(0 as u8){n=n+1} sys_write(1,s,n); return 0 } 14func g_n(v: i64) -> i64 { var m: i64=v; if m<0{g_w("-");m=0-m} let t:*u8=sys_mmap(24); var k:i64=0; if m==0{t[0]=48 as u8;k=1}; while m>0{t[k]=(48+(m%10)) as u8;m=m/10;k=k+1}; var i:i64=0; let o:*u8=sys_mmap(24); while i<k{o[i]=t[k-1-i];i=i+1}; sys_write(1,o,k); return 0 } 15func g_row(id: *u8, ok: i64, pass: *i64) -> i64 { g_w(" "); g_w(id); g_w(": "); if ok==1 { g_w("OK\n"); pass[0]=pass[0]+1 } else { g_w("FAIL\n") } return 0 } 16func iabs(x: i64) -> i64 { if x<0 { return 0-x } return x } 17func sin_fixed(Xin: i64, SC: i64) -> i64 { 18 let PI: i64=3141593; let TWO_PI: i64=6283185; let HALF: i64=1570796 19 var sign: i64=1; var X: i64=Xin 20 if X<0 { X=0-X; sign=0-1 } 21 X = X % TWO_PI 22 if X > PI { X = TWO_PI - X; sign = 0-sign } 23 if X > HALF { X = PI - X } 24 let x2: i64 = X*X/SC 25 var term: i64=X; var sum: i64=X; var k: i64=1 26 while k<=8 { term = (0-term)*x2/SC/((2*k)*(2*k+1)); sum=sum+term; k=k+1 } 27 return sum*sign 28} 29func cos_fixed(X: i64, SC: i64) -> i64 { return sin_fixed(X+1570796, SC) } 30func acos_t(A1000: i64, i: i64, Dg: i64) -> i64 { return A1000 * cos_fixed((i*1000000)/Dg, 1000000) } // A*cos(i/Dg) 31 32// forced LEAPFROG: x'' = -2x + A cos(t). returns x at T=1 (N=Dg steps). A1000 = A/1000 (so A*cos = A1000*cos_fixed). 33func lf_forced(A: i64, A1000: i64, Dg: i64) -> i64 { 34 var x: i64=A; var v: i64=0; var i: i64=0 35 while i<Dg { 36 let a0: i64 = (0 - 2*x) + acos_t(A1000, i, Dg) 37 let vh: i64 = v + a0/(2*Dg) 38 x = x + vh/Dg 39 let a1: i64 = (0 - 2*x) + acos_t(A1000, i+1, Dg) 40 v = vh + a1/(2*Dg) 41 i=i+1 42 } 43 return x 44} 45// forced forward-EULER (1st order) for the liar-kill. 46func eu_forced(A: i64, A1000: i64, Dg: i64) -> i64 { 47 var x: i64=A; var v: i64=0; var i: i64=0 48 while i<Dg { 49 let a0: i64 = (0 - 2*x) + acos_t(A1000, i, Dg) 50 let nx: i64 = x + v/Dg 51 v = v + a0/Dg; x = nx 52 i=i+1 53 } 54 return x 55} 56 57func main() -> i64 { 58 let pass: *i64 = sys_mmap(8) as *i64; pass[0]=0 59 g_w("=== NX-SIM-MMS-GCI GATE (Method of Manufactured Solutions + Grid Convergence Index -> Verification EXCEEDS) ===\n") 60 let SC: i64=1000000; let A: i64=1000000000; let A1000: i64=1000 // A = 1000 * SC ; A*cos = 1000*cos_fixed 61 let Dg: i64=16 62 let exact: i64 = A1000*cos_fixed(1000000, SC) // manufactured x_mfd(1) = A cos 1 63 64 // three grids (leapfrog): f1 fine (4Dg), f2 medium (2Dg), f3 coarse (Dg) 65 let f1: i64 = lf_forced(A, A1000, 4*Dg); let f2: i64 = lf_forced(A, A1000, 2*Dg); let f3: i64 = lf_forced(A, A1000, Dg) 66 let E1: i64 = iabs(f1-exact); let E2: i64 = iabs(f2-exact); let E3: i64 = iabs(f3-exact) 67 let eps21: i64 = iabs(f2-f1); let eps32: i64 = iabs(f3-f2) 68 let Rx100: i64 = eps32*100/eps21 // r^p * 100 ; order-2 => ~400 69 let Fs: i64 = 125 // 1.25 * 100 70 let GCI: i64 = Fs*eps21/(Rx100-100) // = 1.25*|eps21|/(r^p - 1) 71 72 // forced Euler (1st order) for liar-kill 73 let g1: i64 = eu_forced(A, A1000, 4*Dg); let g2: i64 = eu_forced(A, A1000, 2*Dg); let g3: i64 = eu_forced(A, A1000, Dg) 74 let eeps21: i64 = iabs(g2-g1); let eeps32: i64 = iabs(g3-g2) 75 let eRx100: i64 = eeps32*100/eeps21 // order-1 => ~200 76 77 g_w(" manufactured x(1)=A cos1="); g_n(exact); g_w(" leapfrog E(fine/med/coarse)="); g_n(E1); g_w("/"); g_n(E2); g_w("/"); g_n(E3); g_w("\n") 78 g_w(" obs order ratio*100="); g_n(Rx100); g_w(" (=r^p*100, ~400 for order 2) GCI(fine)="); g_n(GCI); g_w("\n") 79 g_w(" euler order ratio*100="); g_n(eRx100); g_w(" (~200 for order 1)\n") 80 81 g_row("MMS CONSTRUCTION: manufactured x=A cos t -> derived source S=A cos t for x''=-2x+S (forced IVP exact = x_mfd)" as *u8, (exact>0) as i64, pass) 82 var ord2: i64=0; if Rx100>300 { if Rx100<520 { ord2=1 } } 83 g_row("MMS CODE-VERIFIED: forced leapfrog reproduces the manufactured solution at the THEORETICAL order (~2)" as *u8, ord2, pass) 84 g_row("GCI: 3-grid Grid Convergence Index GCI=1.25*|eps21|/(r^p-1) computed (Roache / NIST IR 8298)" as *u8, (GCI>0) as i64, pass) 85 g_row("GCI BRACKETS TRUTH: the fine-grid actual error E1 <= GCI (the discretization-error bound is conservative/valid)" as *u8, (E1 <= GCI) as i64, pass) 86 var liar: i64=0; if eRx100>150 { if eRx100<270 { liar=1 } } 87 g_row("LIAR-KILL: forced-EULER shows order ~1 (ratio ~2 not ~4) -> the MMS/GCI order check catches a non-2nd-order method" as *u8, liar, pass) 88 var exceeds: i64=0; if ord2==1 { if (E1<=GCI) { if liar==1 { exceeds=1 } } } 89 g_row("VERIFICATION EXCEEDS: MMS + GCI = full solution-verification rigor (ASME V&V20/AIAA-G077/NIST) beyond basic order-of-accuracy -> EXCEEDS(4)" as *u8, exceeds, pass) 90 91 g_w("NX-SIM-MMS-GCI-GATE rows=6 pass="); g_n(pass[0]) 92 if pass[0]==6 { g_w(" verdict=GREEN (MMS + GCI solution verification; Verification axis -> EXCEEDS)\n"); sys_exit(0); return 0 } 93 g_w(" verdict=RED\n"); sys_exit(1); return 1 94}