diff --git a/docs/shell-transient-heat-transfer/numerical-oracle-scratch.md b/docs/shell-transient-heat-transfer/numerical-oracle-scratch.md new file mode 100644 index 0000000..e13b4a9 --- /dev/null +++ b/docs/shell-transient-heat-transfer/numerical-oracle-scratch.md @@ -0,0 +1,425 @@ +# 독립 numerical oracle scratch + +이 파일은 numerical-review.md의 사전 계산을 재현하는 문서 부록이다. Solver 또는 +production test가 아니며 외부 package·reference·출력 파일을 사용하지 않는다. +Python 표준 라이브러리만 사용한다. Fraction은 oracle에서 입력 double의 exact 산술을 +확인하는 용도이며 production arbitrary-precision dependency를 요구하지 않는다. + +T1 root에서 다음 명령은 각 Python block을 순서대로 메모리에서 실행한다. + +```powershell +@' +from pathlib import Path +import re +p = Path('docs/shell-transient-heat-transfer/numerical-oracle-scratch.md') +blocks = re.findall(r'```python\n(.*?)```', p.read_text(encoding='utf-8-sig'), re.S) +assert len(blocks) == 4 +scope = {} +for block in blocks: + exec(compile(block, str(p), 'exec'), scope) +print('ALL SCRATCH BLOCKS COMPLETED') +'@ | python - +``` + +실행 exit는 0이어야 한다. Assertion 범위는 exact algebra와 독립 oracle 자체에 한정된다. +수렴/이산 오차의 출력 숫자 자체를 production tolerance로 사용하지 않고 본 보고서의 +roundoff·conditioning·오차 전파 기준을 적용한다. NumPy 부재 probe exit1은 실행 실패 +이력으로 보고서에 남겼으며 아래 코드에는 NumPy를 쓰지 않는다. + +## A. Slab spectral·quadrature·공간/시간 oracle + +```python +import math,json +def zeros(n,m=None): return [[0.0]*(n if m is None else m) for _ in range(n)] +def tr(A):return list(map(list,zip(*A))) +def mm(A,B):return [[sum(a*b for a,b in zip(row,col)) for col in zip(*B)]for row in A] +def mv(A,x):return [sum(a*b for a,b in zip(row,x)) for row in A] +def solve(A,b): + a=[list(r)+[v]for r,v in zip(A,b)];n=len(b) + for k in range(n): + p=max(range(k,n),key=lambda i:abs(a[i][k]));a[k],a[p]=a[p],a[k] + d=a[k][k] + for j in range(k,n+1):a[k][j]/=d + for i in range(n): + if i!=k: + d=a[i][k] + for j in range(k,n+1):a[i][j]-=d*a[k][j] + return [r[-1]for r in a] +def chol(A): + n=len(A);L=zeros(n) + for i in range(n): + for j in range(i+1): + v=A[i][j]-sum(L[i][k]*L[j][k]for k in range(j)) + L[i][j]=math.sqrt(v)if i==j else v/L[j][j] + return L +def eig(A): + a=[r[:]for r in A];n=len(a);v=[[float(i==j)for j in range(n)]for i in range(n)] + for z in range(100*n*n): + p,q=max(((i,j)for i in range(n)for j in range(i+1,n)),key=lambda ij:abs(a[ij[0]][ij[1]])) + if abs(a[p][q])<1e-15*sum(abs(a[i][i])for i in range(n)):break + phi=.5*math.atan2(2*a[p][q],a[q][q]-a[p][p]);c=math.cos(phi);s=math.sin(phi) + J=[[float(i==j)for j in range(n)]for i in range(n)] + J[p][p]=J[q][q]=c;J[p][q]=s;J[q][p]=-s + a=mm(mm(tr(J),a),J);v=mm(v,J) + ix=sorted(range(n),key=lambda i:a[i][i]) + return [a[i][i]for i in ix],[[v[j][i]for i in ix]for j in range(n)] +def slabmat(ne): + n=2*ne+1;C=zeros(n);K=zeros(n);l=1/ne + for e in range(ne): + for i in range(3): + for j in range(3): + C[2*e+i][2*e+j]+=l/30*[[4,2,-1],[2,16,2],[-1,2,4]][i][j] + K[2*e+i][2*e+j]+=[[7,-8,1],[-8,16,-8],[1,-8,7]][i][j]/(3*l) + return C,K +C,K=slabmat(2);Cf=[r[1:]for r in C[1:]];Kf=[r[1:]for r in K[1:]] +L=chol(Cf);Li=tr([solve(L,[float(i==j)for i in range(4)])for j in range(4)]) +lam,V=eig(mm(mm(Li,Kf),tr(Li)));W=mm(tr(Li),V) +yp=[-1-x for x in solve(Cf,[r[0]for r in C[1:]])] +coef=mv(tr(W),mv(Cf,yp)) +def spectral(t,n=None): + amp=[a*(math.exp(-l*t)if n is None else math.exp(-n*math.log1p(l*t/n)))for a,l in zip(coef,lam)] + T=[1]+[1+x for x in mv(W,amp)] + if n is None:dT=[0]+mv(W,[-l*x for l,x in zip(lam,amp)]) + else: + oldamp=[x*(1+l*t/n)for l,x in zip(lam,amp)] + old=[1]+[1+x for x in mv(W,oldamp)] + dT=[(a-b)/(t/n)for a,b in zip(T,old)] + E=sum(mv(C,T));R=mv(C,dT)[0]+mv(K,T)[0] + q=[] + for e,ss in [(0,-1),(0,0),(0,1),(1,-1),(1,0),(1,1)]: + q.append(-4*sum(a*b for a,b in zip([ss-.5,-2*ss,ss+.5],T[2*e:2*e+3]))) + return T,q,E,R +def exact(t): + M=80 + T=[1-sum(4/((2*m+1)*math.pi)*math.sin((2*m+1)*math.pi*x/2)*math.exp(-((2*m+1)*math.pi/2)**2*t)for m in range(M))for x in [0,.25,.5,.75,1]] + q=[sum(2*math.cos((2*m+1)*math.pi*x/2)*math.exp(-((2*m+1)*math.pi/2)**2*t)for m in range(M))for x in [0,.25,.5,.5,.75,1]] + E=1-sum(8/((2*m+1)*math.pi)**2*math.exp(-((2*m+1)*math.pi/2)**2*t)for m in range(M)) + return T,q,E,q[0] +def errors(a,b):return [max(abs(x-y)for x,y in zip(a[0],b[0])),max(abs(x-y)for x,y in zip(a[1],b[1])),abs(a[2]-b[2]),abs(a[3]-b[3])] +out={"lambda":lam,"projected_initial_free":[1+x for x in yp],"slab":[]} +for t in [.001,.01,.05,.1,.5,1]: + s=spectral(t);x=exact(t) + out["slab"].append({"t":t,"penetration":math.sqrt(t),"spatial_error":errors(s,x),"BE4096_temporal":errors(spectral(t,4096),s),"BE8192_temporal":errors(spectral(t,8192),s),"spectral":s}) +print(json.dumps(out,indent=2)) + +def gauss(n): + x=[];w=[] + for i in range(1,n+1): + z=math.cos(math.pi*(i-.25)/(n+.5)) + for _ in range(30): + p0,p1=1.,z + for j in range(2,n+1):p0,p1=p1,((2*j-1)*z*p1-(j-1)*p0)/j + dp=n*(z*p1-p0)/(z*z-1);nz=z-p1/dp + if abs(nz-z)<2e-16: z=nz;break + z=nz + x.append(z);w.append(2/((1-z*z)*dp*dp)) + return list(zip(x,w)) +def dot(x,y):return sum(a*b for a,b in zip(x,y)) +def cross(a,b):return [a[1]*b[2]-a[2]*b[1],a[2]*b[0]-a[0]*b[2],a[0]*b[1]-a[1]*b[0]] +def geometry(X,xi,eta): + signs=[(-1,-1),(1,-1),(1,1),(-1,1)] + N=[(1+s*xi)*(1+t*eta)/4 for s,t in signs] + D=[[s*(1+t*eta)/4,t*(1+s*xi)/4]for s,t in signs] + a=[sum(D[i][0]*X[i][j]for i in range(4))for j in range(3)] + b=[sum(D[i][1]*X[i][j]for i in range(4))for j in range(3)] + W=cross(a,b);J=math.sqrt(dot(W,W));norm=[z/J for z in W] + aa,ab,bb=dot(a,a),dot(a,b),dot(b,b);det=J*J + ra=[(bb*x-ab*y)/det for x,y in zip(a,b)];rb=[(aa*y-ab*x)/det for x,y in zip(a,b)] + G=[[D[i][0]*ra[j]+D[i][1]*rb[j]for j in range(3)]for i in range(4)] + return N,G,J,norm +def surface(X,n): + M=zeros(4);S=zeros(4) + for xi,wx in gauss(n): + for eta,wy in gauss(n): + N,G,J,no=geometry(X,xi,eta) + for i in range(4): + for j in range(4): + M[i][j]+=N[i]*N[j]*J*wx*wy;S[i][j]+=dot(G[i],G[j])*J*wx*wy + return M,S +Z0,Z1=slabmat(2) +# h=1: Z0 physical mass already; Z1 physical normal stiffness already. +def mats(X,n): + M,S=surface(X,n);C=zeros(20);K=zeros(20) + for i in range(4): + for j in range(4): + for a in range(5): + for b in range(5): + C[5*i+a][5*j+b]=M[i][j]*Z0[a][b] + K[5*i+a][5*j+b]=S[i][j]*Z0[a][b]+M[i][j]*Z1[a][b] + return C,K +def fn(A):return math.sqrt(sum(x*x for r in A for x in r)) +def rel(A,B):return fn([[a-b for a,b in zip(r,s)]for r,s in zip(A,B)])/fn(B) +def response(X,C,K): + A=[[c/.1+k for c,k in zip(cr,kr)]for cr,kr in zip(C,K)] + T=solve(A,[float(i==0)for i in range(20)]) + flux=[] + for xi,eta in [(-1/math.sqrt(3),-1/math.sqrt(3)),(1/math.sqrt(3),-1/math.sqrt(3)),(1/math.sqrt(3),1/math.sqrt(3)),(-1/math.sqrt(3),1/math.sqrt(3))]: + N,G,J,no=geometry(X,xi,eta) + for off,s in [(0,-1),(0,0),(0,1),(2,-1),(2,0),(2,1)]: + Ls=[s*(s-1)/2,1-s*s,s*(s+1)/2];Ds=[4*(s-.5),-8*s,4*(s+.5)] + flux.append([-sum((G[i][j]*Ls[a]+no[j]*N[i]*Ds[a])*T[5*i+off+a]for i in range(4)for a in range(3))for j in range(3)]) + return T,flux +def rv(a,b):return math.sqrt(sum((x-y)**2 for x,y in zip(a,b)))/math.sqrt(sum(x*x for x in b)) +cases={"affine":[[0,0,0],[1,0,0],[1,1,0],[0,1,0]],"distorted":[[0,0,0],[2,0,0],[1.5,1,0],[0,1,0]]} +for w in [.1,.5,2]:cases["saddle"+str(w)]=[[x,y,w*x*y]for x,y in [(-1,-1),(1,-1),(1,1),(-1,1)]] +ans={} +for name,X in cases.items(): + ds=[] + for n in [3,5,7]: + C,K=mats(X,n);T,q=response(X,C,K);ds.append((C,K,T,q)) + rows=[] + for a,b in [(ds[0],ds[1]),(ds[1],ds[2])]: + rows.append([rel(a[0],b[0]),rel(a[1],b[1]),rv(a[2],b[2]),max(math.sqrt(dot([x-y for x,y in zip(q,r)],[x-y for x,y in zip(q,r)]))for q,r in zip(a[3],b[3]))/max(math.sqrt(dot(r,r))for r in b[3])]) + ans[name]=rows +print(json.dumps({"quadrature":ans},indent=2)) +# Independent uniform Q1 Neumann mode and continuum error by high-order integral. +def spatial(n): + h=1/n;lam=6*(1-math.cos(math.pi/n))/(h*h*(2+math.cos(math.pi/n))) + t=.1;af=math.exp(-lam*t);ae=math.exp(-math.pi**2*t);et=eg=0 + for i in range(n): + for s,w in gauss(12): + x=(i+(s+1)/2)*h + T=af*((1-s)*math.cos(math.pi*i/n)+(1+s)*math.cos(math.pi*(i+1)/n))/2 + G=af*(math.cos(math.pi*(i+1)/n)-math.cos(math.pi*i/n))/h + et+=(T-ae*math.cos(math.pi*x))**2*w*h/2 + eg+=(G+math.pi*ae*math.sin(math.pi*x))**2*w*h/2 + return math.sqrt(et),math.sqrt(eg) +ss=[spatial(n)for n in [8,16,32,64]] +print(json.dumps({"spatial_errors":ss,"orders":[[math.log(a/b,2)for a,b in zip(ss[i],ss[i+1])]for i in range(3)]},indent=2)) +print(json.dumps({"time_orders":[math.log(abs((1+1/n)**(-n)-math.exp(-1))/abs((1+1/(2*n))**(-2*n)-math.exp(-1)),2)for n in [16,32,64,128]]})) + + +def direct54(X): + C=zeros(20);K=zeros(20) + for xi,wx in gauss(3): + for eta,wy in gauss(3): + N,G,J,no=geometry(X,xi,eta) + for off in [0,2]: + for s,ws in gauss(3): + ls=[s*(s-1)/2,1-s*s,s*(s+1)/2] + ds=[4*(s-.5),-8*s,4*(s+.5)] + H=[0.]*20;B=zeros(20,3) + for i in range(4): + for a in range(3): + p=5*i+off+a;H[p]=N[i]*ls[a] + B[p]=[G[i][j]*ls[a]+no[j]*N[i]*ds[a]for j in range(3)] + weight=wx*wy*ws*J/4 + for i in range(20): + for j in range(20): + C[i][j]+=weight*H[i]*H[j];K[i][j]+=weight*dot(B[i],B[j]) + return C,K +for name,X in cases.items(): + C,K=mats(X,3);Cd,Kd=direct54(X) + print("54point",name,rel(C,Cd),rel(K,Kd)) + assert rel(C,Cd)<1e-12 and rel(K,Kd)<1e-12 +# Generalized eigen residual and C orthogonality, independently reassemble. +C,K=slabmat(2);Cf=[r[1:]for r in C[1:]];Kf=[r[1:]for r in K[1:]] +res=fn([[a-b*l for a,b,l in zip(kr,cr,lam)]for kr,cr in zip(mm(Kf,W),mm(Cf,W))])/fn(mm(Kf,W)) +orth=fn([[v-float(i==j)for j,v in enumerate(row)]for i,row in enumerate(mm(mm(tr(W),Cf),W))]) +print("slab eigen residual",res,"C orthogonality",orth) +assert res<1e-12 and orth<1e-12 + +def spatial_full(n,ng): + h=1/n;t=.1;steps=16*n*n + l=6*(1-math.cos(math.pi/n))/(h*h*(2+math.cos(math.pi/n))) + af=math.exp(-steps*math.log1p(l*t/steps));ase=math.exp(-l*t);ae=math.exp(-math.pi**2*t) + et=eg=etime=0 + for i in range(n): + for s,w in gauss(ng): + x=(i+(s+1)/2)*h + v=((1-s)*math.cos(math.pi*i/n)+(1+s)*math.cos(math.pi*(i+1)/n))/2 + G=af*(math.cos(math.pi*(i+1)/n)-math.cos(math.pi*i/n))/h + et+=(af*v-ae*math.cos(math.pi*x))**2*w*h/2 + eg+=(G+math.pi*ae*math.sin(math.pi*x))**2*w*h/2 + etime+=((af-ase)*v)**2*w*h/2 + return math.sqrt(et),math.sqrt(eg),math.sqrt(etime) +sbe=[] +for n in [8,16,32,64]: + a=spatial_full(n,12);b=spatial_full(n,16) + er=spatial(n)[0];ratio=a[2]/er + quadr=max(abs(x-y)/abs(x)for x,y in zip(a,b)) + assert ratio<=1/32 and quadr<=1e-8 + sbe.append(a) + print("spatial BE",n,"time/spatial",ratio,"12/16 relative change",quadr) +for i in range(3): + pt=math.log(sbe[i][0]/sbe[i+1][0],2);pg=math.log(sbe[i][1]/sbe[i+1][1],2) + assert 1.95<=pt<=2.05 and .98<=pg<=1.03 + print("spatial BE order",pt,pg) + + +``` + +## B. Exact Fraction integration·rank·jump + +```python + +from fractions import Fraction as F +def fmat(A,d):return [[F(v,d)for v in r]for r in A] +def mul(a,b): + c=[F(0)]*(len(a)+len(b)-1) + for i,x in enumerate(a): + for j,y in enumerate(b):c[i+j]+=x*y + return c +def der(a):return [i*a[i] for i in range(1,len(a))] +def integ(a):return sum(F(x,i+1)for i,x in enumerate(a)) +sh=[[1,-3,2],[0,4,-4],[0,-1,2]] +assert [[integ(mul(a,b))for b in sh]for a in sh]==fmat([[4,2,-1],[2,16,2],[-1,2,4]],30) +assert [[integ(mul(der(a),der(b)))for b in sh]for a in sh]==fmat([[7,-8,1],[-8,16,-8],[1,-8,7]],3) +M=fmat([[4,2,1,2],[2,4,2,1],[1,2,4,2],[2,1,2,4]],36) +S=fmat([[4,-1,-2,-1],[-1,4,-1,-2],[-2,-1,4,-1],[-1,-2,-1,4]],6) +fac=[([1,-1],[1,-1]),([0,1],[1,-1]),([0,1],[0,1]),([1,-1],[0,1])] +assert [[integ(mul(a,c))*integ(mul(b,d))for c,d in fac]for a,b in fac]==M +assert [[integ(mul(der(a),der(c)))*integ(mul(b,d))+integ(mul(a,c))*integ(mul(der(b),der(d)))for c,d in fac]for a,b in fac]==S +Z0=[[F(0)for _ in range(5)]for _ in range(5)] +Z1=[[F(0)for _ in range(5)]for _ in range(5)] +for off in [0,2]: + for i in range(3): + for j in range(3): + Z0[i+off][j+off]+=F([[4,2,-1],[2,16,2],[-1,2,4]][i][j],30) + Z1[i+off][j+off]+=F([[7,-8,1],[-8,16,-8],[1,-8,7]][i][j],3) +C=[[M[i//5][j//5]*Z0[i%5][j%5]for j in range(20)]for i in range(20)] +K=[[S[i//5][j//5]*Z0[i%5][j%5]+M[i//5][j//5]*Z1[i%5][j%5]for j in range(20)]for i in range(20)] +def pivots(A): + a=[r[:]for r in A];p=[] + for k in range(len(a)): + d=a[k][k];p.append(d) + for i in range(k+1,len(a)): + for j in range(k+1,len(a)):a[i][j]-=a[i][k]*a[k][j]/d + return p +assert all(x>0 for x in pivots(C)) +assert all(x>0 for x in pivots([r[:19]for r in K[:19]])) +assert all(sum(r)==0 for r in K) +assert sum(map(sum,C))==2 +cc=[[F(2),F(1)],[F(1),F(2)]];kk=[[F(1),F(-1)],[F(-1),F(1)]] +old=[F(0),F(0)] +for expected,react in [(F(-1,5),F(24,5)),(F(1,25),F(36,25))]: + new=[(2*(2*old[0]+old[1])-1)/5,F(1)] + r=[2*sum(cc[i][j]*(new[j]-old[j])for j in range(2))+sum(kk[i][j]*new[j]for j in range(2))for i in range(2)] + assert new[0]==expected and r==[0,react] + assert 3*sum(new[j]-old[j]for j in range(2))==react/2 + old=new +print("exact: C SPD 20/20; K rank19; K1=0; 1TC1=2; jump -1/5,1/25; Rc 24/5,36/25") + +``` + +## C. Time predicate와 경계 + +```python + +import math,random,struct +from fractions import Fraction as F +def input_cell(x): + v=F(x);prev=F(math.nextafter(x,-math.inf)) + nxt=math.nextafter(x,math.inf) + hi=(v+F(nxt))/2 if math.isfinite(nxt) else v+(v-prev)/2 + even=not(struct.unpack('>Q',struct.pack('>d',x))[0]&1) + return (prev+v)/2,hi,even +def cell_overlap(t,d,m): + a,b,ca=input_cell(t);c,e,cb=input_cell(d);c*=m;e*=m + lo=max(a,c);hi=min(b,e) + return lo0 and d>0 and type(inc)is int and 1<=inc<=2**32-1):return "input-error" + if d>=t:return 1,False,t + if F(t)/F(d)>inc+1:return "INC" + _,e=math.frexp(t);q=math.ldexp(t,-e)/math.ldexp(d,-e);m=round(q) + close=cell_overlap(t,d,m) + assert close==dyadic_overlap(t,d,m) + if close:n=m + else: + factor=F(2)**(-e+88);X=F(t)*factor;Y=F(d)*factor + assert X.denominator==Y.denominator==1 + X=int(X);Y=int(Y);f=math.floor(q) + if f*Y>X:f-=1 + elif (f+1)*Y<=X:f+=1 + assert f*Y<=X<(f+1)*Y + n=f if X==f*Y else f+1 + if n>inc:return "INC" + if not close:assert n==-(-F(t)//F(d)) + delta=d if close else fm(-(n-1),d,t) + assert 00):continue + t=math.nextafter(t,math.inf if i%2 else -math.inf) + if not math.isfinite(t):continue + q=t/d + if not math.isfinite(q)or q>2**32:continue + mm=round(q) + assert cell_overlap(t,d,mm)==dyadic_overlap(t,d,mm) + count+=1 +print("input-cell/fixed-dyadic independent comparisons",count) +assert grid_cell(5.351455173733616,.17262758624947144)[0]==31 +assert grid_cell(math.nextafter(.9,math.inf),.3)[0]==4 +assert grid_cell(20*math.ulp(0.),6*math.ulp(0.))[0]==3 +assert grid_cell(17*math.ulp(0.),5*math.ulp(0.))[0]==4 +print("INPUT-CELL ASSERTIONS PASSED") + +for t,d in [(1.,.1),(math.ldexp(.9,-1000),math.ldexp(.3,-1000)),(math.ldexp(.9,1000),math.ldexp(.3,1000)),(1e-323,5e-324)]: + print("input-cell scale",repr(t),repr(d),grid_cell(t,d)) +assert grid_cell(10.,3.,3)=="INC" and grid_cell(10.,3.,4)[0]==4 +assert grid_cell(float(2**32),1.,2**32-1)=="INC" +for bad in [0,-1,1.5,2**32]:assert grid_cell(1.,.1,bad)=="input-error" +for bad in [0.,-1.,math.inf,math.nan]:assert grid_cell(bad,.1)=="input-error" +print("CURRENT TIME GRID ASSERTIONS PASSED") + +``` + +## D. 입력 rounding-cell의 독립 반증 검색 + +고정 seed의 검색은 폐기된 arithmetic-only predicate와 입력 rounding-cell의 차이를 +검출한다. 이 출력은 predicate 실행 성공과 승인 시간 의미의 정합을 구별하기 위한 +반증 evidence다. 상세 판정은 numerical-review.md §8을 따른다. + +```python + +import math,random +from fractions import Fraction as F +def rad(x):return max(x-math.nextafter(x,-math.inf),math.nextafter(x,math.inf)-x)/2 +def pred(t,d,m): + q=t/d;p=m*d + B=math.nextafter(d*rad(q)+rad(p),math.inf) + return abs(float(F(t)-m*F(d)))<=B,B +rng=random.Random(41);found=[] +for z in range(20000): + d=rng.uniform(.125,.25);m=rng.randint(2,100);t=m*d + for dr in [-1,1]: + tt=math.nextafter(t,math.inf if dr>0 else -math.inf) + loT=(F(tt)+F(math.nextafter(tt,-math.inf)))/2 + hiT=(F(tt)+F(math.nextafter(tt,math.inf)))/2 + loD=(F(d)+F(math.nextafter(d,-math.inf)))/2 + hiD=(F(d)+F(math.nextafter(d,math.inf)))/2 + overlap=max(loT,m*loD)<=min(hiT,m*hiD) + ok,B=pred(tt,d,m) + if overlap and not ok: + found.append((tt,d,m,float(abs(F(tt)-m*F(d))),B));break + if found:break +print("input-cell overlap rejected by retired arithmetic predicate:",found) +t=math.nextafter(.9,math.inf);d=.3;m=3 +print("specific .9 upper:",float(F(t)-m*F(d)),rad(t)+m*rad(d),"input-cell overlap",abs(F(t)-m*F(d))<=F(rad(t))+m*F(rad(d))) + +``` diff --git a/docs/shell-transient-heat-transfer/numerical-review.md b/docs/shell-transient-heat-transfer/numerical-review.md index 3b9f5a5..ffdd13a 100644 --- a/docs/shell-transient-heat-transfer/numerical-review.md +++ b/docs/shell-transient-heat-transfer/numerical-review.md @@ -1,121 +1,495 @@ # Shell 과도 열전달 독립 수치 검토 -- 검토일: 2026-09-28. 상태: **Proposed 계획의 독립 검토 완료 / 구현 검증 미실행**. -- 판정: 제안된 약형·행렬·BE partition·반력·한쪽 heat flux는 수학적으로 일관된다. - 아래 미결정 사항 때문에 구현 착수 계약이나 제품 정확도가 승인되었다는 뜻은 아니다. -- 범위: [PRD](PRD.md)의 FEM-REQ-SHELL-THERMAL-001–018, - [정식화](formulation.md) §§1–10, [상세 계획](design-plan.md) §§1–12. -- 기준: [프로젝트 지침](../../AGENTS.md), [프로젝트 PRD](../PRD.md), - [ARCHITECTURE](../ARCHITECTURE.md), [ADR-028](../ADR.md#adr-028), - [공통 TOLERANCE](../TOLERANCE.md). 상위 계약·소스·테스트·reference는 수정하지 않았다. +- 검토일: 2026-09-29. **정식화 일관성 및 구현 전 수치 검증 설계 PASS / + solver 구현·정확도·reference·physics 검증 미실행**. +- 독립 검토 대상: formulation commit e052b6e, source 접속 범위 보완 efa5403, + 입력 rounding-cell 보완 dd409b5의 실제 diff 및 §7.1–7.2. + [PRD STH-C07](PRD.md#sth-c07--독립-검증과-수치-계약의-선행-조건)의 위임 범위에서 + 아래 식·상수·fixture·판정 기준을 production kernel 구현 전에 고정한다. +- 상태: **DONE_WITH_CONCERNS**. 고정 5점의 초기 침투층 및 warped 3×3 적분 오차를 + 정량 기록했다. 지원 방법·geometry·외부 tolerance 변경을 제안하거나 승인하지 않는다. +- 기준: [정식화](formulation.md), [승인 ADR](ADR.md), [공통 수치 기준](../TOLERANCE.md), + [개발 규칙](../standards.md), [수치 주의점](../engineering-notes.md). + 재현 계산은 [scratch appendix](numerical-oracle-scratch.md)에 있다. + Reviewer는 정식화 작성자가 아니며 production·test·reference를 수정하지 않았다. -## 1. 영향이 큰 판단과 필요한 결정 +## 1. 판정과 이전 검토의 정리 -| ID / 판정 | 근거·영향 | 필요한 조치 / owner | -| --- | --- | --- | -| NR-01 / 물리 결정 대기 | Midsurface metric을 두께 전체에 쓰는 product-domain은 실제 offset volume의 곡률 항을 생략한다. 같은 normalized ζ를 연결한 두께 단차·fold는 실제 3D 접합부가 아니다. | 메인/사용자: 이 이상화를 승인하고 지원 한계를 명시. Numerical/I/O: 곡률·왜곡·퇴화·orientation 판정과 허용 geometry 범위를 구현 전에 확정. 높은 적분 차수로 모델링 오차가 제거된다고 주장하지 않음. | -| NR-02 / 초기 상태 계약 대기 | Raw 초기 온도를 첫 BE 증분의 이전 상태로 쓰는 식은 맞지만 t=0 즉시 projection과 첫 RHS·에너지가 다르다. Reference 입력에는 step 전 `*Boundary / Set-10,11,11`과 step 안의 `Set-10,11,11,1000.`이 모두 있다. | 메인/사용자: 첫 증분 jump 의미 결정. I/O: 생략 초기값, 생략 boundary 값, pre-step boundary와 step 재정의의 의미를 명시. 단일 transient step 지원을 이유로 pre-step 데이터를 무시하지 않음. | -| NR-03 / 공개 입력 결정 대기 | Edge kernel의 면적당 입력과 `q_in hL`은 일관되지만 DS4 `*DFLUX, S1…S4` 대응은 확인되지 않았다. 기존 reference에는 edge 하중이 없다. | I/O/메인: 공개 syntax·단위·대상 edge·중복 입력 의미 확인. 입증되지 않은 keyword 지원을 선언하거나 line-rate에 h를 중복 곱하지 않음. | -| NR-04 / 구현 전 수치 계약 부족 | 3×3 표면 적분은 일반 warped Q4에서 exact가 아니다. Consistent C와 P2 두께 보간은 maximum principle을 보장하지 않는다. 고정 5점은 초기 두께 경계층을 임의로 해상할 수 없다. | Numerical owner: 아래 geometry·적분·수렴·내부 residual 판정을 고정하고 독립 사례로 검증. 실패 뒤 clamp/lumping/인공 확산 또는 tolerance 완화 금지. 방법 변경은 새 검토 대상. | -| NR-05 / 문서 정합 보완 반영 | 최초 계획의 24개 integration Gauss 위치와 정식화의 54점 직접 적분·24개 trace recovery가 달랐다. 최신 계획에서 구분됨. | 유지할 계약: production 3×3 surface×exact Z, 직접 검산 54점, affine 검산 24 Gauss점, 공식 recovery 24 trace점. 숫자 24가 같다는 이유로 위치를 대체하지 않음. | -| NR-06 / 검증 범위 보완 반영 | 직렬 strip 저항식은 단면 전체 end 온도 조건을 전제로 한다. 내부 station 지정 온도와 all-constrained는 현재 공개 입력 범위 밖이다. | 최신 계획처럼 exact strip 및 all-constrained를 assembly/algebra fixture로 분류. CLI bottom/top 지정 사례는 해당 경계층을 포함한 oracle로 검증. | +2026-09-28 Proposed 검토는 당시 evidence다. NR-01의 product-domain 및 normalized +station 접속, NR-02의 raw 초기값과 BC 상속/override는 현재 승인 계약으로 해결되었다. +NR-03의 edge 입력 제안은 **011 전체 연기**로 대체된다. 과거 edge weight는 이번 +구현·검증 요구가 아니다. 종료 시간의 정수배 제한 제안도 폐기되었다. -NR-01의 물리 지원 제한, NR-02의 시작 의미, 종료 시간이 Δt의 정수 배가 아니면 거부하는 -제안은 사용자 요구에 영향을 주는 결정이다. Z 계수·scatter 순서·구간 Jacobian·동일 -물리량의 검증 metric 고정은 합의된 물리를 구현하기 위한 상세 계약이다. 영속 flux -좌표계·trace identity·t=0 validity는 I/O 계약이며 이 보고서가 대신 승인하지 않는다. - -## 2. 독립 수식·정확 산술 검산 - -Python 표준 라이브러리 `fractions.Fraction`을 명령 내 scratch로 사용했다. 별도 계산 파일, -production/test 변경 또는 solver 실행은 없었다. 아래 입력과 결과로 검산을 재현할 수 있다. - -| 항목 | 검산과 결과 | +| 항목 | 독립 판단 / 현재 경계 | | --- | --- | -| DOF/shape | Node-major `p=5(i−1)+a`, a=1…5. 각 길이 1인 두께 구간의 좌표 u∈[0,1]에서 shape `[1−3u+2u², 4u−4u², −u+2u²]`를 곱하여 항별 정확 적분. | -| 두께 질량 | `Z0_seg=(1/30)[[4,2,−1],[2,16,2],[−1,2,4]]` 재현. 조립 뒤 `Z0·1=(1,4,2,4,1)/6`, 합 2. | -| 두께 전도 | 위 다항식을 u로 미분하여 적분하면 `Z1_seg=(1/3)[[7,−8,1],[−8,16,−8],[1,−8,7]]`. 조립 뒤 `Z1·1=0` 정확 성립. | -| 체적/계수 | `dz=(h/2)dζ`, `∂z=(2/h)∂ζ`이므로 `C=ρc(h/2)M⊗Z0`, `K=k[(h/2)S⊗Z0+(2/h)M⊗Z1]`. 접선과 normal의 직교성으로 교차항 0. C는 energy/temperature, K는 power/temperature. | -| 평판의 rank/에너지 | 단위 정사각형, `h=2`, `k=ρc=1`의 수치 단위 fixture. `M=(1/36)[[4,2,1,2],[2,4,2,1],[1,2,4,2],[2,1,2,4]]`, `S=(1/6)[[4,−1,−2,−1],[−1,4,−1,−2],[−2,−1,4,−1],[−1,−2,−1,4]]`. C의 정확 LDL pivot 20개 모두 양수, K의 leading 19×19 pivot 모두 양수, `K·1=0`, `1ᵀC1=2` 정확 성립. 즉 C SPD, K PSD/rank 19. | -| Edge weight | 한 endpoint당 `q_in hL(1,4,2,4,1)/24`; 두 endpoint의 합은 `q_in hL`. 두께 station에 1/5씩 나누면 잘못된 일관 하중이다. | -| 적분 차수 | Affine Q4 면내 M/S는 2×2 exact. 두께 mass·면내 전도는 4차여서 구간별 3점, 두께 전도는 2차. Production 3×3×2×3=54점은 exact Z 분리식과 같은 quadrature를 표현한다. 실제 54점 코드의 동등성은 미검증. | +| 약형·분리 행렬·부호 | Q4×두 C0 P2, 실제 midsurface metric, exact Z, consistent C와 BE는 일관된다. 단위·node-major 순서·법선 방향을 확인했다. | +| Whole patch | affine normal image의 separating vector certificate 및 bilinear midpoint identity의 injectivity 유도가 성립한다. 적분점만 검사하는 것은 불충분하다. | +| Source 접속 | 전역 node-edge scan은 disconnected/contact 범위를 침범할 수 있었다. efa5403의 source endpoint를 공유하는 coarse/partial edge 후보 제한을 독립 확인했다. | +| 초기 jump | raw t=0 보존 및 첫 RHS의 C_fc T_c_old 항이 필요하다. Slab의 δ→0 자유 초기 극한은 분석용 oracle이며 raw 초기 dataset을 바꾸는 projection 절차가 아니다. | +| 시간 | 표현오차 수준 endpoint-close와 실제 shortened δ를 구별하는 predicate를 검산했다. 저장 시각의 subtraction과 연산 δ를 혼동하지 않는다. | +| Recovery | 4 면내점×6 두께 trace의 global q와 R=[e1,n×e1,n], 0−/0+ 분리, full residual 반력 부호가 일관된다. t=0 flux/reaction은 없다. | +| 남는 위험 | consistent C의 초기 undershoot, 고정 두께 해상도, metric conditioning, 3×3 유한 적분 오차. Clamp/lumping/차수 변경/지원 축소로 숨길 수 없다. | -유효 geometry와 양의 적분 weight에서는 `vᵀCv>0`, `vᵀKv≥0`의 에너지 형태가 유지된다. -전역 K의 상수 null mode는 **연결된 thermal component마다 하나**이며, 분리된 mesh를 -항상 rank N−1로 검사하면 안 된다. 모든 활성 DOF가 양의 체적 요소에 연결되면 -`A=C/Δt+K`는 구속 없이도 SPD다. 극단적인 h 또는 Δt는 conditioning을 악화시킬 수 -있으므로 factorization 성공만으로 정확도를 판단하지 않는다. +PASS는 **설계의 수학적 일관성과 독립 사전 oracle**에 한정한다. +실제 MSVC solver·HDF5·외부 reference 결과가 통과했다는 뜻은 아니다. -첫 증분 partition의 독립 반례 fixture는 다음과 같다. 일관된 단위로 -`C=[[2,1],[1,2]]`, `K=[[1,−1],[−1,1]]`, `Δt=1/2`, `F=0`, -`T_raw⁰=(0,0)`, 두 번째 DOF의 지정 온도 `g=1`을 둔다. +## 2. Geometry와 mesh 판정 -- 올바른 `A=[[5,1],[1,5]]`; `T_f¹=−1/5`, `T_f²=1/25`. -- `R_c¹=24/5`, `R_c²=36/25`; `E=1ᵀCT`에서 두 증분 모두 `ΔE=Δt R_c` 정확 성립. -- 최초 RHS를 일정 g의 이후 증분 식으로 바꾸면 `T_f¹=+1/5`가 되어 검출된다. - 이 fixture는 signed 온도와 consistent C의 undershoot가 에너지 안정성과 모순되지 - 않음을 보이며 실제 20 DOF shell 해나 외부 reference 통과 증거는 아니다. +### 2.1 Whole-patch certificate와 uncertainty -`r=CΔT/Δt+KT−F`의 constrained 성분은 모델에 공급되는 양의 반력 열률이다. -전역 수지는 `ΔE=Δt(1ᵀF+1_cᵀR_c+1_fᵀr_f)`이며, free residual을 포함한 대수적 -동일성 확인과 free residual 자체가 작은지의 확인을 분리한다. `q=−kBT` 경계 적분은 -일반적으로 이 반력과 같지 않다. 입력 없는 단열계의 열에너지 보존과 BE의 C-norm -소산도 서로 다른 검증이다. 초기 t=0 반력은 미정의이고 material history update는 없다. +정식화 §2.1의 r=a+bξ+cη+dξη, w=w0+uξ+vη를 사용한다. +Convex set W={w(ξ,η)}의 원점에 가장 가까운 w*는 w*·w≥||w*||²를 만족한다. +따라서 네 corner에 대한 **같은** 고정 double 벡터 s의 양의 interval dot product가 +전체 patch를 증명한다. 후보 최소점 계산이 부정확해도 최종 certificate가 양수이면 +안전하다. 후보 탐색 실패를 수학적 퇴화의 증명으로 보고하면 안 된다. +r(p)−r(q)=Dr((p+q)/2)(p−q)이므로 전체 patch에서 rank(Dr)=2이면 단사다. +이 증명은 bilinear map 전용이며 일반 고차 curved element로 확대하지 않는다. -## 3. 구현 후 반증 가능한 검증 계획 +**고정 규칙:** binary64, u=2^-53, epsilon=2^-52, round-to-nearest, +gradual underflow; 입력 double은 singleton interval이다. 각 사칙연산 endpoint를 +양쪽 nextafter로 한 번씩 넓힌다. 0을 포함한 divisor는 사용하지 않는다. +sqrt는 비음수 endpoint를 계산한 뒤 같은 방식으로 넓힌다. 재척도는 exact +power-of-two이며 손실된 subnormal도 enclosure에 포함한다. Fast-math/FTZ와 +증명되지 않은 재결합은 이 bound의 전제가 아니다. -| 사례·수준 | 독립 oracle와 판정 의도 | +모든 lower(s·w_corner)>0이면 +J_lower=min lower(s·w_corner)/upper(sqrt(s·s))>0를 보존한다. +증명되지 않으면 다른 유한 candidate를 시도하고 모두 실패하면 uncertain geometry로 +fail closed한다. 정확히 중복된 좌표·source node, 비유한 값, 증명된 w=0은 invalid다. +고정 Jacobian epsilon·warp/aspect cutoff는 **0개**다. Certificate 뒤에도 각 사용점의 +metric inverse/gradient가 유한하고 interval에서 0 divisor가 배제되는지 확인한다. +κ(g)는 진단과 오차 전파에 쓰며 지원 cutoff가 아니다. + +### 2.2 Source topology와 T-edge uncertainty + +Edge incidence는 1 또는 반대 순서의 2; 3 이상·같은 방향·중복 patch는 거부한다. +각 vertex link는 하나의 path 또는 cycle이며 끊어진 fan·branch·중복 link는 거부한다. + +T-edge 후보는 **같은 instance에서 source endpoint를 공유한** +(p,q)와 (p,v) 또는 (v,q)이고 p,q,v는 서로 다른 node다. Coarse 방향을 양쪽으로 +시도한다. a=q−p, b=v−p의 interval에서 |a_j| lower가 가장 큰 성분 j를 택하고 +λ=b_j/a_j, 잔여 성분 b_i−λa_i를 계산한다. + +- 잔여 성분 interval 하나가 0을 제외하거나 λ가 [0,1] 밖이면 비접합이다. +- 공선 및 0<λ<1이 증명되면 nonconforming T-edge다. +- 공선/열린 구간 포함이 불확실하면 uncertain 접속으로 실패한다. 정확한 끝점 일치는 + 별도 중복 좌표 진단이며 열린 구간으로 clamp하지 않는다. +- Interval이 0을 포함한다고 exact 공선을 주장하지 않는다. Exact double-rational + 검산은 선택 사항이며 새 arbitrary-precision runtime dependency를 요구하지 않는다. + Interval fail-closed만으로 계약을 충족한다. + +모든 node 대 모든 edge scan, 좌표 welding, source 후보가 없는 component 접촉 검출은 +하지 않는다. 같은 component의 멀리 떨어진 부분끼리의 접촉도 대상이 아니다. +이는 global contact 검증을 완료했다는 뜻이 아니다. + +### 2.3 고정 geometry fixture군 + +| 군 | 좌표/변환과 필수 판정 | | --- | --- | -| Shape/행렬, kernel | 위 exact M/S/Z를 기대값으로 사용. Node-major Kronecker 순서, symmetry, 상수 null, SPD, 양의 에너지, 직접 54점/분리식 동등성 확인. h를 바꿔 C 및 면내 K의 h scaling과 두께 K의 1/h scaling을 각각 검사. | -| 온도 patch/방향, kernel | 평판에서 `T=A+b·r+c z`를 주고 `q=−k(b+c n)` 확인(b는 평판 접선). 임의 rigid rotation에서 T 불변, q 회전. Connectivity 반전 시 station a↔6−a와 면 조건도 함께 반전. | -| 한쪽 trace, recovery | 면내 상수, 두께 값 `(0,0,0,1,2)`의 임의 kernel state는 `q3(0−)=0`, `q3(0+)=−4k/h`. 24개 trace 중 0−/0+가 별도 identity로 저장되는지 확인. 이 임의 state는 공개 초기조건 fixture가 아님. | -| 열용량/하중/수지, CLI+assembly | 단열·균일 초기 온도는 불변. Edge만 주면 총 입력 q_in hL; 집중 입력은 P 그대로. Prescribed DOF의 P도 full F에 남아 reaction을 바꾸는지 확인. 반력 부호, free residual, ΔE와 BE 소산을 각각 확인. | -| BE, algebra | 위 두 증분 fixture, no-constrained 및 generic all-constrained. 이전 constrained state를 raw 초기값으로 유지하고 실패 candidate가 이전 state/output을 덮어쓰지 않는지 확인. | -| 면내 해석해, CLI | 단열 직사각판 `T=T0+A cos(πx/L)exp[−α(π/L)²t]`, `α=k/(ρc)`. 초기값은 station 전체 동일. 두께 균일해가 유지되는지와 면내 temperature/flux 오차 확인. | -| 시간 수렴, 독립 semi-discrete oracle | 같은 고정 mesh의 generalized eigenpair `Kv=λCv`로 `T(t)=exp(−λt)v`, BE 증폭률 `(1+λΔt)^−n` 대조. dt, dt/2, dt/4에서 1차 오차를 확인해 공간오차와 분리. Eigenvector 초기값은 algebra fixture로 한정. | -| 두께 slab, CLI | x=z+h/2, bottom Ts, top 단열, 초기 Ti. `θ=(T−Ts)/(Ti−Ts)=Σ[m≥0]4/[(2m+1)π] sin[(2m+1)πx/(2h)] exp[−α((2m+1)π/(2h))²t]`. 양의 시간에서 5 station T와 미분한 q, 전체 energy 및 reaction 비교. Series 절단 오차를 별도 제한. | -| Mesh/두께 해상도 | Smooth affine mesh에서 면내 temperature L2의 기대 2차, gradient L2의 기대 1차를 확인. Slab은 고정 5점이므로 dt→0의 두께 공간오차 plateau와 열침투 깊이 √(αt)를 보고. 면내 refinement로 두께 오차가 사라진다고 판정하지 않음. | -| 재료·fold·두께 접속 | 직렬 저항 ΣL/(khw)는 full-end prescribed assembly fixture. CLI에서는 bottom/top 조건의 실제 oracle을 별도로 사용. Fold 등온은 q=0; 비균일 경우 global q의 동일성을 강요하지 않고 공유 DOF에 조립되는 열률 수지를 검사. | -| Warped/왜곡 geometry | 허용될 geometry군에 대해 3×3→5×5→7×7에서 K/C 및 대표 T/q의 변화 확인. Geometry mesh 오차, quadrature 오차, offset-volume 생략 오차를 분리. Degenerate·자기 교차·orientation 불일치는 음성 사례. | +| Affine | (0,0,0),(1,0,0),(1,1,0),(0,1,0); 평행사변형 (0,0,0),(2,0,0),(3,1,0),(1,1,0). Certificate PASS. | +| 평면 왜곡 | (0,0,0),(2,0,0),(1.5,1,0),(0,1,0). J>0; M exact, S의 3×3 exactness 주장 금지. | +| Saddle | r=(ξ,η,aξη), a=.1,.5,2; 자연 corner를 node로 사용. w=(-aη,-aξ,1), s=(0,0,1), J_lower=1의 exact 근거. 모두 지원. | +| 긴 aspect | rectangle 1×2^-40; 유효 certificate 필요. 조건수와 metric 계산을 별도 기록하며 aspect만으로 거부 금지. | +| Scale/이동/회전 | 유효군의 2^-200,1,2^200배; unit patch에 (2^20,-2^20,2^20) 이동. Signed axis permutation 및 z축 3-4-5 회전. 입력 반올림을 oracle interval에 포함. | +| 적분점 검사의 반례 | r=(ξ,(ξ−1/4)η,0). 3×3 모든 점의 J는 비영이지만 ξ=1/4에서 J=0. Whole-patch 거부. | +| Invalid | bow-tie 순서, 반복 source/좌표, collinear 네 node, NaN/Inf. 자동 재순서/normal flip 금지. | +| Uncertain | affine b=(1,1,0),c=(1,1+2^-52,0): cancellation interval이 0을 포함할 수 있다. 증명된 certificate 또는 uncertain 실패만 허용; 임의 epsilon/형상 보정 금지. | +| 접속 양성 | 반대 순서 edge 공유 Q4 두 개, 90도 fold, h=1/2 두께 단차, k=1/10 재료 접속; disconnected 두 patch. | +| 음성/대조 | edge 3 incidence, 동일 directed edge, node 하나만 공유한 두 fan; coarse(p,q)/partial(p,v),v=(p+q)/2. 각 음성의 별도 source identity/instance 좌표 접촉 대조군은 자동접속·T-edge 거부가 없어야 함. | -내부 roundoff 판정은 exact oracle를 기준으로 행렬 norm·연산 수·machine epsilon과 -conditioning을 반영하여 구현 전에 고정한다. 예를 들어 free residual은 -`||r_f||/(||A_ff|| ||T_f||+||b_f||)`, 수지 defect는 저장·입력·반력 항의 절댓값 합으로 -정규화하되 분모 0은 별도 처리한다. Zero floor `max(1,scale)`는 추가하지 않는다. -Roundoff 한계, geometry margin, quadrature 변화 한계와 관찰 order의 허용 구간은 아직 -확정되지 않았다. 위 기대 order 자체를 임의 5% 통과 기준으로 바꾸거나 공통 외부 -tolerance를 이 검사들에 전용하지 않는다. 이 미결정 때문에 검증 실행 준비 완료 판정은 보류한다. +## 3. 내부 roundoff, 행렬 및 solve 기준 -## 4. 기존 reference의 비교 계약과 관찰 +### 3.1 공통 오차 언어 -Read-only 원본: `reference/shell-thermal/shell-thermal.inp`, -`reference/shell-thermal/shell-thermal NodalTemperature.csv`. -이번 검토에서 실제로 읽은 위치는 `C:/git/FESADev/reference/shell-thermal`이며 -현 worktree로 복사하지 않았다. 별도 `reference-model.md` 없이 이 절에 필요한 계약을 둔다. +외부 5%/1%와 무관하게 γ_m=mu/(1−mu),u=2^-53이며 mu<1을 요구한다. +행렬/shape는 Frobenius norm, solve는 ∞ norm, 물리 integral은 아래 L2 norm이다. +분모가 0이면 numerator도 정확히 0일 때만 통과한다. max(1,scale)은 금지다. -- 제안 family: `(shell-thermal, Step-1, t=37, bottom temperature, NT11, blocking)`. - 온도 단위는 원본 입력의 일관된 단위이며 scalar이므로 좌표 회전과 무관하다. -- Row identity: `(instance=Part-1-1, source node label, source DOF=11)`와 step/time. - CSV에 instance 열이 없으므로 입력의 단일 instance와 source label 집합을 검증한 뒤 복원. - HDF5는 모든 station을 보존하되 비교 view는 bottom만 명시적으로 선택한다. -- Snapshot 정적 검사: 729행·729개 고유 node label, 모두 finite, 단일 Step-1/t37 frame. - min `−15.5742216`, max `1000`; `S=1000`, near-zero band `B=10`. - 음수 전부가 near-zero는 아니다. `−15.5742216`은 상대오차 5% 대상이다. -- [공통 정책](../TOLERANCE.md): `|r|≤B`이면 `|f−r|≤B`, 나머지는 상대오차≤0.05, - 모든 행의 relative RMS≤0.01을 함께 만족해야 한다. S=0은 exact-zero 규칙 적용. - Reference의 음수·미소값을 clamp하거나 단위/온도 영점을 이동하여 scale을 바꾸지 않는다. -- Missing/extra/duplicate/nonfinite, frame/instance/station mismatch를 tolerance 전에 거부. - 현재 수행한 CSV 자체 검사는 FESA HDF5와의 일대일 대응 검증을 대신하지 않는다. +P_ij는 계산한 signed 합의 절댓값이 아니라 **각 quadrature/분리 항의 절댓값 합**이다. +Cancellation 때문에 P와 |X|를 바꿔 쓰지 않는다. 단위 affine fixture의 closed-form +rational oracle에 대한 고정 roundoff bound는 다음과 같다. -원본 NT11 한 frame은 내부 4 station, flux, reaction, edge load, 다중 재료/두께 또는 -전체 시간 이력을 검증하지 못한다. 초기값·pre-step 경계·jump·consistent C·시간 적분 -차이가 reference 비교에 미치는 영향은 아직 측정되지 않았다. 비교 실패는 이 차이들을 -분리 조사할 근거이며 FESA가 Abaqus 내부 알고리즘과 같아야 한다는 요구가 아니다. +| 검사 | 고정 bound / 연산 수 근거 | +| --- | --- | +| M,S,Z 및 분리 C,K | \|\|X_hat−X_exact\|\|F ≤γ_4096\|\|P_X\|\|F. 점당 shape·metric·gradient·weight·entry scalar path128 이하×9+분리/scatter64 이하를 포괄한다. | +| 직접54점 대 분리 | \|\|X_54−X_sep\|\|F ≤γ_8192\|\|P_54\|\|F+γ_4096\|\|P_sep\|\|F. 54×128+reduction<8192. | +| Affine gradient·H·R | component 절댓값 전파 scale P에 error norm≤γ_2048\|\|P\|\|. Partition sum은 P=Σ\|H\|, gradient sum은 P=Σ\|\|B_p\|\|; 비영 patch는 exact gradient 크기도 포함. | +| Symmetry | \|\|X−Xᵀ\|\|F≤2E_X. E_X는 위 bound. 강제 대칭화로 검사를 통과시키지 않는다. | +| 상수 null | \|\|K1\|\|2≤sqrt(20)E_K; \|1ᵀC1−2\|≤20E_C. Fraction 검산은 exact equality. | +| h scaling | h=1/2,1,2, k/rho/c=1; C와 면내 K는 h, normal K는 1/h. 두 계산의 E를 해당 배율로 더한다. | -## 5. 수행 범위와 잔여 위험 +Operation budget은 단순 loop의 최장 path 상한이다. 실제 경로가 이보다 길면 결과에 맞춰 +상수를 늘리지 않고 경로/interval을 다시 검토한다. 일반 warped/translated geometry에는 +affine bound를 그대로 쓰지 않는다. §2 interval을 같은 quadrature에 전파하여 X^I를 +얻고 E_X=||max(|lower(X^I)−X_mid|,|upper(X^I)−X_mid|)||F를 사용한다. +각 point의 κ2(g)와 κ2(g)γ_64를 기록한다. Inverse perturbation 상대 증폭은 +κ2(g)γ_64/(1−κ2(g)γ_64)이고 분모가 양수가 아니면 이 상계는 불확실하다. +직접 interval inverse가 유한 enclosure를 증명하면 사용하고 아니면 uncertain이다. +κ를 geometry 지원 조건으로 바꾸지 않는다. Roundoff와 quadrature 절단 차이는 별개다. -수행: 관련 계약의 정적 대조, 위 exact 유리수 행렬/두 증분 검산, 원본 입력·CSV의 -read-only 정적 확인, 문서 정합 확인. 미수행: solver 구현·build/CTest, FESA HDF5 생성, -외부 reference 수치 비교, 곡면/접합부/수렴 실험, 실제 결과의 physics evaluation. -초기·geometry·edge syntax·내부 수치 판정·영속 결과 계약을 확정하고 위 검증을 실행해야 -한다. 본 보고서의 판정은 **정식화 초안의 독립 일관성 검토**로 한정한다. +### 3.2 Exact fixture와 rank + +단위 square,h=2,k=rho=c=1: +M=1/36 [[4,2,1,2],[2,4,2,1],[1,2,4,2],[2,1,2,4]], +S=1/6 [[4,-1,-2,-1],[-1,4,-1,-2],[-2,-1,4,-1],[-1,-2,-1,4]]. + +길이1 ζ 구간의 Z0=1/30 [[4,2,-1],[2,16,2],[-1,2,4]], +Z1=1/3 [[7,-8,1],[-8,16,-8],[1,-8,7]]을 station(1,2,3),(3,4,5)에 +조립한다. C=M⊗Z0,K=S⊗Z0+M⊗Z1, node-major 순서다. +정확 유리수 LDL에서 C의20 pivot>0,K의 leading19 pivot>0,K1=0,1ᵀC1=2를 +확인했다. C SPD,K PSD/rank19가 독립 증명된다. + +Floating 결과는 matrix bound 외에 C의 λ_min>E_C, K를 1⊥에 제한한 +λ_min>E_K와 null eigenvalue 크기≤E_K를 확인한다. Eigen 연산 자체의 backward +error enclosure도 더한다. Separation이 불가능하면 rank는 inconclusive이며 +clamp로 PASS를 만들지 않는다. 전역 K는 component마다 null mode 하나다. + +Gradient fixture는 단위 affine 평판의 접선 b에 대해 T=A+b·r+cz, +A=2,b=(3,-2,0),c=4로 q=−k(b+cn)이다. Warped patch의 임의 global b에 +같은 식을 강제하지 않으며 그 경우 실제 reciprocal tangent의 projection을 사용한다. +Constant T=2, 면내/두께선형도 각각 검사한다. 회전 후 q는 같은 회전, +normal 반전 후 station a↔6−a와 조건까지 재매핑한다. RᵀR=I,detR=+1. +Trace state(0,0,0,1,2)는 q3(0−)=0,q3(0+)=−4k/h이며 평균하지 않는다. + +### 3.3 Free solve와 conditioning + +n=free DOF 수,A=A_ff,b=같은 δ와 full old state로 만든 RHS다. +Full residual의 free 성분 및 독립 A T_f−b를 확인한다. + +β=||r_f||∞/(||A||∞||T_f||∞+||b||∞), +β_limit=nγ_(n+1)+γ_(2n+4)+γ_(4n+8). + +첫 항은 Cholesky의 |L||L|ᵀ bound와 || |L||L|ᵀ ||∞≤n||A||∞에서, +나머지는 triangular solves 및 residual 평가의 dot-product bound에서 나온다. +Backend가 전제를 충족하지 못하면 더 엄밀한 backward-error evidence가 필요하며 +β_limit를 늘리지 않는다. n=0이면 free solve/β를 생략하고 full reaction·수지를 검사한다. + +κ∞(A)β_limit<1일 때 forward relative-error bound는 +κ∞(A)β_limit/(1−κ∞(A)β_limit). +κ는 작은 fixture의 독립 dense inverse norm 또는 검증된 condition estimator의 +상계로 구한다. 이 상계가 유효하지 않으면 backward PASS를 온도 정확도 PASS로 +바꾸지 않는다. 참값 비교에도 이 forward bound와 coefficient interval을 전파한다. + +## 4. 에너지, 반력과 BE 소산 + +C=[[2,1],[1,2]],K=[[1,-1],[-1,1]],δ=1/2,Q=0,rawT0=(0,0),g2=1에서 +T_f=(-1/5,1/25),Rc=(24/5,36/25)를 exact 확인했다. Raw constrained 항을 빼거나 +먼저 g로 projection하면 실패한다. Signed undershoot를 삭제하면 안 된다. +이 fixture와 all-constrained 검사는 algebra다. Full-end strip ΣL/(khw)도 assembly +전용이며 공개 내부 station BC 지원을 뜻하지 않는다. + +dT=Tnew−Told,ΔE=1ᵀC dT,Rc=(C dT/δ+K Tnew−Q)_c. +온도 영점을 뺀 두 큰 U를 차분하는 대신 ΔE를 직접 계산해도 같은 물리량이다. + +**물리 수지 defect** +D_E=|ΔE−δ(ΣQ+ΣRc)|, +S_E=Σ_ij|C_ij dT_j|+δΣ_i|Q_i|+δΣ_c|Rc_i|. + +m=전체 DOF 수,D_f=||A_ff||∞||T_f||∞+||b_f||∞라 하면 + +B_E=δnβ_limit D_f+δ||1ᵀK||1||Tnew||∞+γ_(4m+12)S_E ++δγ_(2m+6)(|| |C||dT| ||1/δ+|| |K||Tnew| ||1+||Q||1) ++γ_1Σ_ij|C_ij|(|Tnew_j|+|Told_j|). + +두 번째 줄은 full residual 평가의 dot-product/division/subtraction 오차이며, +마지막 항은 온도 subtraction 오차다. D_E≤B_E와 free β≤β_limit를 **각각** +만족해야 한다. K 상수 mode 항은 §3 null bound도 만족해야 하므로 K 오류를 면제하지 +않는다. 보고 metric D_E/S_E에서 S_E=0이면 exact-zero 분기다. Actual free residual을 +수지에 다시 넣어 0을 만드는 대수 항등식은 물리 평형 판정을 대신하지 않는다. + +단열·무입력·무구속의 별도 C-norm 소산 defect는 + +D_C=|0.5(TnewᵀCTnew−ToldᵀCTold+dTᵀCdT)+δTnewᵀKTnew|. + +S_C는 앞의 세 C quadratic form의 **성분별 절댓값 합**에 각각 0.5를 곱하고, +마지막 K quadratic form의 성분별 절댓값 합에 δ를 곱하여 모두 더한 값이다. +B_C=δ||Tnew||1β_limit D_f+γ_(6m+20)S_C ++0.5||dT||2||C−Cᵀ||2||Told||2 ++γ_1(|Tnew|+|Told|)ᵀ|C||dT|. + +Residual work, 합/subtraction, symmetry defect, dT 오차 전파를 구분한 식이다. +D_C≤B_C,||Tnew||C²≤||Told||C²+2B_C를 검사하고 S_C=0은 exact-zero다. +이 제곱 norm은 실제 열에너지와 단위가 다르다. + +P/−P,constrained DOF의 P,disconnected component를 별도 실행한다. 집중 입력은 +signed power 그대로 Q에 가산한다. 단열 균일 T=2 보존,no-BC transient,component별 +수지,정상상태 ΣQ+ΣRc=0을 검사한다. Raw flux 경계 적분과 FE 반력의 exact equality를 +강제하지 않는다. + +## 5. Warped quadrature의 사전 진단 + +§2 affine/distorted/saddle a=.1,.5,2에서 k=rho=c=h=1. +각 3/5/7 surface rule과 exact Z로 C_r,K_r를 독립 계산한다. +대표 응답은 old=0,Q=e_(node1,station11),δ=.1,noBC의 한 BE 증분이다. + +d_X(r,s)=||X_r−X_s||F/||X_s||F(X=C,K), +d_T=||T_r−T_s||2/||T_s||2, +d_q=max_trace||q_r−q_s||2/max_trace||q_s||2. + +| fixture | d_C(3,5)/(5,7) | d_K(3,5)/(5,7) | d_T(3,5)/(5,7) | d_q(3,5)/(5,7) | +| --- | --- | --- | --- | --- | +| distorted | roundoff/roundoff | 2.189065e-7/5.970737e-12 | 2.385149e-6/6.505495e-11 | 1.070060e-6/2.918586e-11 | +| saddle .1 | 3.944257e-7/7.758135e-13 | 3.862043e-7/7.500898e-13 | 7.916236e-7/1.286406e-12 | 9.283663e-7/1.601583e-12 | +| saddle .5 | 1.667355e-4/1.431780e-7 | 1.643476e-4/1.408553e-7 | 3.301794e-4/2.150240e-7 | 3.703751e-4/2.649276e-7 | +| saddle 2 | 5.631484e-3/2.213607e-4 | 5.618988e-3/2.217451e-4 | 6.690109e-3/1.033767e-4 | 5.483294e-3/8.641481e-5 | + +표의 반올림 값을 tolerance로 쓰지 않는다. 판정은 (1) production3점이 독립3점 +oracle와 §3 bound로 일치,(2) 각 fixture d(5,7)≤d(3,5)+양쪽 roundoff bound, +(3) affine 모든 차이는 roundoff 이내다. Appendix에서 재계산한 metric과 interval을 +함께 보고한다. d가 roundoff 이하이면 ratio/order를 만들지 않는다. 독립54점/분리 +C,K의 scratch 상대 norm 차이는 최대7.73e-16이었다. + +이는 **선정 fixture의 적분 수렴 진단**이다. 7점을 exact해로 보거나 모든 허용 +geometry에 같은 d 상한을 요구하지 않는다. 일반 모델의 진단이 비단조이면 조사·보고하며 +지원 geometry cutoff,production 차수 변경/adaptive quadrature로 자동 교체하지 않는다. + +## 6. 시간·면내 수렴 + +### 6.1 고정 spatial eigen oracle의 시간1차 + +Affine 단열 unit square,균일8×8 mesh,h=1,k=rho=c=1, +station 동일 v_ij=cos(πi/8)를 쓴다. Consistent Q1의 정확한 discrete eigenvector이며 +λ_h=6(1−cos(π/8))/[(1/8)²(2+cos(π/8))]. +독립 assembly로 Kv−λCv 및 C normalization을 §3 bound로 확인한다. +τ=1/λ_h,δ=τ/N,N=16,32,64,128,256,INC≥256. + +T_exact=e^-1v,T_BE=(1+1/N)^-Nv, +e_N=||T_BE−T_exact||C/||v||C,p=log2(e_N/e_2N). +모든 pair의 p 허용구간은 **[.97,1.02]**. +독립 기대값은 .9817124,.9907338,.9953356,.9976599. +log(1+x)의 전개에서 e_N=e^-1/(2N)+O(N^-2)이며 N≥16의 독립 값으로 구간을 +고정했다. Roundoff propagation이 e_N/1000을 넘으면 order는 inconclusive다. +공간오차와 시간오차를 분리하는 고정 spatial oracle다. + +### 6.2 Fourier rectangle의 면내2차/gradient1차 + +Unit square,h=1,k=rho=c=1,noBC/noQ,T0=cos(πx)를 모든 station에 동일하게 둔다. +τ=.1에서 연속 T=cos(πx)exp(−π²t),∇T=(-πsin(πx)exp(−π²t),0,0)와 비교한다. +Affine n×n,n=8,16,32,64; norm은 product-volume temperature 및 global gradient L2. +평가 적분은 cell별12×12 Gauss×exact thickness다. 12→16점 oracle 변화가 +해당 오차의 1e-8 이하인지 확인한다(관찰 order 영향1e-7 미만의 적분 budget). + +기본 공간 oracle는 독립 mode의 exp(-λ_hτ)다. 실제 BE run은 +N=16n²,δ=τ/N,INC=N으로 고정하고 같은 mode의 BE-exponential temperature L2 차이가 +continuum 온도 오차의 **1/32 이하**임을 확인한다. 시간오차를 공간수렴으로 오인하지 +않으며 동일 mesh의 finite-N modal oracle도 함께 비교한다. + +| n | Semidiscrete temperature L2 | Gradient L2 | +| ---: | ---: | ---: | +| 8 | 6.832505333e-3 | 9.419553867e-2 | +| 16 | 1.720379027e-3 | 4.697214677e-2 | +| 32 | 4.308592308e-4 | 2.346997504e-2 | +| 64 | 1.077625665e-4 | 1.173296293e-2 | + +모든 pair에서 T p는 **[1.95,2.05]**,gradient p는 **[.98,1.03]**. +독립 T order=(1.98969,1.99744,1.99936),gradient=(1.00385,1.00099,1.00025). +Q1 smooth interpolation O(h²)/O(h),λ_h=π²+O(h²)의 전개 및 위 fixture의 +독립 계산에 근거한다. 집중 source의 특이 해에 이 smooth order를 요구하지 않는다. + +## 7. 고정5점 slab: 허용 이산 오차와 plateau + +Unit area,h=1,k=rho=c=1,Ti=0,Ts=1,bottom step/top 단열이며 면내 node 조건은 같다. +x=(0,1/4,1/2,3/4,1) T,six trace x=(0,1/4,1/2−,1/2+,3/4,1)의 q3, +E=1ᵀCT,bottom 전체 Rc를 비교한다. Unit square1요소 및2×2 mesh다. +t=0은 rawTi이며 flux/reaction을 생성하지 않는다. + +무차원 시간 t*=αt/h²,T/ΔT,q/(kΔT/h),E/(ρc area hΔT),R/(k areaΔT/h). +이 fixture의 척도는 모두1이다. T/flux는 위 위치의 max absolute error, +E/R은 absolute error다. Station 상대오차나 near-zero floor를 만들지 않는다. + +### 7.1 Fourier 및 독립 P2 spectral oracle + +정식화 §10 series에서 M=80,최소 t*=.001. +β=π²t*/4,r=exp[-8β(M+1)]이면 온도 tail≤ +4exp[-β(2M+1)²]/[π(2M+1)(1-r)],flux/reaction tail≤ +2exp[-β(2M+1)²]/(1-r),energy tail≤temperature tail. +최소 시간에서도 모든 무차원 tail은 **1e-27 미만**이다. +1−r는 -expm1(-8β(M+1))로 계산한다. Summation roundoff는 +γ_(8M+16)×absolute term sum으로 전파하며 transcendental 평가의 interval 또는 +확인한 library error bound도 더한다. + +독립1D P2 FE는 길이1/2인 두 segment를 조립한다. +C_seg=(1/2)/30 [[4,2,-1],[2,16,2],[-1,2,4]], +K_seg=2/3 [[7,-8,1],[-8,16,-8],[1,-8,7]]. +Bottom elimination 뒤4×4 generalized eigenproblem의 Cholesky+Jacobi 계산은 +3D product kernel을 호출하지 않는다. + +Raw 초기 정책에서 δ→0 자유 초기 극한은 +T_f(0+)=-C_ff^-1C_fc=(-5/34,3/17,-1/34,1/17). +y0=T_f(0+)−1,C-orthonormal v_j의 a_j=v_jᵀC_ff y0로 +T_f(t)=1+Σa_jv_j exp(-λ_jt). +이는 첫 finite-step RHS의 극한이며 FESA t=0 state/dataset을 변경하지 않는다. + +λ=(2.4686647564,22.9461660098,77.0631371712,198.6985026508). +Scratch eigen residual3.994e-16,C-orthogonality8.069e-16. +Solver 허용치는 이 숫자 자체가 아니라 §3의 coefficient/eigen/residual enclosure다. +BE oracle는 exp(-λt)를 (1+λt/N)^-N으로 대체한다. +N>1의 마지막 반력은 modal derivative -λ×final amplitude로 구할 수 있어 +작은 δ의 subtraction cancellation을 피한다. + +### 7.2 사전 고정 물리 오차와 허용치 + +각 t*마다 별도 N=4096,8192,δ=t*/N,INC=N run. +아래 첫 값은 **δ→0 P2 공간오차**, 괄호는 위로 반올림한 **D_space**다. +초기의 큰 허용치는 method의 제한이며 제품 전 범위 정확도 보증이 아니다. + +| t*/sqrt(t*) | T∞ | q∞ | E | Rc | +| --- | --- | --- | --- | --- | +| .001/.03162 | .1399996(.141) | 10.64818(10.7) | .03031345(.0304) | 10.86927(10.9) | +| .01/.10000 | .01524952(.0153) | 1.274464(1.28) | .00467600(.00468) | .823687(.824) | +| .05/.22361 | .00597851(.00598) | .252470(.253) | .00129690(.00130) | .00631488(.00632) | +| .1/.31623 | .001387325(.00139) | .135427(.136) | .000665986(.000666) | .009524106(.00953) | +| .5/.70711 | .000141091(.000142) | .0283969(.0284) | .0000944866(.0000945) | .000525742(.000526) | +| 1/1 | .0000413942(.0000414) | .00815759(.00816) | .0000158234(.0000159) | .0000478014(.0000479) | + +Continuum 오차 허용치는 **D_space+D_time(N)+B_oracle+B_solver**다. +D_time은 독립 modal 차이의 max/absolute norm이다. 온도/flux sampling L, +energy row1ᵀC,reaction row K_cf−C_cf C_ff^-1K_ff 각각에 대해 +LΣa_jv_j[(1+λ_jt/N)^-N−exp(-λ_jt)]를 계산한다. +Cancellation이 불확실하면 Σ|La_jv_j|×|amplitude차이|로 enclosure한다. +B_oracle은 series tail+summation/eigen interval,B_solver는 §3 행렬/solve 오차를 +L에 전파한 값이다. **B_oracle+B_solver>D_space/1000이면 plateau 판정은 +inconclusive**다. 이 조건은 geometry 지원 cutoff가 아니다. + +N=4096의 독립 D_time은 아래와 같다. N=8192에서는 대략 절반이며 appendix가 재현한다. + +| t* | T∞ | q∞ | E | Rc | +| ---: | ---: | ---: | ---: | ---: | +| .001 | 7.5820e-7 | 8.7770e-6 | 4.6839e-8 | 6.0181e-6 | +| .01 | 1.6648e-5 | 1.9198e-4 | 1.8013e-6 | 1.5383e-4 | +| .05 | 2.9817e-5 | 2.2173e-4 | 7.3403e-6 | 1.9484e-4 | +| .1 | 3.0985e-5 | 1.9790e-4 | 1.1177e-5 | 1.6441e-4 | +| .5 | 6.8903e-5 | 1.1406e-4 | 4.3931e-5 | 1.0879e-4 | +| 1 | 8.0305e-5 | 1.3224e-4 | 5.1125e-5 | 1.2621e-4 | + +물리 오차 허용치 외에 **동시에** FESA와 finite-N BE oracle가 +B_oracle+B_solver 내에서 일치해야 한다. 큰 초기 공간오차가 잘못 구현된 방법까지 +허용하지 않도록 분리했다. Table rounding은 production 결과와 무관하게 고정했다. + +### 7.3 Plateau 정량 보고 + +E_space=||oracle_P2−Fourier||와 FESA error e_N에 대해 +|e_N−E_space|≤D_time+B_oracle+B_solver를 확인한다. +면내1×1과2×2의 slab 차이는 §3 roundoff 이내여야 한다. +두께 지배 해가 면내 refinement로 좋아진다고 판정하면 실패다. + +관측 plateau의 기준은 D_time+B_oracle+B_solver≤E_space/10. +N=8192에서 .001/.01/.05/.1의 네 quantity는 충족한다. 뒤 시각에서 미충족하는 +T/E/R은 **미도달**로 보고하고 δ→0 spectral 공간오차를 따로 제시한다. +예컨대 t*=1의 E는 시간오차가 공간오차보다 커서 두 N만으로 plateau가 관측되었다고 +할 수 없다. 침투 깊이 sqrt(αt),station 간격 h/4,segment 길이 h/2,E_space,D_time을 +함께 보고한다. 초기 큰 flux/reaction 오차는 음수 온도 제거로 고칠 수 없다. + +## 8. 시간격자 predicate, count와 재분해 + +### 8.1 독립 반증과 수정 범위 + +Arithmetic-only D≤d rad(q)+rad(p)는 입력 표현오차를 모두 포함하지 않았다. +τ=5.351455173733616,d=.17262758624947144,m=31에서 +exact |τ−md|=8.326672684688674e-16 > 기존 Bhat=7.507374033534679e-16인데, +입력 rounding-cell I(τ)와 mI(d)는 겹친다. 이 반례를 blocking으로 보고했고 +승인된 입력 표현오차 의미에 맞춰 **input-cell overlap**으로 보완했다. +이전 predicate의 실행 성공을 의미 정합 승인으로 대체하지 않았다. +Appendix D는 이 실패 이력을 보존하며 현재 판정은 C를 따른다. + +### 8.2 현재 고정 predicate + +입력 double x의 원래 rounding-cell은 +L_x=(prev(x)+x)/2,U_x=(x+next(x))/2다. Significand LSB가 짝수면 두 endpoint +포함, 홀수면 두 endpoint 제외(nearest-even). 최대 finite의 U_x는 next=Inf를 +평균하지 않고 x+(x−prev(x))/2라는 유한 real dyadic threshold다. +최소 subnormal의 half-ULP도 binary64로 먼저 반올림하지 않는다. + +τ,d를 공통 power-of-two로 τ∈[1/2,1)에 옮겨 q=fl(τ/d), +m=nearest-even integer(q)를 만든다. **원래 입력의 cell도 동일 scale로 옮긴다.** +Normal화된 값의 이웃으로 cell을 다시 만들면 subnormal 입력의 표현오차를 축소하므로 +잘못이다. Endpoint-close는 I(τ)와 mI(d)의 exact 교집합이 비어 있지 않은 경우다. +두 구간이 한 점에서만 닿으면 그 점이 양쪽에 모두 포함되는지도 확인한다. + +INC≤2^32−1,q≤INC+1 및 위 scale에서 cell endpoint는 고정 dyadic 격자 2^-88에 +정확히 놓인다. Endpoint 정수와 m곱의 보수적 최대122bit는 두 uint64 또는 네 uint32 +limb 안에 들어간다. Carry/shift를 검사하는 bounded 정수 산술이며 새 arbitrary-precision +runtime dependency가 아니다. Reviewer의 Fraction은 이 bounded 비교의 독립 oracle다. +나눗셈·곱셈 roundoff는 후보/연산 bound 설명에 쓰고 cell을 더 넓히는 별도 epsilon으로 +중복 가산하지 않는다. + +q overflow 또는 INC+1 초과는 preflight 실패; close면 N=m이다. Nonclose에서는 +X=τs2^88,Y=ds2^88,f=floor(q)로 fY≤X<(f+1)Y를 exact 비교하고 필요하면 f를 +한 번 ±1 교정한다. X=fY이면 N=f, 아니면 N=f+1이다. Division 오차를 count에 +그대로 전달하지 않는다. d≥τ는 한 증분 δ=τ의 선행 분기다. +INC 기본100,양의 uint32,N≤INC. 정상 stored time=fl(jd),마지막=입력τ, +strict increase/finite를 검사한다. Close의 마지막 δ는 **d 그대로**다. +실제 shortened만 fma(-(N−1),d,τ)를 쓰며 A/RHS/reaction/energy/dissipation은 같은 δ다. + +| fixture | 필수 결과 | +| --- | --- | +| τ=10,d=3 | N=4,stored3,6,9,10;δ=3,3,3,1;마지막 refactor | +| τ=2,d=3 및 d=τ | 한 증분 δ=τ | +| τ=.9,d=.3 | N=3,stored 마지막 .9;세 δ=.3,factorization 재사용 | +| τ=1,d=.1 | N=10;stored .30000000000000004 등을 보존,δ=.1 유지 | +| nextafter(.9,-∞),d=.3 | .8999999999999999:close,N=3,δ=.3 | +| 두 ULP 아래 τ,d=.3 | .8999999999999998:nonclose,N=3,finalδ=.2999999999999998 | +| nextafter(.9,+∞),d=.3 | .9000000000000001:nonclose,N=4,finalδ=1.6653345369377348e-16 | +| 기존 반례 τ=5.351455173733616,d=.17262758624947144 | input cells overlap,N=31,모든 δ=d;가짜 tiny 증분 제거 | +| Subnormal midpoint | a=최소 positive subnormal:τ=20a,d=6a는 even/even 닫힌 접촉,N=3;τ=17a,d=5a는 odd/odd 열린 접촉,N=4,마지막δ=2a | +| Scale/최대값 | (.9,.3)×2^-1000 및×2^1000은 같은 count/close;τ=1e-323,d=5e-324는 N=2;τ=maxfinite,d=fl(maxfinite/3)는 N=3 | +| INC | N/N−1 경계,생략100;0,-1,1.5,4294967296 실패;τ=2^32,d=1,INC=2^32−1은 사전 count 실패 | +| 비유한/비표현 | Inf/NaN/0/음수,binary64에서0으로 underflow하는 token 거부;거대한 τ/d는 cast/allocation 전 실패 | +| 저장 가능성 | finite-positive/uint32 count 제한 아래 정상 간격 d는 ulp(jd)보다 충분히 크다. Generator의 finite/strict-increase 보호 분기는 중복/Inf 시간을 주입한 내부 fixture로 검사;불가능한 정상 공개 입력을 만들지 않음. | + +특히 nextafter(.9,+∞)/.3의 tiny remainder를 유지하는 근거는 단순 predicate 결과가 +아니다. exact |τ−3d|=1.6653345369377348e-16이고 두 입력 cell 반폭 합은 +1.3877787807814457e-16이므로 **실제 rounding cells가 분리된다**. 이를 가짜 표현오차 +증분으로 흡수하면 입력 표현오차 범위를 넘는다. 반면 위 m=31 반례는 cells가 겹치므로 +추가 tiny 증분을 만들면 안 된다. 작은 실제 δ의 나쁜 conditioning은 §3과 candidate +보존으로 처리하며 nominal δ로 바꾸지 않는다. + +현재 predicate의 independent seed901 10,000쌍은 원래 cell Fraction overlap과 +2^-88 정수 비교가 모두 일치했고 최대122bit를 넘지 않았다. 원래 subnormal,maxfinite, +열린/닫힌 접촉,반례 및 scale 경계를 따로 검산했다. FMA oracle는 exact Fraction +product-sum을 한 번 float로 반올림한다. Production dependency 요구가 아니다. + +## 9. 외부 reference와 evidence 경계 + +필수 원본은 reference/shell-thermal/shell-thermal.inp와 +reference/shell-thermal/shell-thermal NodalTemperature.csv다. +이 검토는 원본 복사·수정·생성 또는 Abaqus 실행을 하지 않았다. +2026-09-28 read-only snapshot(729row,min=-15.5742216,max=1000)은 과거 관찰이며 +이번 solver comparison을 대신하지 않는다. + +Family=(shell-thermal,source Step-1,exact t=37,NT11 bottom,scalar, +사용자 일관 온도 단위,blocking),identity=(Part-1-1,source node label,DOF11). +CSV instance 생략은 입력 단일 instance와 label 집합 검증으로만 복원한다. +Missing/extra/duplicate/nonfinite/frame/instance/station 오류는 tolerance 전에 실패한다. + +공통 TOLERANCE 그대로 S=max|reference|=1000,B=.01S=10; +|r|≤B이면 |f−r|≤B,그 밖에는 |f−r|/|r|≤.05, +모든729row RMS/S≤.01을 동시에 만족한다. S=0이면 모든 값 exact0이다. +−15.5742216은 relative 대상이며 clamp/온도 영점 이동 금지. +추가 Abaqus CSV는 필수가 아니다. 한 frame으로 내부 station/flux/reaction/시간 이력을 +검증했다고 할 수 없다. + +## 10. 실행 기록과 완료 경계 + +git rev-parse --show-toplevel로 T1 worktree를 확인했다(exit0). +NumPy/SciPy import 환경 probe는 NumPy 부재로 exit1;설치 없이 표준 라이브러리 +scratch로 전환했다. [Appendix](numerical-oracle-scratch.md)의 명령은 Python blocks를 +지정 순서로 메모리에서 실행한다. 실제 scratch exit0의 핵심 출력: + +- Fraction:C positive pivot20,K positive leading pivot19/K1=0,1ᵀC1=2; + jump−1/5,1/25 및 reaction24/5,36/25. +- 직접54점/분리식:다섯 geometry 상대 norm 차이 최대7.73e-16. +- Slab eigen residual3.994e-16/C-orthogonality8.069e-16,§7 공간/시간 오차. +- §5 quadrature,§6 예상 차수,§8 time predicate/INC assertions. + +최종 재현 명령은 appendix 첫 PowerShell block이며 exit 0, +마지막 출력은 ALL SCRATCH BLOCKS COMPLETED였다. Input-cell 독립 비교는 +10,000/10,000 및 CURRENT TIME GRID ASSERTIONS PASSED를 출력했다. +상대 링크 파일 존재/Markdown 표 열 수 정적 검사도 exit 0(8 links,2 files), +두 소유 파일만 대상으로 한 git diff --check는 exit 0이었다. + +수학 scratch는 production correctness를 증명하지 않는다. **미실행:**MSVC build/CTest, +구현 단위·통합 검사,FESA 결과 생성,NT11 comparison,실제 solver convergence/energy/ +physics 검토. 필요한 수치 기준은 고정되었으며 구현 결과 뒤에 선택할 미정 tolerance는 +없다. 필수 실패 시 구현·conditioning·oracle 전제를 조사하고 physics/method/acceptance +충돌은 메인에 보고한다.