1// Headless validation of the solver — runs the same MATLAB script the
2// browser worker runs, in Node, against the installed numbl. Each solve
3// fills matlab/solve_template.m with the parameters (mirroring
4// buildSolveScript in src/engine/engine.ts), runs it standalone (as the
5// browser does in a fresh session), and reads result.json back from the
6// VFS. The VFS is shared across solves, standing in for numbl/browser's
7// IndexedDB-persisted /system, so `mip load --install surfacefun` only
8// downloads once. Node has no synchronous XMLHttpRequest, so
9// websave/webread are shimmed with curl (responses cached in .cache/ keyed
10// by URL, so repeat runs are offline).
11//
12// npm run engine-test
14import {
15 executeCode,
16 VirtualFileSystem,
17 BrowserFileIOAdapter,
18 BrowserSystemAdapter,
19} from 'numbl'
20import {unzipSync} from 'fflate'
21import {execFileSync} from 'node:child_process'
22import crypto from 'node:crypto'
23import fs from 'node:fs'
24import path from 'node:path'
25import {fileURLToPath} from 'node:url'
27const root = path.dirname(path.dirname(fileURLToPath(import.meta.url)))
28const cacheDir = path.join(root, '.cache')
30const MIP_MHL_URL =
31 'https://github.com/mip-org/mip-core/releases/download/mip-numbl/mip-numbl-any.mhl'
32const MIP_SYSTEM_PREFIX = '/system/mip/packages/gh/mip-org/core/mip/'
33// In the browser this is passed by numbl/browser's session worker; executeCode
34// scans searchPaths directories since the same numbl change.
35const MIP_SEARCH_PATH = MIP_SYSTEM_PREFIX + 'mip'
37function curlCached(url) {
38 fs.mkdirSync(cacheDir, {recursive: true})
39 const key = crypto.createHash('sha1').update(url).digest('hex').slice(0, 16)
40 const cached = path.join(cacheDir, key)
41 if (!fs.existsSync(cached)) {
42 console.log(`fetching ${url}`)
43 execFileSync('curl', ['-sfL', '-o', cached, url], {stdio: 'inherit'})
44 }
45 return fs.readFileSync(cached)
46}
48// numbl's BrowserFileIOAdapter implements websave/webread with synchronous
49// XHR (fine in a web worker, absent in Node); override with curl.
50class NodeFileIOAdapter extends BrowserFileIOAdapter {
51 constructor(vfs) {
52 super(vfs)
53 this.nodeVfs = vfs
54 }
55 websave(url, filename) {
56 this.nodeVfs.writeFile(this.nodeVfs.normalizePath(filename), new Uint8Array(curlCached(url)))
57 }
58 webread(url) {
59 return curlCached(url).toString('utf8')
60 }
61}
63const template = fs.readFileSync(path.join(root, 'matlab', 'solve_template.m'), 'utf8')
65// Mirrors buildSolveScript in src/engine/engine.ts.
66const buildSolveScript = params => {
67 const fills = {
68 MESHFILE: 'mesh.msh',
69 PDE: params.pde,
70 F_EXPR: params.f,
71 C_EXPR: params.pde === 'helmholtz' ? params.c : '0',
72 ORDER: String(params.p),
73 CLOSED: params.closed ? 'true' : 'false',
74 }
75 return template.replace(/\{\{(\w+)\}\}/g, (token, key) => fills[key] ?? token)
76}
78async function main() {
79 const vfs = new VirtualFileSystem()
80 const enc = new TextEncoder()
82 // Bootstrap mip into the system VFS, as the browser worker does.
83 const mipEntries = unzipSync(new Uint8Array(curlCached(MIP_MHL_URL)))
84 let nMip = 0
85 for (const [name, content] of Object.entries(mipEntries)) {
86 if (name.endsWith('/')) continue
87 vfs.writeFile(MIP_SYSTEM_PREFIX + name, content)
88 nMip++
89 }
90 console.log(`mip core: ${nMip} files into VFS`)
92 vfs.writeFile(
93 '/project/mesh.msh',
94 fs.readFileSync(path.join(root, 'public', 'samples', 'sphere.msh'))
95 )
96 vfs.setCwd('/project')
98 const decoder = new TextDecoder()
100 const solve = params => {
101 const script = buildSolveScript(params)
102 vfs.writeFile('/project/solve_pde.m', enc.encode(script))
103 vfs.writeFile('/project/result.json', enc.encode('')) // no stale reads
104 const t = Date.now()
105 executeCode(
106 script,
107 {
108 onOutput: text => process.stdout.write(`[numbl] ${text}`),
109 onDrawnow: () => {},
110 displayResults: false,
111 maxIterations: 1e9,
112 optimization: '1',
113 fileIO: new NodeFileIOAdapter(vfs),
114 system: new BrowserSystemAdapter(vfs),
115 },
116 [{name: 'solve_pde.m', source: script}],
117 vfs.normalizePath('/project/solve_pde.m'),
118 [MIP_SEARCH_PATH]
119 )
120 const result = JSON.parse(decoder.decode(vfs.readFile('/project/result.json')))
121 console.log(`solve [${params.pde}] in ${(Date.now() - t) / 1000}s`)
122 return result
123 }
125 // 1. Poisson on the closed sphere
126 let d = solve({pde: 'poisson', f: 'x.*y.*z', c: '', p: 6, closed: true})
127 console.log(` npatches=${d.npatches} n=${d.n} u in [${d.umin.toFixed(6)}, ${d.umax.toFixed(6)}]`)
128 if (d.npatches !== 216 || d.n !== 7 || d.ptype !== 'quad')
129 throw new Error('unexpected solution shape')
130 if (!isFinite(d.umin) || !isFinite(d.umax) || d.umin === d.umax)
131 throw new Error('degenerate solution values')
133 // Eigenfunction check: x*y*z is a degree-3 solid harmonic, so on the unit
134 // sphere lap_S (x*y*z) = -12 * (x*y*z). Solving with f = -12*x*y*z must
135 // reproduce u = x*y*z, whose max on the sphere is 1/(3*sqrt(3)).
136 d = solve({pde: 'poisson', f: '-12*(x.*y.*z)', c: '', p: 8, closed: true})
137 const expected = 1 / (3 * Math.sqrt(3))
138 console.log(` eigencheck: umax=${d.umax.toFixed(6)} expected~${expected.toFixed(6)}`)
139 if (Math.abs(d.umax - expected) > 0.01) throw new Error('eigenfunction check failed')
141 // 2. Helmholtz with a variable coefficient
142 d = solve({pde: 'helmholtz', f: '1 + 0*x', c: '100*(1 - z)', p: 6, closed: true})
143 console.log(` u in [${d.umin.toFixed(6)}, ${d.umax.toFixed(6)}]`)
145 // 3. Bad expression errors out of the run (the host surfaces the message)
146 let err = null
147 try {
148 solve({pde: 'poisson', f: 'this is not matlab', c: '', p: 4, closed: true})
149 } catch (e) {
150 err = e
151 }
152 if (!err) throw new Error('expected an error for bad expression')
153 console.log(` error path OK: ${String(err.message).slice(0, 100)}`)
155 // 4. A later solve is unaffected (fresh run per solve)
156 d = solve({pde: 'poisson', f: 'x', c: '', p: 4, closed: true})
157 if (d.type !== 'solution') throw new Error('solve after error failed')
159 // 5. Triangle mesh: the same eigenfunction check on an icosahedral sphere.
160 // Flat triangles hug the sphere worse than the cubed sphere's bilinear
161 // quads, so geometry error dominates: umax is ~0.0107 low at this
162 // subdivision level and shrinks 4x per level (O(h^2)) — hence the wider
163 // tolerance.
164 vfs.writeFile(
165 '/project/mesh.msh',
166 fs.readFileSync(path.join(root, 'public', 'samples', 'sphere-tri.msh'))
167 )
168 d = solve({pde: 'poisson', f: '-12*(x.*y.*z)', c: '', p: 6, closed: true})
169 console.log(` tri eigencheck: umax=${d.umax.toFixed(6)} expected~${expected.toFixed(6)}`)
170 if (d.ptype !== 'tri' || d.npatches !== 320 || d.n !== 7)
171 throw new Error('unexpected triangle solution shape')
172 if (Math.abs(d.umax - expected) > 0.015)
173 throw new Error('triangle eigenfunction check failed')
175 // Helmholtz with a variable coefficient also works on triangle patches
176 d = solve({pde: 'helmholtz', f: '1 + 0*x', c: '100*(1 - z)', p: 6, closed: true})
177 console.log(` tri helmholtz: u in [${d.umin.toFixed(6)}, ${d.umax.toFixed(6)}]`)
178 if (d.ptype !== 'tri' || !isFinite(d.umin) || d.umin === d.umax)
179 throw new Error('triangle helmholtz solve failed')
181 console.log('engine-test: all checks passed')
182}
184main().catch(err => {
185 console.error(err)
186 process.exit(1)
187})