| """Independent checks, mandatory before data publication.""" |
| import sys,json |
| from pathlib import Path |
| sys.path.insert(0,str(Path(__file__).resolve().parents[1]/'src')) |
| from rod import Rod,Controls |
| import numpy as np |
| from scipy.integrate import solve_bvp |
| from scipy.optimize import brentq |
|
|
| def run(): |
| checks={} |
| r=Rod(24); c=Controls(chi=.004,fx=.005,fy=.002,h=.001,k1=.3,k2=-.2,k3=.6) |
| z=.4*np.sin(np.linspace(0,np.pi/2,r.d)); e,g,K=r.terms(z,c) |
| eps=1e-6; eye=np.eye(r.d) |
| gn=np.array([(r.terms(z+eps*v,c)[0]-r.terms(z-eps*v,c)[0])/(2*eps) for v in eye]) |
| Kn=np.column_stack([(r.terms(z+eps*v,c)[1]-r.terms(z-eps*v,c)[1])/(2*eps) for v in eye]) |
| checks['gradient_fd_max_abs']=float(np.max(np.abs(g-gn))) |
| checks['hessian_fd_max_abs']=float(np.max(np.abs(K-Kn))) |
| assert checks['gradient_fd_max_abs']<1e-8 and checks['hessian_fd_max_abs']<1e-8 |
| axial=r.solve(Controls(fx=.03)) |
| checks['axial_tip_error']=float(np.linalg.norm(axial['tip']-[1.03,0])) |
| assert checks['axial_tip_error']<1e-12 |
| circular=[] |
| for n in [24,48,96]: |
| rr=Rod(n); cc=Controls(k1=1.2,k2=1.2,k3=1.2) |
| ans=rr.solve(cc,np.linspace(0,1.2,n+1)[1:]) |
| exact=np.array([np.sin(1.2)/1.2,(1-np.cos(1.2))/1.2]) |
| circular.append(float(np.linalg.norm(ans['tip']-exact))) |
| checks['circular_tip_errors_n24_48_96']=circular |
| assert circular[2]<circular[1]<circular[0] and circular[-1]<1e-5 |
| rr=Rod(96); cc=Controls(chi=1e-4,h=.01,delta=.1) |
| pc=(-1+np.sqrt(1+16*cc.h*cc.delta-np.pi**2*cc.chi))/2-2*cc.h |
| def eig(p): |
| cc.fx=p |
| return rr.solve(cc)['lambda_min_over_chi'] |
| pn=brentq(eig,pc-.001,pc+.001,xtol=1e-13) |
| checks['critical_load_analytic']=pc; checks['critical_load_discrete_rod96']=pn |
| checks['critical_load_error']=abs(pc-pn) |
| assert abs(pc-pn)<1e-7 |
| |
| errors=[]; chain_errors=[] |
| for p in [-.017159,-.016428,-.016316]: |
| chi,h,d=1e-4,.01,.1; a=p+2*h |
| s=np.linspace(0,1,101) |
| def fun(s,y): |
| th=y[0]; stretch=1+a*np.cos(th) |
| return np.vstack([y[1],(a*stretch-4*h*d)*np.sin(th)/chi,stretch*np.cos(th),stretch*np.sin(th)]) |
| def bc(l,r): return [l[0],r[1],l[2],l[3]] |
| candidates=[] |
| for amp in [-.3,-1.,-2.]: |
| y=np.vstack([amp*np.sin(np.pi*s/2),amp*np.pi/2*np.cos(np.pi*s/2),s,amp*s*s/2]) |
| b=solve_bvp(fun,bc,s,y,tol=1e-9,max_nodes=10000) |
| if b.success and b.y[3,-1]<-1e-3: candidates.append(b) |
| assert candidates |
| b=min(candidates,key=lambda b:abs(b.y[0,-1])) |
| target=b.sol(np.linspace(0,1,65))[2:4].T |
| c=Controls(chi=chi,h=h,delta=d,fx=p) |
| arod=Rod(96); out=arod.solve(c,b.sol(np.linspace(0,1,97))[0,1:]) |
| assert out['accepted'] and out['tip'][1]<0 and out['stability']=='stable' |
| errors.append(float(np.sqrt(np.mean(np.sum((out['centerline']-target)**2,axis=1))))) |
| one=[] |
| for n in [15,30,60,120]: |
| rc=Rod(n,'chain'); v=rc.solve(c,b.sol(np.arange(1,n)/n)[0]) |
| assert v['accepted'] and v['tip'][1]<0 |
| one.append(float(np.linalg.norm(v['tip']-target[-1]))) |
| assert one[-1]<one[0] |
| chain_errors.append(one) |
| checks['fig9_rod96_centerline_rmse_over_L']=errors |
| checks['fig9_chain_tip_errors_n15_30_60_120']=chain_errors |
| assert max(errors)<.001 |
| |
| rc=Rod(15,'chain'); c=Controls(chi=1e-4,h=.01) |
| omega2=4*np.sin(np.pi/(2*(2*15-1)))**2 |
| expected=(-1+np.sqrt(1-4*c.chi*15**2*omega2+16*c.h*c.delta))/2-2*c.h |
| c.fx=expected; v=rc.solve(c) |
| checks['chain_critical_eigenvalue']=v['lambda_min_over_chi'] |
| assert abs(v['lambda_min_over_chi'])<1e-8 |
| checks['passed']=True |
| return checks |
|
|
| if __name__=='__main__': |
| checks=run() |
| p=Path(__file__).resolve().parents[1]/'reports/validation.json' |
| p.write_text(json.dumps(checks,indent=2)) |
| print(json.dumps(checks,indent=2)) |
|
|
