diff --git a/docs/shell-transient-heat-transfer/numerical-oracle-scratch.md b/docs/shell-transient-heat-transfer/numerical-oracle-scratch.md index e13b4a9..7397a3d 100644 --- a/docs/shell-transient-heat-transfer/numerical-oracle-scratch.md +++ b/docs/shell-transient-heat-transfer/numerical-oracle-scratch.md @@ -13,7 +13,7 @@ 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 +assert len(blocks) == 5 scope = {} for block in blocks: exec(compile(block, str(p), 'exec'), scope) @@ -172,23 +172,6 @@ for name,X in cases.items(): 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): @@ -219,33 +202,43 @@ orth=fn([[v-float(i==j)for j,v in enumerate(row)]for i,row in enumerate(mm(mm(tr 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))) + +def spatial_bounded(n,steps=32,t=.001,ng=12): + h=1/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 + et=eg=etime=ese=0. for i in range(n): for s,w in gauss(ng): - x=(i+(s+1)/2)*h + x=(i+(s+1)/2)*h;measure=w*h/2 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=[] + g=(math.cos(math.pi*(i+1)/n)-math.cos(math.pi*i/n))/h + et+=(af*v-ae*math.cos(math.pi*x))**2*measure + eg+=(af*g+math.pi*ae*math.sin(math.pi*x))**2*measure + etime+=((af-ase)*v)**2*measure + ese+=(ase*v-ae*math.cos(math.pi*x))**2*measure + return math.sqrt(et),math.sqrt(eg),math.sqrt(etime),math.sqrt(ese) +rows=[] 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) + a=spatial_bounded(n);b=spatial_bounded(n,ng=16) + ratio=a[2]/a[3];quad=max(abs(x-y)/abs(x)for x,y in zip(a,b)) + assert ratio<1/32 and quad<1e-8 + rows.append(a) + fluxbytes=n*n*24*3*32*8 + tempbytes=(n+1)**2*5*33*8 + print("bounded space",n,"T,grad,time,semi",a,"time/space",ratio,"quad",quad,"historybytes",fluxbytes+tempbytes) +for a,b in zip(rows,rows[1:]): + pt=math.log(a[0]/b[0],2);pg=math.log(a[1]/b[1],2) assert 1.95<=pt<=2.05 and .98<=pg<=1.03 - print("spatial BE order",pt,pg) + print("bounded order",pt,pg) +time_orders=[] +for n in [16,32,64,128]: + e1=abs(math.exp(-n*math.log1p(1/n))-math.exp(-1)) + e2=abs(math.exp(-2*n*math.log1p(1/(2*n)))-math.exp(-1)) + order=math.log(e1/e2,2) + assert .97<=order<=1.02 + time_orders.append(order) +print("fixed spatial time orders",time_orders) ``` @@ -423,3 +416,89 @@ 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))) ``` + +## E. Library ULP 가정 없는 Fourier interval oracle + +고정 dyadic 320-bit outward enclosure, Machin π, Taylor sin/cos/exp와 analytic +remainder를 사용한다. Python Fraction은 검증 generator 전용이다. 출력 interval의 +outward binary64 hex endpoints를 C++ fixture에 고정할 수 있다. FESA solver와 +C++ std::sin/std::exp의 ULP 보장을 가정하지 않는다. + +```python + +from fractions import Fraction as F +import math +BITS=320;SCALE=1<0 or b[1]<0 + return mul(a,(F(1)/b[1],F(1)/b[0])) +def scale(a,c):return mul(a,interval(c)) +def atan_bounds(q): + a=sum((F(-1)**j)*q**(2*j+1)/(2*j+1) for j in range(64)) + tail=q**129/129 + return a,a+tail +a,b=atan_bounds(F(1,5)),atan_bounds(F(1,239)) +PI=interval(16*a[0]-4*b[1],16*a[1]-4*b[0]) +assert PI[1]-PI[0]=0 + if x[0]>=512:return interval(0,F(2)**-512) + assert x[1]<=512 + y=scale(x,F(1,1024));term=interval(1);s=term + for j in range(1,65): + term=scale(mul(term,neg(y)),F(1,j));s=add(s,term) + tail=F(1,2)**65/math.factorial(65) + s=add(s,interval(-tail,tail)) + for _ in range(10):s=mul(s,s) + return s +def sincos_pi8(k): + k%=16;quad=(k+2)//4;rem=k-4*quad + x=scale(PI,F(rem,8));x2=mul(x,x) + st=x;ss=x;ct=interval(1);cs=ct + for j in range(1,32): + st=scale(mul(st,neg(x2)),F(1,2*j*(2*j+1)));ss=add(ss,st) + ct=scale(mul(ct,neg(x2)),F(1,(2*j-1)*2*j));cs=add(cs,ct) + ss=add(ss,interval(-F(1,math.factorial(65)),F(1,math.factorial(65)))) + cs=add(cs,interval(-F(1,math.factorial(64)),F(1,math.factorial(64)))) + q=quad%4 + return [(ss,cs),(cs,neg(ss)),(neg(ss),neg(cs)),(neg(cs),ss)][q] +def series_interval(t): + Ts=[interval(1)for _ in range(5)];qs=[interval(0)for _ in range(5)];E=interval(1) + for m in range(80): + odd=2*m+1 + x=scale(mul(PI,PI),F(t)*odd*odd/4);dec=exneg(x) + for j in range(5): + si,co=sincos_pi8(odd*j) + Ts[j]=sub(Ts[j],scale(div(mul(si,dec),PI),F(4,odd))) + qs[j]=add(qs[j],scale(mul(co,dec),2)) + E=sub(E,scale(div(dec,mul(PI,PI)),F(8,odd*odd))) + assert max(hi-lo for lo,hi in Ts+qs+[E])2)을 사용한다. +5. M=80 series를 위 interval로 합한다. Tail을 더하기 전 각 quantity enclosure 폭은 + **2^-270 미만**이어야 한다. 여기에 앞의 analytic Fourier tail을 더한다. + 필요한 binary64 fixture는 lower를 float 변환 후 nextafter(-∞),upper를 변환 후 + nextafter(+∞)하여 export한다. 이 endpoint 차이가 B_oracle이며 별도 floor가 없다. + +실제 generator에서 최소 시간의 tailT=1.658661e-30,tailq=4.194724e-28, +최종 enclosure 최대 폭8.389448e-28이었다. 나머지 시각의 폭은8.93e-89 이하였다. +이 방식은 테스트 oracle generator 전용이며 production arbitrary-precision dependency가 +아니다. Spectral oracle의 exp도 동일 enclosure를 적용하고 eigenpair/coefficients의 +backward-error interval을 함께 전파한다. Appendix A의 통상 math 계산은 표 수치의 +간단한 재현용이며 interval E가 transcendental correctness의 근거다. + +T10은 generator의 T/q/E/R interval을 outward binary64 fixture로 고정하여 사용할 수 +있다. 또는 binary64 outward Taylor가 동일한 최종 B_oracle 및 §6–7의 오차 budget을 +증명하면 그 경로도 가능하다. 2^-270은 고정 고정밀 generator의 내부 검사이며 production +정밀도 요구가 아니다. 무차원 interval error를 실제 단위로 옮길 때 T는 |ΔT|, +q는 k|ΔT|/h,E는 rho c area h|ΔT|,R은 k area|ΔT|/h를 곱한다. +온도 기준 이동은 interval 폭을 바꾸지 않는다. 선형 sampling/trace/energy/reaction row L은 +성분별 |L| times input error로 전파하며 signed 값 자체를 절댓값으로 바꾸지 않는다. 독립1D P2 FE는 길이1/2인 두 segment를 조립한다. C_seg=(1/2)/30 [[4,2,-1],[2,16,2],[-1,2,4]], @@ -485,6 +577,8 @@ scratch로 전환했다. [Appendix](numerical-oracle-scratch.md)의 명령은 Py 최종 재현 명령은 appendix 첫 PowerShell block이며 exit 0, 마지막 출력은 ALL SCRATCH BLOCKS COMPLETED였다. Input-cell 독립 비교는 10,000/10,000 및 CURRENT TIME GRID ASSERTIONS PASSED를 출력했다. +실행 예산 보완 뒤 five-block 재검산도 exit 0이었다. Bounded-space order/time ratio, +FOURIER INTERVAL ASSERTIONS PASSED 및 최종 완료 출력을 확인했다. 상대 링크 파일 존재/Markdown 표 열 수 정적 검사도 exit 0(8 links,2 files), 두 소유 파일만 대상으로 한 git diff --check는 exit 0이었다.