Fix memory leak in radTSend::OutRelaxResultsInfo#3
Conversation
TotOutArray was not properly deallocated.
| double *t = TotOutArray; | ||
| std::vector<double> TotOutArray; | ||
| TotOutArray.resize(TotOutElem); | ||
| double *t = TotOutArray.data(); |
There was a problem hiding this comment.
Nice catch, and very good hint how to work with vector data via double* !
In this case, double* RelaxStatusParamArray can probably be directly submitted to MultiDimArrayOfDouble(..); I'll check this further and correct. Thanks!
There was a problem hiding this comment.
Glad to be of help.
If you want to give ownership of an object when calling a function, C++11 makes it possible to give objects like a vector to a function without expensive copying. For example:
void StoreMultiDimArrayOfDouble(std::vector<double> Array, int* Dims, int NumDims);
// And then call the function like this
StoreMultiDimArrayOfDouble(std::move(param), dims, n_dims);Compared to passing a raw pointer with the soft contract that the ownership is given to the function, with this approach there is no ambiguity and much lower risk of mistakes.
If you want to go a little deeper, effectively the same thing could be done, ever so slightly more efficiently, with
void StoreMultiDimArrayOfDouble(std::vector<double>&& Array, int* Dims, int NumDims);This is an "rvalue reference", a new kind of reference that is similar to a normal reference but it also means that the recipient of it is allowed to consume the object, the caller promises that the next method call is either the destructor or assignment operator. (std::move doesn't actually move, what it does is to return an rvalue reference to its parameter. Code that uses std::move should make sure to not use the moved object afterwards)
…crash
rkleaf_recompress assumed the stacked rank k_in <= min(m, n), so its QR
step used dorgqr(m, k_in, k_in) which fails ("Parameter 2 incorrect")
when k_in > m. This happens after rk-into-rk stacking: the combined
rank kw = kt_c + kt_inc can exceed the leaf row/col count (e.g., a
24x24 leaf accumulating two rank-20 increments -> kw=40 > 24).
Symptom: nx=5 leaf=10 hex cube (750 DOF) returned -4.002 (ERR_LAPACK)
with "Intel oneMKL ERROR: Parameter 2 was incorrect on entry to DORGQR."
Fix: use economy QR dimensions
ru = min(m, k_in) columns of Q_U / rows of R_U
rv = min(n, k_in) columns of Q_V / rows of R_V
S = R_U R_V^T is ru x rv (not k_in x k_in)
SVD of S -> p = min(ru, rv); k_new bounded by p <= min(m, n)
This is the mathematically correct economy SVD recompression for an
over-complete rk representation -- an m x n matrix has rank <= min(m,n),
so any kt > min(m,n) is redundant and the recompression collapses it.
Results (leaf=10, mu_r=10, hex cube):
nx=3 (162 DOF) 2.5e-14 OK
nx=4 (384 DOF) 5.5e-14 OK
nx=5 (750 DOF) 1.1e-13 OK (was -4.002 LAPACK error)
nx=6 (1296 DOF) 1.0e-13 OK
nx=7 (2058 DOF) 1.8e-13 OK
The DORGQR warning is gone. All synthetic regression unchanged at
machine precision (dense 0-5, rk, rkrk, rkdeep, radia_exact).
Resolves open item ochubar#3 from project_hlu_phase4_complete memory
(Intel MKL DORGQR warning) and the nx=5 leaf=10 failure. Phase 5
full pivoting still pending for genuinely non-pivot-amenable matrices,
but the recompression path is now robust to rank overflow.
Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
…e + BH-curve Picard
First step of the nonlinear extension (user: "make hdiv-type radia work for nonlinear too").
VERIFY-FIRST formulation result: the demag-factor eigenvalue framing A_eig = (1/chi)M_mass - N is NOT
the applied-field system. The physical applied-field weak form is
A+ m = M_mass h_ext , A+ = (1/chi) M_mass + N (PLUS N)
(M = chi(H_ext + H_demag), H_demag,weak = -N m). Probed both signs on a tet sphere: A+ reproduces the
analytic M/H = chi/(1+chi D) to <=2.5% for mu_r<=100; the MINUS system gives nonsense (negative/divergent).
Nonlinear (secant-chi Picard with a saturating BH curve) CONVERGES (15-25 iters) and SATURATES
(M -> M_sat). The ochubar#3 near-field correction closes the accuracy gap: the centroid-monopole under-
estimates D (0.315 vs 1/3) -> high-chi M ~12% off analytic; folding the near correction into N raises
D to 0.328 -> nonlinear M within ~2% of the analytic uniform-sphere answer
(H0=0.1: monopole 0.334, near-corrected 0.305, analytic(1/3) 0.299; residual ~2% is RT0/mesh).
- examples/feec_vim/hdiv_demag_tet_nonlinear.py -- applied-field A+ solve + Picard + correct analytic
reference (bisection) + the near-correction comparison.
- tests/feec/test_hdiv_vim_tet_nonlinear.py (feec suite 47/47) -- Picard converges + saturates;
near-corrected nonlinear M closer to analytic (D=1/3) than monopole + within ~8% at strong drive.
- radia-mcp knowledge hdiv_vim("nonlinear") updated: plan -> VALIDATED FOUNDATION + the A+ gotcha +
the next steps (per-element non-uniform field reconstruction + rad.MatSatIsoTab MMM/MSC cross-check
+ Hantila factor-once + the scalable C++ near-correction).
NEXT: per-element (non-uniform body) chi_e(H_e) via RT0->cell field reconstruction, cross-checked
against rad.Solve + rad.MatSatIsoTab on a cube/C-yoke.
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
The earlier sessions found simple per-element Picard / Hantila could not
robustly solve the nonlinear (BH-curve / saturation) HDiv-VIM demag --
NaN / wrong root at deep saturation. This adds the robust solver:
damped Newton-Raphson, the Newton counterpart of Radia's MMM/MSC
newton_damping path (the outer tangent + damping machinery is SHARED;
only the demag operator N = B^T G B differs).
solve_nonlinear_newton (examples/feec_vim/hdiv_demag_tet_nonlinear.py)
needs three ingredients, each isolated verify-first:
(a) the CONSISTENT TENSOR tangent dM/dH = chi_diff Hhat(x)Hhat
+ chi_sec (I - Hhat(x)Hhat) -- the scalar chi_diff*I stalls at
moderate drive (chi_sec != chi_diff);
(b) the ochubar#3 near-field correction on N -- without it the per-element
Newton converges to a WRONG root (M~0.09 vs 0.30 on the sphere at
H0=0.1), because the centroid-monopole N is accurate in average
but poor per-element (the scalar Picard hid this);
(c) line-search damping + a scalar-chi Picard warmstart (BH knee).
Results (sphere, chi0=1000): fast (4-11 iters) and accurate (rel 1e-4
to 1e-6 vs the analytic uniform sphere) at deep saturation -- exactly
where per-element Picard/Hantila failed. The ~2% at low drive is the
operator demag-D systematic (0.328 vs 1/3), shared with Picard, not a
Newton failure. Honest limitation: the BH knee is stiff (Newton
converges to the correct answer but slowly; scalar Picard handles it).
Golden-locked: tests/feec/test_hdiv_vim_tet_newton.py (3 tests; feec
suite 50/50). Knowledge updated: hdiv_vim("nonlinear") now records the
SOLVED status + ingredients + the remaining non-uniform / C++ next steps.
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…nest open boundary)
Cross-checked the nonlinear HDiv-VIM Newton against Radia MMM/MSC on a
CUBE (non-spherical, sharp corners, non-uniform local M). Honest finding,
not a clean win:
- at SATURATION the cube agrees with Radia to ~1% (consistent with the
sphere; M -> Msat is demag-D-independent);
- at MODERATE drive (H0=2e5) there is a ~13% gap, and Radia shows the
cube M is genuinely very non-uniform there (std = 22% of mean).
Diagnosed (NOT convergence, NOT knee): maxit 120==400 give the same
M_avg; the cube monopole+nearcorr D=0.3125 is near-corr-insensitive
(near_factor 2/3/4) vs analytic 1/3. The HDiv-VIM cube is demag-limited
at M_avg ~ H0/D_monopole; Radia's non-uniform solution reaches a lower
effective demag. So the gap is the centroid-MONOPOLE operator's
per-element field accuracy on a sharp-cornered body -- the same operator
accuracy limit as the linear demag, exposed by non-smooth geometry. Fix
= higher-order Wilton Gram (the ochubar#3 production direction), NOT a solver
change.
Recorded in hdiv_vim("nonlinear"/"status") so the next session does not
re-walk it: the nonlinear Newton is production-grade at saturation and on
smooth bodies; the sharp-geometry moderate-drive non-uniform case waits
on the Wilton Gram.
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
… Wilton Gram ochubar#3 (more validation): a prolate spheroid has an ANALYTIC demag factor != 1/3, so it checks the Wilton surface Gram on a non-1/3 smooth body (the sphere/cube fix it at the isotropic 1/3). 2:1 prolate (long axis z), maxh=0.6: monopole N_z = 0.1693 (2.5%) Wilton N_z = 0.1742 (0.3%) analytic = 0.1736 -> the Wilton analytic surface integral is correct for a general curved-surface triangulation, not just 1/3. Golden: tests/feec/test_hdiv_vim_ellipsoid.py. Also VECTORIZE the Wilton Gram: tri_potential now accepts a batch of observation points (r as (M,3) -> (M,)) with full numpy broadcasting over the 3 edges, and wilton_surface_block loops O(n_bf) (one batched potential eval per SOURCE triangle) instead of O(n_bf^2) scalar calls. The wilton_gram golden test went 73s -> 2s (no regression), and the prolate demag is now practical. feec 56/56, suite 105s -> 42s. Knowledge hdiv_vim records ochubar#3 (ellipsoid) + the vectorization; remaining ochubar#3 options noted (real rad.MatSatIsoTab steel TABLE, C-yoke). Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…ine*.md Sugahara: paper drafts are Japanese too (many target IEEJ / domestic venues). - CLAUDE.md "Repository Language" policy: add exception ochubar#3 -- docs paper drafts / outlines (paper_outline*.md) MAY be Japanese, header marker `> 論文ドラフト (paper draft)`. Reference/public docs (README, api.md, theory/method/showcase) and commit messages stay English. - paper_outline_sheet_metal.md (Paper 1) + paper_outline.md (Paper 2): translated to Japanese; code identifiers, golden test names, numbers, references, section numbers kept verbatim. Also fixed the leftover Bambu->DIY gantry references and marked the printable-former deliverable as shipped (--former-stl/--former-output). Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
TotOutArraywas not properly deallocated. This PR changes the code to usestd::vector, which is a more idiomatic way to handle dynamically sized arrays, and it automatically manages the memory.If the memory allocation fails, an exception is thrown. (But: C++ code is almost always written with the assumption that memory allocations do not fail. It is very hard to write programs that use dynamic memory allocation and that work when memory allocations fail. It is more common to avoid dynamic memory allocation altogether if that is a concern.)
The
.resize()call will zero-initialize the vector, which carries a performance penalty. This performance penalty will be tiny if the vector is small enough to fit on a cache page aka 4K. If you worry about this, it is possible to useunique_ptrinstead ofvector.I found this with Address Sanitizer. The change removes the warning and seems to work. The change is tested on Linux/Clang only.