baram2584/FESADev
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

docs(thermal): bound convergence cost and certify oracle errors

Browse Source
This commit is contained in:
KOKO\Mimi
2026-09-29 18:20:01 +09:00
parent c698bf7c97
commit 6df5a17ca8
2 changed files with 242 additions and 69 deletions
Download Patch File Download Diff File Expand all files Collapse all files
docs/shell-transient-heat-transfer/numerical-oracle-scratch.md
+117 -38
View File
@@ -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<<BITS
def down(x):return F((x*SCALE).numerator//(x*SCALE).denominator,SCALE)
def up(x):return -down(-x)
def interval(a,b=None):return (down(F(a)),up(F(a if b is None else b)))
def add(a,b):return down(a[0]+b[0]),up(a[1]+b[1])
def neg(a):return -a[1],-a[0]
def sub(a,b):return add(a,neg(b))
def mul(a,b):
v=[x*y for x in a for y in b];return down(min(v)),up(max(v))
def div(a,b):
assert b[0]>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]<F(2)**-295
def exneg(x):
assert x[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])<F(2)**-270
# Uniform analytic tail bound, using safe r upper at minimum time.
beta=scale(mul(PI,PI),F(t)/4)
first=exneg(scale(beta,161*161))
rr=exneg(scale(beta,8*81))
tailq=div(scale(first,2),sub(interval(1),rr))[1]
tailt=div(scale(first,4),scale(mul(PI,sub(interval(1),rr)),161))[1]
Ts=[add(v,interval(-tailt,tailt))for v in Ts]
qs=[add(v,interval(-tailq,tailq))for v in qs]
E=add(E,interval(-tailt,tailt))
return Ts,qs,E,tailt,tailq
for t in [.001,.01,.05,.1,.5,1.]:
Ts,qs,ee,tt,tq=series_interval(t)
assert tt<F(1,10**27) and tq<F(1,10**27)
# Certified brackets can be exported as outward binary64 hex endpoints for C++ fixtures.
width=max(float(hi-lo)for lo,hi in Ts+qs+[ee])
assert width<1e-26
print("Fourier certified",t,"max-width",width,"tailT",float(tt),"tailq",float(tq),
"E bracket",math.nextafter(float(ee[0]),-math.inf).hex(),math.nextafter(float(ee[1]),math.inf).hex())
print("FOURIER INTERVAL ASSERTIONS PASSED")
```