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

docs(thermal): account for input rounding cells in time grid

Browse Source
This commit is contained in:
KOKO\Mimi
2026-09-29 18:03:02 +09:00
parent efa5403b4a
commit dd409b57a1
1 changed files with 150 additions and 31 deletions
Download Patch File Download Diff File Expand all files Collapse all files
docs/shell-transient-heat-transfer/formulation.md
+150 -31
View File
@@ -310,41 +310,57 @@ raw 초기 상태의 projection으로 바꾸지 않는다.
기본값 100이며 허용 증분 수의 상한이다. `d≥τ`이면 한 증분, `δ_0=τ`, 저장 시각 `τ`다.
그 외에는 **반복 덧셈 없이** 정수 `j`의 곱으로 정상 시각을 만든다.
정수비 근방은 다음 연산 bound로 한 번만 판정한다. 먼저 비율이 `INC+1` 이내인지
확인하고 `τ`를 `[1/2,1)`에 놓는 같은 exact power-of-two로 `τ,d`를 재척도한다.
이 범위에서는 재척도된 `d`도 normal이며 아래 `p`의 overflow를 피할 수 있다.
아래 predicate는 재척도된 변수로 평가하며 결과와 시각·δ는 원래 단위를 유지한다.
`fl`은 binary64 nearest-even이고,
`rad(x)=max(x-nextafter(x,-∞), nextafter(x,+∞)-x)/2`는 해당 값의 보수적인
rounding-cell 반폭이다. 반폭이 subnormal로 표현되지 않으면 powers-of-two scale에서
계산한다. 여기서는 `q=fl(τ/d)`, `m=nearest_integer(q)`, `p=fl(m d)`를 쓴다.
`q`가 overflow하거나 `q>INC+1`이면 증분 수 초과다. 나머지 경우 `m≤2³²` 부근의
정수는 double에 정확히 표현된다. Tie의 정수 선택도 nearest-even으로 고정한다.
정수비 근방의 동등성은 **입력 binary64로 반올림될 수 있는 실제 시간 구간**으로 판정한다.
입력 중심값을 오차 없는 물리량으로 간주하지 않는다. `fl`은 nearest-even이며 먼저
`q=fl(τ/d)`를 구한다. `INC+1`은 uint64 범위에서 계산한다.
`q`가 overflow하거나 `q>INC+1`이면 증분 수 초과다.
나머지는 `m=nearest_integer(q)`를 후보로 쓴다. Tie는 nearest-even으로 고정하며
`m≤2³²`는 정확히 표현된다. `INC`에 맞추려고 다른 정수 후보로 바꾸지 않는다.
원래 입력 `x>0`의 predecessor/successor를 `x−,x+`라 하면 rounding cell은
\[
B=d\,\operatorname{rad}(q)+\operatorname{rad}(p),\qquad
\widehat B=\operatorname{nextafter}(\operatorname{fl}(B),+\infty),\qquad
D=|\operatorname{fma}(-m,d,\tau)|,\qquad
\text{endpoint-close}\iff m\ge1\ \land\ D\le\widehat B.
I(x)=\left\langle (x^-+x)/2,\ (x+x^+)/2\right\rangle,
\qquad \text{endpoint-close}\iff m\ge1\ \land\ I(\tau)\cap mI(d)\ne\varnothing.
\]
이는 임의 시간 epsilon이 아니다. 나눗셈 반올림의 dimensional uncertainty는
`d rad(q)`, 곱셈의 uncertainty는 `rad(p)`이므로 합이 두 연산에 대한 보수적 bound다.
특히 quotient가 정수 `m`으로 반올림되는 경우
`|τ-md|≤d rad(q)`이고, 정상 시각 곱 자체의 표시 불확실성은 `rad(p)`다.
Fused remainder로 곱의 선행 반올림을 제거하므로 `τ-fl(md)`의 이중 오차를 비교하지 않는다.
이 bound 안의 endpoint는 정상 마지막 증분과 같은 계수를 사용하는 명시적 정책이다.
실제 remainder가 이 범위 밖이면 작더라도 보존한다. 외부 입력 숫자의 임의 오차나 물리
시간 tolerance를 추가하지 않는다. 재척도된 `rad(q)`는 power-of-two이고 `d rad(q)`는
normal 범위의 exact binary scaling이다. 따라서 `fl(B)`는 덧셈 한 번의 반올림뿐이며
`nextafter` 한 번이 보수적 상계를 준다. FMA 한 번의 잔여 반올림까지 포함하면
accepted endpoint의 정확한 거리 bound는 `|τ-md|≤Bhat+rad(D)`다. 이렇게 **표현 가능한
threshold에 대한 `≤` 한 번**으로 판정하여 경계에서도 결과를 고정한다. FMA를 두 개의
일반 연산으로 분해하거나 compiler 재결합에 맡기지 않는다.
꺾쇠의 두 끝점은 `x`의 binary64 significand 최하위 bit가 0(even)이면 포함하고,
1(odd)이면 제외한다. Subnormal도 같은 규칙이며 최소 양수의 predecessor는 0이다.
최대 finite의 successor는 부동소수점 `+∞` 대신 **수학적인 다음 lattice 값 `2^1024`**로
둔다. 최대 finite significand는 odd이므로 overflow midpoint를 cell에 포함하지 않는다.
Midpoint를 double로 반올림하지 않고 아래 정수·지수 표현으로 정확히 구성한다.
입력 변환은 correctly-rounded nearest-even을 전제로 하며 임의 입력 오차는 더하지 않는다.
1. Endpoint-close이면 `N=m`. 아니면 `N=floor(q)+1=ceil(τ/d)`다. Quotient 반올림으로
정수 경계를 넘을 수 있는 경우는 위 endpoint-close bound 안에 이미 포함되므로
이 분기에서 floor는 안전하다. 먼저 nearest-integer branch를 판정하는 순서를 유지한다.
`I(τ)=(a_L,a_U)`, `mI(d)=(b_L,b_U)`의 끝점 포함 flag를 보존하여
`L=max(a_L,b_L)`, `U=min(a_U,b_U)`를 구한다. `L<U`면 close, `L>U`면 nonclose다.
`L=U`면 그 접점이 **두 cell 모두에 포함될 때만** close다. 나눗셈·곱셈 오차와 입력 cell
폭을 중복 가산하는 epsilon은 없다. `q`는 정수 후보만 만들고, 겹침 판정에는 반올림된
`fl(md)`를 사용하지 않는다. 같은 endpoint를 허용하지 않는 cell은 잔여시간이 작아도 구별한다.
정확한 predicate는 임의정밀도 dependency 없이 다음 고정 폭 연산으로 구현할 수 있다.
- `e=frexp(τ).exponent`, `s=2^(−e)`로 두어 `sτ∈[1/2,1)`로 만든다. **원래 입력의
cell을 먼저 구성한 뒤** `s`를 곱한다. Normal화된 값의 prev/next로 cell을 다시 만들면
원래 subnormal 입력의 표현 범위를 줄이므로 금지한다. `s` 자체를 double로 만들 필요는 없다.
- `d<τ`, 유한 `q≤INC+1≤2³²`이면 정확한 `τ/d`는 `2³²+2^(−21)` 이하이므로
`sd>2^(−34)`이다. Scaled 중심값과 원래 cell의 midpoint는 모두 lattice `2^(−88)`의
정수배다. Normal 입력의 최소 반간격은 이 lattice보다 작지 않으며, subnormal은 더
성긴 원래 간격을 유지한다. Binade 경계의 비대칭 간격도 따로 보존한다.
- `X=τs2^88`, `Y=ds2^88`와 네 cell endpoint를 원래 binary64 significand/exponent에서
shift·덧셈으로 정확히 구성한다. Scaled cell 끝점은 2보다 작으므로 각 정수는 89 bit 이하다.
`m` 및 floor 교정용 `floor(q)+1`은 최대 `2³²+1`이므로 곱은 **122 bit 이하**다.
- 두 `uint64` 상당의 high/low word로 저장·비교한다. High word를 먼저 비교하고, 같으면
low word를 비교한다. 덧셈/뺄셈의 carry/borrow와 shift의 word 경계를 명시적으로 처리한다.
곱은 네 base-`2³²` digit와 정수 승수의 최대 두 digit로 계산할 수 있다. 각
`a_i*b_j+result_digit+carry≤2^64−1`이므로 중간값도 uint64에 들어가며 다음 digit으로
carry를 넘긴다. 최상위 128 bit 밖 carry가 없음을 검사한다. `__int128`, extended float,
compiler 연산 재결합 또는 arbitrary-precision library에 의존하지 않는다.
1. Endpoint-close이면 `N=m`. 아니면 중심값의 정확한 `N=ceil(τ/d)`를 쓴다.
`f=floor(q)`에서 정수 비교 `fY≤X<(f+1)Y`를 확인하고 필요하면 `f`를 한 번 ±1
교정한다. `q≤2³²`의 division 오차는 1보다 작으므로 한 번이면 충분하다.
`X=fY`이면 `N=f`, 아니면 `N=f+1`이다. Division 반올림을 증분 수에 그대로 전달하지 않는다.
`N`은 `INC`와 비교하기 전 uint32로 축소하지 않는다.
2. `N>INC`이면 solve 전에 실패한다. 정상 저장 시각은 `t_j=fl(j d)` (`j=1,…,N−1`),
마지막은 **입력 그대로** `t_N=τ`다. 모든 시각이 유한하고 엄격히 증가하며
`t_{N−1}<τ`인지 preflight에서 확인한다. `t_0=0`이다.
@@ -352,7 +368,9 @@ threshold에 대한 `≤` 한 번**으로 판정하여 경계에서도 결과를
마지막 증분만 `δ=fl(τ-(N−1)d)`로 만든다. 이 remainder는
`fma(-(N−1),d,τ)`처럼 곱-차를 한 번에 반올림한다. `δ>0`·유한성을 검사한다.
저장된 `fl(jd)`끼리 다시 빼서 정상 `δ`를 만들지 않는다. 개념적인 격자 `j d`와
저장 표현의 차이 및 endpoint-close 조정은 위 rounding bound로 제한된다.
저장 표현의 차이는 각 곱의 rounding cell로 제한된다. Endpoint-close 조정은 입력 cell
overlap으로 제한된다. 원래 단위의 `|τ−md|`는 해당 방향의 `τ` 반간격과 `m`배 `d`
반간격 합 이내이고, 계산된 FMA remainder에는 한 번의 최종 rounding만 남는다.
4. 해당 증분의 `A`, RHS, full residual, reaction, `ΔE` 수지, BE 소산은 모두 **같은
실제 `δ`**를 사용한다. Endpoint-close는 동일 계수를 재사용하며 실제 shortened
마지막 `δ`가 `d`와 다를 때 refactorize한다.
@@ -363,6 +381,107 @@ threshold에 대한 `≤` 한 번**으로 판정하여 경계에서도 결과를
`nextafter`로 endpoint-close 경계를 양쪽에서 이동한 사례, 아주 작고 큰 시간 척도,
`INC=N`/`N−1`, uint32 overflow 방지를 독립 검증한다. 이 격자는 adaptive time stepping이 아니다.
### 7.2 입력 cell 수정의 실제 검산 기록
2026-09-29 Python 3.11 `Fraction`을 독립 exact-dyadic oracle로 사용한 아래 scratch를
실행했다(exit 0). 이는 문서 수식 검산이며 solver build/test나 수치 correctness 승인이 아니다.
기존 `Bhat=nextafter(fl(d rad(q)+rad(fl(md))),+∞)`는 입력 cell을 누락하여 첫 사례를
nonclose로 오판했다. Scratch의 임의정밀도는 검산용이며 production dependency가 아니다.
```text
old-failure: m=31, exactD=8.326672684688674e-16 > Bhat=7.507374033534679e-16
corrected: close=True, N=31; INC=31 accept, INC=30 reject
tau=nextafter(.9,+inf), d=.3: close=False, N=4, delta=1.6653345369377348e-16
tau=.9, d=.3: close=True, N=3
tau=20*minsub, d=6*minsub: even/even midpoint included, close=True, N=3
tau=17*minsub, d=5*minsub: odd/odd midpoint excluded, close=False, N=4
power-of-two scales [-900,-400,0,400,900]: 5/5
maxfinite, original-subnormal, binade cases: 5/5
seed41 input-cell versus fixed-lattice predicate: 10000/10000
uint32 INC boundary: N=4294967295 accept; N=4294967296 reject
```
<details>
<summary>실행한 재현 scratch (Python 표준 라이브러리)</summary>
```python
import math, random, struct
from fractions import Fraction as F
def cell(x):
a=F(x); p=F(math.nextafter(x,0.0))
nx=math.nextafter(x,math.inf)
n=F(nx) if math.isfinite(nx) else a+(a-p)
even=(struct.unpack('>Q',struct.pack('>d',x))[0]&1)==0
return (a+p)/2,(a+n)/2,even
def overlap(a,b):
lo=max(a[0],b[0]); hi=min(a[1],b[1])
return lo<hi or (lo==hi and a[2] and b[2])
def predicate(t,d):
m=round(t/d); a=cell(t); b=cell(d)
return overlap(a,(m*b[0],m*b[1],b[2]))
def fixed(t,d):
e=math.frexp(t)[1]; s=F(2)**(-e+88)
a=cell(t); b=cell(d); m=round(t/d)
vals=[a[0]*s,a[1]*s,b[0]*s,b[1]*s,F(t)*s,F(d)*s]
assert all(v.denominator==1 for v in vals)
al,ah,bl,bh,x,y=map(int,vals)
assert max(al,ah,m*bl,m*bh,x,y).bit_length()<=122
return overlap((al,ah,a[2]),(m*bl,m*bh,b[2]))
def count(t,d):
if d>=t: return 1
if predicate(t,d): return round(t/d)
z=F(t)/F(d); return -(-z.numerator//z.denominator)
def count_fixed(t,d):
if d>=t: return 1
if fixed(t,d): return round(t/d)
e=math.frexp(t)[1]; s=F(2)**(-e+88)
x,y=int(F(t)*s),int(F(d)*s); f=math.floor(t/d)
if f*y>x: f-=1
elif (f+1)*y<=x: f+=1
assert f*y<=x<(f+1)*y
return f if x==f*y else f+1
def rad(x):
return max(F(x)-F(math.nextafter(x,-math.inf)),
F(math.nextafter(x,math.inf))-F(x))/2
t,d=5.351455173733616,.17262758624947144
m=round(t/d); q=t/d; p=m*d
old=math.nextafter(float(F(d)*rad(q)+rad(p)),math.inf)
distance=abs(F(t)-m*F(d))
assert distance>F(old) and predicate(t,d) and fixed(t,d)
print('old-failure:',m,float(distance),old,distance<=F(old))
print('corrected:',predicate(t,d),'N=',count(t,d),'INC31=',count(t,d)<=31,'INC30=',count(t,d)<=30)
t2=math.nextafter(.9,math.inf)
print('real-remainder:',predicate(t2,.3),'N=',count(t2,.3),'delta=',float(F(t2)-3*F(.3)))
print('decimal:',predicate(.9,.3),count(.9,.3))
tiny=math.ulp(0.0)
for a,b in [(20,6),(17,5)]:
print('midpoint:',a,b,cell(a*tiny)[2],cell(b*tiny)[2],predicate(a*tiny,b*tiny),count(a*tiny,b*tiny))
for e in [-900,-400,0,400,900]:
tt,dd=math.ldexp(t,e),math.ldexp(d,e)
assert predicate(tt,dd) and fixed(tt,dd) and count(tt,dd)==31
print('power-of-two-scales: 5/5')
maxf=float.fromhex('0x1.fffffffffffffp+1023')
for tt,dd in [(maxf,maxf/31),(20*tiny,6*tiny),(17*tiny,5*tiny),
(1.,math.nextafter(.25,0.)),(1.,math.nextafter(.25,math.inf))]:
assert predicate(tt,dd)==fixed(tt,dd)
assert count(tt,dd)==count_fixed(tt,dd)
print('extreme/binade cases: 5/5')
rng=random.Random(41)
for _ in range(10000):
dd=math.ldexp(rng.uniform(.5,1),rng.randint(-900,900))
mm=rng.randint(2,2**32)
tt=mm*dd
for _ in range(rng.randint(0,2)):
tt=math.nextafter(tt,rng.choice([0.,math.inf]))
assert predicate(tt,dd)==fixed(tt,dd)
print('seed41 cell-vs-fixed: 10000/10000')
assert count(float(2**32-1),1.)==2**32-1
assert count(float(2**32),1.)==2**32
print('uint32 INC limit: N=4294967295 accept; N=4294967296 reject')
```
</details>
## 8. 반력 열률, 에너지와 안정성
반력은 하중 제거 전 full residual에서 구한다.