From 179073924bc203c0938d4623f3e9042e62512611 Mon Sep 17 00:00:00 2001 From: John Date: Mon, 20 Jul 2026 16:00:29 -0400 Subject: [PATCH 1/5] heat source calibration --- .gitignore | 1 + .pre-commit-config.yaml | 32 +- Allwmake | 54 + README.md | 15 +- .../meltPoolDimensions/meltPoolDimensions.C | 2 +- .../absorptionModels/Kelly/KellyAbsorption.C | 17 +- .../calibrateHeatSourcecpython-310.pyc | Bin 0 -> 56901 bytes bin/calibrateHeatSource | 2617 +++++++++++++++++ etc/bashrc | 20 + .../heatSources/nLightAFX-1000.cfg | 35 +- requirements.txt | 10 + tutorials/heatSourceCalibration/README.md | 106 + tutorials/heatSourceCalibration/config.yml | 33 + .../heatSourceCalibration/experiments.yml | 25 + tutorials/heatSourceCalibration/template/0/T | 50 + tutorials/heatSourceCalibration/template/0/U | 45 + .../heatSourceCalibration/template/0/p_rgh | 50 + .../heatSourceCalibration/template/Allclean | 7 + .../heatSourceCalibration/template/Allrun | 13 + .../heatSourceCalibration/template/README.md | 29 + .../template/constant/dynamicMeshDict | 39 + .../heatSourceCalibration/template/constant/g | 22 + .../template/constant/heatSourceDict | 66 + .../template/constant/scanPath | 3 + .../template/constant/transportProperties | 18 + .../template/system/blockMeshDict | 85 + .../template/system/controlDict | 69 + .../template/system/decomposeParDict | 27 + .../template/system/fvSchemes | 50 + .../template/system/fvSolution | 76 + tutorials/nLightAFX/README.md | 15 +- tutorials/nLightAFX/constant/heatSourceDict | 8 +- 32 files changed, 3595 insertions(+), 44 deletions(-) create mode 100644 bin/__pycache__/calibrateHeatSourcecpython-310.pyc create mode 100755 bin/calibrateHeatSource rename tutorials/nLightAFX/constant/nLightAFX.cfg => etc/heatSources/nLightAFX-1000.cfg (84%) create mode 100644 requirements.txt create mode 100644 tutorials/heatSourceCalibration/README.md create mode 100644 tutorials/heatSourceCalibration/config.yml create mode 100644 tutorials/heatSourceCalibration/experiments.yml create mode 100644 tutorials/heatSourceCalibration/template/0/T create mode 100644 tutorials/heatSourceCalibration/template/0/U create mode 100644 tutorials/heatSourceCalibration/template/0/p_rgh create mode 100755 tutorials/heatSourceCalibration/template/Allclean create mode 100755 tutorials/heatSourceCalibration/template/Allrun create mode 100644 tutorials/heatSourceCalibration/template/README.md create mode 100644 tutorials/heatSourceCalibration/template/constant/dynamicMeshDict create mode 100644 tutorials/heatSourceCalibration/template/constant/g create mode 100644 tutorials/heatSourceCalibration/template/constant/heatSourceDict create mode 100644 tutorials/heatSourceCalibration/template/constant/scanPath create mode 100644 tutorials/heatSourceCalibration/template/constant/transportProperties create mode 100644 tutorials/heatSourceCalibration/template/system/blockMeshDict create mode 100644 tutorials/heatSourceCalibration/template/system/controlDict create mode 100644 tutorials/heatSourceCalibration/template/system/decomposeParDict create mode 100644 tutorials/heatSourceCalibration/template/system/fvSchemes create mode 100644 tutorials/heatSourceCalibration/template/system/fvSolution diff --git a/.gitignore b/.gitignore index 60ab4ee4..4c855fce 100644 --- a/.gitignore +++ b/.gitignore @@ -20,6 +20,7 @@ # Derived files lex.yy.c +.venv/ # Corefiles core diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index f87ebde5..12440dee 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -2,27 +2,45 @@ minimum_pre_commit_version: "3.7.0" repos: - repo: https://github.com/pre-commit/pre-commit-hooks - rev: v6.0.0 + rev: v3.4.0 hooks: - id: trailing-whitespace - files: &openfoam_cpp_files '\.(C|H|I|cc|cpp|cxx|hh|hpp|hxx)$' - exclude: &openfoam_excludes '(^|/)(lnInclude|ThirdParty|platforms|build|processor[0-9]+|postProcessing)/' + files: &python_files '(\.pyi?$|(^|/)bin/calibrateHeatSource$)' - id: end-of-file-fixer - files: *openfoam_cpp_files - exclude: *openfoam_excludes + files: *python_files - id: mixed-line-ending args: [--fix=lf] - files: *openfoam_cpp_files - exclude: *openfoam_excludes + files: &openfoam_cpp_files '\.(C|H|I|cc|cpp|cxx|hh|hpp|hxx)$' + exclude: &openfoam_excludes '(^|/)(lnInclude|ThirdParty|platforms|build|processor[0-9]+|postProcessing)/' - id: check-merge-conflict - id: check-yaml + - repo: https://github.com/psf/black + rev: 24.3.0 + hooks: + - id: black + args: [--line-length=80] + files: *python_files + + - repo: https://github.com/codespell-project/codespell + rev: v2.2.6 + hooks: + - id: codespell + args: ["-L", "mater,nd,coefficent,ot"] + files: *python_files + - repo: local hooks: + - id: python-max-80-columns + name: "Python style: max 80 characters" + language: pygrep + entry: '^.{81,}$' + files: *python_files + - id: openfoam-no-tabs name: "OpenFOAM style: no tabs" language: pygrep diff --git a/Allwmake b/Allwmake index 830f8594..591905b9 100755 --- a/Allwmake +++ b/Allwmake @@ -27,3 +27,57 @@ $ADDITIVEFOAM_BUILD_FLAGS \ applications/Allwmake $targetType $* #------------------------------------------------------------------------------ +# Build the Python environment + +pythonCommand= + +if command -v python3 >/dev/null 2>&1 +then + pythonCommand=python3 +elif command -v python >/dev/null 2>&1 && \ + python -c 'import sys; raise SystemExit(sys.version_info.major != 3)' +then + pythonCommand=python +fi + +if [ -z "$pythonCommand" ] +then + echo "WARNING: Python 3 was not found; skipping the Python environment." >&2 +elif ! "$pythonCommand" -c \ + 'import sys; raise SystemExit(sys.version_info < (3, 10))' +then + echo "WARNING: Python 3.10 or newer is required; skipping the Python environment." >&2 +elif [ ! -f "$ADDITIVEFOAM_PROJECT_DIR/requirements.txt" ] +then + echo "WARNING: requirements.txt was not found; skipping the Python environment." >&2 +else + additivefoamVenv="$ADDITIVEFOAM_PROJECT_DIR/.venv" + + if [ ! -x "$additivefoamVenv/bin/python" ] || \ + [ ! -r "$additivefoamVenv/bin/activate" ] + then + echo "Creating the AdditiveFOAM Python environment" + "$pythonCommand" -m venv "$additivefoamVenv" || exit 1 + fi + + if ! "$additivefoamVenv/bin/python" -c \ + 'import sys; raise SystemExit(sys.version_info < (3, 10))' + then + echo "ERROR: $additivefoamVenv uses Python older than 3.10." >&2 + echo " Remove .venv and run Allwmake again." >&2 + exit 1 + fi + + echo "Installing AdditiveFOAM Python dependencies" + "$additivefoamVenv/bin/python" -m pip install \ + -r "$ADDITIVEFOAM_PROJECT_DIR/requirements.txt" || exit 1 + "$additivefoamVenv/bin/python" -m pip check || exit 1 + + echo "Source $ADDITIVEFOAM_PROJECT_DIR/etc/bashrc to activate the environment" + + unset additivefoamVenv +fi + +unset pythonCommand + +#------------------------------------------------------------------------------ diff --git a/README.md b/README.md index 2d0272c6..af14a54d 100644 --- a/README.md +++ b/README.md @@ -13,6 +13,8 @@ The documentation for `AdditiveFOAM` is hosted on [GitHub Pages](https://ornl.gi |-----------------------------------------------------------|------------------------------------------| | [solver](applications/solvers/additiveFoam) | Development version of the solver | | [tutorials](tutorials) | Tutorial cases | +| [calibration utility](bin/calibrateHeatSource) | Heat-source calibration command | +| [calibration tutorial](tutorials/heatSourceCalibration) | Projected heat-source calibration workflow | ## Installation [![OpenFOAM-14](https://img.shields.io/badge/OpenFOAM-14-blue.svg)](https://github.com/OpenFOAM/OpenFOAM-14) @@ -40,10 +42,21 @@ script from the repository root: ```sh ./Allwmake +source etc/bashrc +``` + +The second `source` activates the newly created Python environment. The master +build script uses Python 3.10 or newer to create `.venv` and install +the dependencies in `requirements.txt`. If Python is unavailable, the compiled +applications are still built and `Allwmake` prints a warning. + +```sh +calibrateHeatSource --help ``` For regular use, source both environments in each new shell or add them to your -shell startup file: +shell startup file. The AdditiveFOAM environment activates `.venv` when it is +available: ```sh source /path/to/OpenFOAM-14/etc/bashrc diff --git a/applications/solvers/additiveFoam/functionObjects/meltPoolDimensions/meltPoolDimensions.C b/applications/solvers/additiveFoam/functionObjects/meltPoolDimensions/meltPoolDimensions.C index 855062ec..a8a5b692 100644 --- a/applications/solvers/additiveFoam/functionObjects/meltPoolDimensions/meltPoolDimensions.C +++ b/applications/solvers/additiveFoam/functionObjects/meltPoolDimensions/meltPoolDimensions.C @@ -247,7 +247,7 @@ bool Foam::functionObjects::meltPoolDimensions::execute() // physical boundary : take face point if above iso value const vectorField& Cf = mesh_.Cf().boundaryField()[patchi]; - const scalarField& pif(TPf.patchInternalField()); + const scalarField pif(TPf.patchInternalField()); forAll(faceCells, facei) { diff --git a/applications/solvers/additiveFoam/movingHeatSource/absorptionModels/Kelly/KellyAbsorption.C b/applications/solvers/additiveFoam/movingHeatSource/absorptionModels/Kelly/KellyAbsorption.C index 9d83eb56..05ff6f84 100644 --- a/applications/solvers/additiveFoam/movingHeatSource/absorptionModels/Kelly/KellyAbsorption.C +++ b/applications/solvers/additiveFoam/movingHeatSource/absorptionModels/Kelly/KellyAbsorption.C @@ -54,10 +54,14 @@ Foam::absorptionModels::Kelly::Kelly etaMin_("etaMin", dimless, absorptionModelCoeffs_), aspectRatioSwitch_ ( - absorptionModelCoeffs_.lookupOrDefault + max ( - "aspectRatioSwitch", - 1.0 + absorptionModelCoeffs_.lookupOrDefault + ( + "aspectRatioSwitch", + 1.0 + ), + scalar(1e-10) ) ) { @@ -116,12 +120,15 @@ bool Foam::absorptionModels::Kelly::read() absorptionModelCoeffs_.lookup("geometry") >> geometry_; absorptionModelCoeffs_.lookup("eta0") >> eta0_; absorptionModelCoeffs_.lookup("etaMin") >> etaMin_; - aspectRatioSwitch_ = + aspectRatioSwitch_ = max + ( absorptionModelCoeffs_.lookupOrDefault ( "aspectRatioSwitch", 1.0 - ); + ), + scalar(1e-10) + ); return true; } diff --git a/bin/__pycache__/calibrateHeatSourcecpython-310.pyc b/bin/__pycache__/calibrateHeatSourcecpython-310.pyc new file mode 100644 index 0000000000000000000000000000000000000000..da12bd018118d9fddacb41d9d226b85edb63044a GIT binary patch literal 56901 zcmb5X34C1FecwC#j+Gz??junPA(7xB$+jel6jw=(1Ti8l$B>xu5O)Xy94x#u6p8U5 zwnf>NQ!{oJw@oCapc5xiv$jc^rtdXP)7IJ7=B4dTx+G4cFi8_TPGXy}#P|LE&z->x zNZTp!!NED_o_p?D|L4D-`JSF!gr7fbe)G)7o`^<1t3>Bdg2bbdNVN8?#JtEv#EW{d zYIHU>5sNAwuf`|h_CGO^u>Z-4r2S7#r0josB5nUO6B+xToygk%+(gd)=O^;~PgD!D zJrg}vXR_Km+dI)~@l>^Mwtu4E;_2$Z?BK+p#WU44vqKX@7LQlgP7GUJ>n29*|N7Yt z6C1;4nV32KPuxHIz{CTV?yJ6S_U#jI zw|IZ`9kUNkJUIK%#6zy=+hPP~hFpLg1ut`#sOA5PvuE0q?9gP5hgP4|?akD)ILa zU*pYr^Ta2J4|(T3pZIgx8mYDX4SS1TgPNbuzR&X3kvHNkdEY|bDdOwB4|pHsjc+Et z!Mp6eM7%_NqxT{2!^Gd~eZ+g29+%16W8|GCzS;XW@8iTLy-#@G&U3HQ<9&zs zomZldMZE9wKKW9_`{d==<>*AEQn?(NIKvxTyno{TQ{I>&zSa94?-k-_i5L01&HG;O zpV6Xe;@A3Z7`yJ-T^1j9U zJMW9+y}$AS?@Q-mpNRQ~EAPJ?VRyt&*x&ZIMI(`NBo>L}4wtLbr~PuhQmmh?6c_yY znaX6n;uXEhLj7#fo37RU>C+eL)AMu13pMgjFBT7a-gJHXh03v~4?a~qTPfFzwfPJF zWTiHi%N?8di0 zO}EF2C#GjFxK?TCY=|@*^w#E5lU1v0PRl-H%5=DiV|JOP+RAi{<*w%5j!&MQUU*`T zjxNm8x_O_>o+?AmxSYR0Q&d-Ka=tq6*G6MO_T<6R;iKaxk3JKmi9UJs*hxLFR_fHi z^WNiSzdYra7tWrjFIFqTnyE_tM0s|hTCwzrvz1Dn7ZTHRlV?ZMLGp>&@>C_r9w$8H zSIXxU4J!x+7%^TvGCz5;LWJ%FsS^ujrXon5ET5+Koc(iOvd)*Ptoa><*Xs?_SEF@Hp_^to3`nwkA5MU#%Jm6=e3&JU5E}w{X128 z=FD`pvcQDi)V~amt$$AvsYhm_%dsXss%56=Q9a%SkGUNCcpv>uEGL`MtI>s_smQwl zE*?XZh+K|GDYqw5+j=n?1qMYU?-3*-Z!|U@ME%iJ5Syp6z4 z*eq({PHOI2KXD@Metr{)dZfu5&O~Qooc^5pSHRckD`!4^U9%LcCuWk(NHfN{pWwVt zT}eji^Gc7^Kl*&0D*87Nn-JnNoHgf)$ z9%|8aj$N(%LS<#_-pSLN#$t$Z8;;(7!+(TCkTfV8b3m5)>I=ZQ%G?Xn{`{Oy9Zm%f z2X;2cL9r2gs5qMSx6>?tlpx4)JOME|(fn4N@-??LbI0LcbhcPF;=6q?ZyegTAhhEG50nk;ti@rBpN3 zOui8HpX9x4GsXOYy1f`Zf1sXeCcWJHNHgK(FUObCmM&PjhjgZ%w47ec>AS4i&QfW+#jeIa8mZ@Aj@0v)WB!Zv!ppSDTYE01TRVGRj(F?V zlV4YCM6q*8{~z;wy`F=lCSHnIog1ej-lk7QueyF0nuXfyl-g{y6zaXrf&h4(d}@gM zKf|*v_N?zC5yot*S7giz8neq$|1&(j&7SsK-M4$&t}+HoJU zdUsiw-umG38t?9A@6}lC&nUB9WxUa5!FsT&FVxp#Z{OqHtGC0{4sWNW+}FpvT}pYo zU(GD_)rT0#JsO9l{`x>O>+Rhb@$TD5y4mmDzc12EH#5qoE-g{nHUE)$bj6wTGSh=(vA3%XTywB(>^+?~~

9q=%>PiOw*{OrR1 z5@qciNk#gju^X{y>PBO$9W`zbkL@!Qv$mf|>z&5vVFQsDK&pj^`IYk*rrG+%+48~y zM{uzQf-N-bq2l$}j-oGkOFQ!`2s5Wby@jP@Cc${h(^GSe9ZyZy=*m>9LeRag?*Q}s z#fOSEzGMEosI@WR5*;=B&nmyy<+ncKzq=!|^*#R^$nFx)Jo!q&mb#$q^f`$Nl%%mk1(}!AJcG{xx=V_NDs;*B4OiBK=`<>x+%I z{XcOAI*XloKO7i!Aj0d(UVjoB_JX>z$ZZx*sdQ-Gs zX*}}v`dXc?7;(pP=`8A&Q?P}F18vnJdTJ$T%hJlSzLus*jc@LJ@oT@pH$?7eE^4#o zYPI6K4g>=WqPLuchzK&TMTb!*n)^@&RM5(fa zH1Ks6WT}xE+jX#7^)Ji?edj8bg;MRpq-fML7pljCL6Pb7XTDmRuFX3HG*g`ho`Sjg zPpNTPs2X-9KXq z>T_a-88+tF6$oSWg--fUD_Bnu#H*D#|C=pVtJH&dd2TU?vp5`@8}p$m80YG2 zxjuQ;&cz^cn!z3&SRodIf$;qDO09D%=$~~bQ9H-sz&gn6NcAr%wy2Ur5wL1Noe)PZuG`ee^cmFPKNK%A}iAy* zf;V;qYVDX`o~?92dLcie$_({iAkqW~i>LyXjebfLDkRZV1jru&F~4I8lFEyNM1#~X zL#1AgepJ+IlAMGgu_-SN5}Po@AF3TrMZB44CPiR@K1N^V4bF((dWc+cQe?(Q6G0AQ zlcvF;sRelhTs6li8J?T>0j2Q2yi#X=zJu++`OChqUYMP$`JYh#2Tq7-c7m5$Kra}6 zD&*0Xj)()<%I_KrA%i;g3>n-+#7_KN>%`Z|cR9K-(sAsAgq zJNy3&MQZxOrby#XdRwb>&d{weKG53Of)`}l6KDG{NP+KF<~(DsK&yY7J`F6LocFz2 z5C>Elu%a$QWS!dhd~7Mblo2;{qO7m&{`9DAqY#3izmaAT3PF_*DR}~0y2Jy!9g8!P5 z0-7)kzF+b0)uS}yRpJ!zmzB!Qg=eULRVBVhfxa0eZ6@8UZBwqnu5i{mhWw9tRuhvZ z!8|2m4pwQ_m@|Vau^Xuvv&kM7uaxtvzA$R^+MUL@vj~ zEMbhrmxMv!Qs+b~Y0A~FP{1}y zG3x(>3NI@-qrlFb9vh3+2JSv>MSx6d3%7m>M3lDXBbNkqU4rl6cnl!&2t(l%z-aUb znyDq2ZFAA{EG{o{fkvupmjQ%;!q}^Ey$2^GsU>9mjR@rcd$)SF>pviPXGpJka{ z$^rCp&DaP7<)zL={6o#`QpUcpx0$goOt>20VkMhUzp9xW_#l=8??O3q@HE7=a!;C> zeUW;qo_24#yvs4feOUxtxht`yJYqfLN2pyYwJ(BjuQMl24H=k>iOUiei(%U1wg9sx z@VTygKBRks=nJc6)*#M*9NFZK4_;RD?Endf&;J=^ z*q~aVJ#GVdS8D)AzWT@fpQZHcQ6snU@ezQI7|B%W&qL&H|dMzu%Hs!xi>|YWDDcHecg}knW<3NLa26;mSOX4nH3Zm0S z9gATz4P&&sx53C)8+{1NNNpeWCeoP};`3G}k#fh-o4KR_TtB~t67ggZgsY!Zs$anF ztwGYgQU5r~o6L@TH+aF7)P^pymX2&59>db*w%+p~L$U1F%PfRRBt zrgFdHtNKznqa&j3Ao>oqKSe1Wpc%M(M1w@Va;YAq9O>T$b+yk_+S)Zn1Q;8+DzGKc z6618_{DRRf_IB*GIBgWLh%H3`UfR&{@$0dD`>x0K@29mtu5XFjs!A}WFOMetU*w&*AJHB!VlHj?85fNI@1iwB$JFlNZ4%mlFbf0?u}71C(u72hn?w?Ie9kG=nrO zbV+7{ED|Y5SPkglyXK}P>N)~WbH%+sA4uCtER9u0_x{G zGy-A+^)i1kll}<#W{rIyaozu29$e2r4Lye>=v3Dl(_FU4LwvG& z%J{#j;I|a~wu0Xw2vQdo5a?9=k0_A_z4qW>gyR$xZ1J3(ud*O&P7ZOa0o$P07-c65 zEp@(hNaHt5-6?jWMsQ^ehuFLRe^Wz@xcdOHt`#M`Q6m`FF8vE4ux8l}wS1isd|ARB z4$l{fd9jV+PQ}6Y5bK1oWIL9!SemC9i)A-Prcwwj#C*?r`7VLQ!}S=)W{>EXmcXJJ zoswQm-&71Om;u5;kOB#g61Fo(*_p|u)XN;<8zaj}W-#4MvF0~{pQe{n-kN5H-8%Fd z^(|##i%ZDh^OfOdlKge{bflT3H{h#?I-WiMPFwNo#~WKth9e4t32g>%9=5OeOeW-e zs2D_dvln3tz>cm5y;7ZUlG-B(dV)+1uGDc5qhL)PMiBIgUq9|62=(jJ6^Mk|B02=K zFwJ4B`twzLC*^o2UF%p!n_!YrN4FEPqpe3gL7R=3!vtSOhTFv*1p?heX~!I3M@hi| z-b8&qq?>FtadO$*+XlQt#W_JEr?h5tup7TL8hf)KPWe3m%xNvvm4&8ZxfKS{*w8(4 zXhwughKGv&Nh;S?;b9m5wpCA zLP83VCpOc-G2q_#eEpanIH}?hRsNy^!F=1>L2{Dk_0inQ-U@o#ycc8G2kCH3g0&qR z7Rfk#VKJM6p7x}vJanEDJ&0on=d$KNXGLH9GX*w>pH)l~84TZgx}MhD9-@F?ftV>N zSSjuovhsQ^1Bd*rL~<~ih#SH3rel!|8tj1C5lg&f44AhPZ~332no0iJKh$lw>R~|) z2XiGX7QFRh?7IR&F-}Z`VsI)>^(_ID_1yFiHcR(_9hnS8j=j=9@KU@fqF|IaLGpQ1 zr%`4+7<~6rrQ^>${m!F@PnMoN`pk(Zo*oZkd+u)>dm1WTs!Tig+%Jus)}WzfV_Fo& zj_127^~qhQQP}V&;T?lFREp4q;@piPQW`mU>5{asIvo;=)m zSf3JT;3#T;sZ+7)YkVpwlv^5D^X1vbot>#xwdEJe(^V6HBZ4}9@Z{rAyW1ITZ5xx3 zJBqjP3~P=}WXn=2C@5xCgXG1t)9BQo_|p=iuf~~fYp<2o$IsA|QiuHg1uYw4Y6#Bw z6~9e9G{=n-cP??m|9drv!?=n>Y~H^<;FL2=LW7b9jLNRlB~R%Lf3{Xo8x;;1tt0k3 zgW%8Sj1duN8Y%o2tT)&mYb1HXx$|z;0oY4w@wV0}Jr@=B;zbc(-&0R{b~e@1%Nb9Q zpYz1l(c9};&r}k?GdV9j&L-3lax;0Jr$x?rIq^9lKU0**^Q6FEk5{;ogb+>6^m;w4 z<5SeqORCR_)V;o?l+LD?BK7{|0dG)(bIu)wL$B$rzECusrhRMA#r#X&@LbeerxPl~ zEk(8tLbhg%Y#lil^*`3k#z4PMM9w$)+Iqt|H$cE{d?n$ODn>Vr@7Lh~Gh-gW3=91Q zQjNz9D6ep@TdelB@Z#>U_}x2-cel9j#8GA|7deqQpvq`RRf1lv0Ub*Kqq4u)ct>}- zb8{8z^U_?C)&?>?kQwRmwLV<)f1h^w|3i&S?}Sk%IYDZt2{?lsG8%t&UgsovEW-8vloe6>K-VE0r}5%vh`#KQSA%*)Q63T}L*hx6 z+?m)-Y|6M}Q|3FzHgm@3kYQj~?L<3%|CRQdwni!f9~*IHV`od(q&2v0LBuUj-_-Nl zpRe7sy^$P!zIN~SMv8y;JU?p7!{nTAP&YrHJMhN9pY8nMcmI>~k2LP=7{WH5qf=B8 zU&LK@C>9$dMH5Fk#Z41|joE5?B#2IpuIV7qK1}@%^+*B+rWb-Fx=+lX9hG0w!!IlN zB0-R=T{s=;G9tW3bd)&$#POqO^qoBN^t(<5sEtZnPeg?e=Okj7vj1Q7`p*%J<{e2P z{#DS^f<8ua1bqg^q<7=ia4^(H#Y%j!g^OjS9e|=DVzrc$7sAUmv0KpRXzNmYz7lEy z46B>y&h=k)jy|mUoS}KWaby&b&N;{!5iuE)a3pTzfS(Sk77|E3jJS4u-q`h6258G~ z>ThzrR)+TyHxfw>>lm_)8>zSYQ<9rd0}=_~FpwBnypC@My~2x1LLQ7*>srx$v?tA< zf!|?DJ7%y@N)d{HrLh)yr2qaV2*bDz@|t1-$8irb$lrV(eIobd1Spq011{2DZd1Y{Lk7ikun z*;brGq8A2tZ?k7(fcH)NyXJ%;slH?OD6k=C%#&*iF%iQ0L;KHBGxsC>eo zzlfrbsti&{y3zE8*A%wf{{>aEUgct!b_cOV{u=yU+C!Yb#vcEFtKN(h&?OXKH7>!{ zuQL+ANYynVJmNNGqMHo$5CXD=Z7}GTKOB-_Hr74Hr}D?tm`v*o>cTvv9}jw^MgT=A`MUUhCSE^|!1qX7UfoFEGd5X!Vd|5= zch{fh#|HLyt=-m{wnk@~NCqP;dsoi1mm@+(UXFr{00BDAH8^hudELq8m{_A!07G?` zhCSO>esJ@w%I##6j!;c2HWk?tZR|C!#R{g6WWsB%I^mij2IeKk8FW{Sm5a z$}<_)pJr?(5~|5C*iKl2viTgsPY7ivM=ae+(bGRoO@GyDckBZrk;a~XaJyRyGkj=% zuF|;sC{`3^j-gI>PgT)rT|jU8MIUiRq|sOWhKGx7#2<~xx%xA;-yrBGAlh>)@|2LM8zW_%{Y;gXaL9aK%RgkoNDC}OEp zSBiJtGcWyW%IM9;T~Dppn6&p-YtkV8ObP6vVk2W>GEC+gX+iYuwb4cvy^X2*S@Ndn zP~-zIjU8-cE>3%v1y2>7MyCx^`UX^^+B!0;wnIs46<)uNi5eR6R%crr*JJgylD>ok z@E(K(*ch8=G{_jU7Jv{Wt-0aKf#U#o!@HyvbidyCX1$$SL95)Z&hVlk6G03up*ZKk zYAC+JiL-Cjg}>9%(5B3xvynqM(2HyoW&&|(RP;tu*ZptMw*^PKwl$e6D?*|O@wMS( zBPp1H4+#v-OY^Ofk;GsIWWqG}v5Yd__U0CJR+)=ZGWP@Oz&0h}|0$7`4%=V=+0=m5 z9UFoxGk;F15$E83X$IIqIS&6?Kc(T7u)5VL_++4A{*$zs?Sx2XdXnvA3PzK#%DNQm z1qO9)8v^zW5s?w@B$9v-lq)Y@ha?V$j;#|7E6%-t7ppnG%nAn?+vp8dv@D9d-o?T!bT&fb^akD>bwN%41mZ`I&Eqa- zE4@3&omi-#f-wsQYMK{8^mq^}J4V`nm>AE6zW|%l?w85xRTR%4#1#yLl=WrQ4VugF7S37bu@A2rx9bh#XLr znKA7I{%c4R<0mmD0Ssw~&Q{QdwD~vEHBN3LjiKQf?c!+8t04si7<)&2th^)se^$G! zbRS6;l1a?0wQKZz)k0dU9(J+qw4~ml{z_~n(ByuNioD*RNk#p$)Fab9+5U6^B`v=! zf3~zfKq(2$FqL9=gbOTgaksz{PDl&gCUK$YH2}~?6T~#dSW6>-fl~Ch*OSYsW?GaV zu{2Z`JBv8l5G%(p%&?fpF-$IXK^f|)(ScZ^WTxXwSq|zn)(v_SYR=1nvgc*lu#|(k z%X$SYAoRZnJ&To+d8{_lQkjwciGRM{BT>?FZ!_N#FV*|JexSwx2=<`+A{zvukvBxn z+ShoO@AdO#L`H0-7_?=F(*TUOJSai?`DF%3cmrncrwJ~Q=|*FXO|Pn|w~ z!2bYGY)AUfYk7#!avZn+uGpU_c!Qv^`x%ofV^^hGTI&3=M-4d_Ua>!w`R{c<7Pq@v-~QG!i?C#~)Gn!XphbrW%Rj?5>*IUNi{FXy<$4uout3 z&88i&mi+fq0$k%0YNl;$!4}tvcGLZDA@7(&YKD!??9^oiK-^w7Ntdyeo*_?CxrIfY zQBIZOcU7ebLZ?*m-;@;R$z(Q&qGUs9d%~e0zo$n!h)gZZIP0Um4)9qEr}e0zLN+Nq zU0@E0zk@ioZG%Efz*Q1d{UyC5wUQvWGU@@sIWf2ZmvFf3G+njDx2j5x!0v6nqJMT& zn|+SDQ{dTBXce1DfDXn`tPm46o_I5t7LEer9vnPmKmi!UHbS0*3!z^=Pd^F> z1=pDYhe34;BOfu4L~TQ?Nw$k5XeI?8|91C=|E?ADE0GeD8-f&^^E4l zNwh$drtAVHb5v9G3LiFgKv;>eq}e!GFKdNVTdas@b zNB)RRg5d5{rhco-|4?(H~+sJnHC0xPDcusnq>agaCE zq@-EsO966n-V|Baxi)?k%~czZJ9t0Xh{91RG+F_Nh~w^D9LaV37L>)4VkDaD_S-!c)? zGesf}t;mrLAZEQ?VWgzlYB|x4oFq9=k{%k?S7NNd6bBdjtdx}sAX))l0+!|&rjZ$q zG(&JjLo+d=7cp(G< zwtlZ^DI6mgBQp5aW_@27Ed;T-g&>Wpmxi54%?2+-Jul2*M^Q1_KftB}bpw=+f;a{t z{(~C$4^yQtyvA>mglDcv@LyDLnIK3q`;y1cg_9a&7TYmbcWF*QOJ6BLD+ZY`j^ut0 z4RO3Z1nET=S&gA=Fvk!jXVZei27i7^S_HCgd98(}ph*VKI|$SgQLq8Mv5~0%Pt`oF zXTYc=Lo8Bhn|DpUz)YKr628{nI>J;Ka5f9Ot z^i&rZZS{@y5UstFz6l#c@t^Hq>z5fj%tAx(XJTR*@`RXriB3OD)u*zeCm2Q{WS~CZVWya`N3UriVfLQH^9LI70w}lQ zb_dyg0(7#x{7>q=zOf^v&_d`bI*U95AsxKGRuvAqeNf-hQRkNl zZrKInXWqqA<7di|K}Q!}j&Q7)%1m2X<~2>?KU46Nyu>p2dH%Us6{s`D6NE0 zcZ4#oo{flumhyGT4*HiDpZ_vS=b-y}de;Nk$Sz~Oz==|*=a1-a+8TmW7vr~ zd!xBGhf>4QF>*El`v=7iLIY`tGBV^s((|{8g8##OdQ#x9^`n*AUacP^;Xpld8(XI} z2>|V+W%dLT6M|xzQ;1c~E1X}M2C!vTQCDN}cq}g62&s6< zj|D1TD#r=Tb^(E9%Y*2X1G0~&kc1a!E5%JshgcIKxRg414{0zVgHDc2j+i`hVqp$y zGk^^P1P}*Ec-hW?myw9JGpAev_Ld#2Z@b?~jkJZ#zzbhMkiqq!P{UYr@~lH;#jq(a zRCt6>sL_0oG{c}E3x*>nJqu`gIuxXxjvFpyPvOW1S(5OIrP^rG3)qbL!jN1Z5tow% zpJ|m-)N)tPV84W|Wu1(0WEc2)lnF8RX^v^*_svaUO)%2=HGTUV{=Z||g1&hO3>8LSk$k=1Od5d4pgZ z!PB`K`S?8BCOMOqw-+J9!oIGu<*vv_2q06ceW1~bG4r6v$>qG46i#pGXO+wv?qBkvbUT7sgm4x2CqU#BC5I+mpd ztt4JYEVaf`L*CkslEdCQrMwY%dihR2E6uz){rPC8rxoRmt{j#{P@29n@ zVxeH&&2ZlWj>%9wY0W~AlklKAeRh7{!#_l9_5zD-dTO>DWLncDR*QMWg1AJ_rsZPH zoTN3@<26mQ>U2Yk+wI*F?xe(pxeqSP=-Aj86a^-zwdvnVvkVgT3v(4aM5K%gkH_Dr z)sY7(Pm!vfE@Q3U6XFBq$w}-pr+)p{^z#=7;7K4Ilr91+-fQBnz@zzPUiW`R!LKU# zH3h$};5P_HhuoQ>lLR`tg5wX|Yk5vDF}5%pB%c!X8>HVgH!aq;JBEH-um3*P2R%n* zg(}w)wK@j#wvbcK+RvX=smm&rs+q?ycdmVu9EX?}I>?!yS7F1Zvh6VnlID#BpGsNtyaQCpkprf!`rqOa@C4sh=vA^=wTZ=4T#C21n^sR`4Q8z0^VEDOtXTbvTvCgc|BytX{sQ*q8?#hY16^)$2Y ziCK*%WHm~S*;g}cuAJjz1hMMEHII%$2Pg1|R`Y*I10(qThTDAKec+Az?B}lza8klG z&|KPoEj}pOSvDZCoPanQiNHPqbQYx^KKOv(uRuGov(r;&QKPiARU7RMvhCexEIXYc zW!%`!aavdzDk;}PJ1txADadlfr+9f|`Leyrv$jNo43?E>|3;t)Vr0I~?kk)!Qj&oR7;{5Kt4+M$x55D~Y|4%5i`q0=|yMSK< zUpCn0%%+3_SPPUqq(yoMBM4h2(GO)X80!VY8X}At=TNpGW~3TTeaZhawRahH4)JCi z=aL8kz|4x0UThoh>JttGM+(9^#Jo~qNF3K3A8BD;lr#Yw>qLf|N{2K9hl5L_gDZe+ zhMch9n;9^*A8$ryl3pG>+&C{_>KRrcOnjpWjzK%d$>F&1a7txsxIkK^IBE=k2VZyC zx?DMlz5v3F<8&zv9-jc?>uvU$hH!Q6FfOA~D1s@sk7FoqsXls(dn>~O`@^_)!+>nS zOgo{~%R#Mb_C^2F04)Cg=0MnQ1Pyl33>a-RxU>e=|6mg?rZ>nwOI%KXO(Gg`J+8fN zWh6(}+>B37G8w-1hx6|-qu7*!vkG(|`QNU<$O++{D>Q}wJz6*4LEh-v4qVh&L`FXd z)cSv}LcgTgUnr(!?*Fy|(_b?p!pI4s=C{IL9mCJqC<1Am4?m)kTJEc{+3!={vOZ-@ z$xP@`0t<-t(v_DKyjYqK0iymF^ol5zAa|NW(8;ia)FGoae8Y#0!}Wd@$+Np!Y^R{l zRnpE_y@3#2E}pJn9t*PrA*u4KOyo+x0_Z zB997J`#U0)$ONpJiP*a$mAG6^SMa{K51vpZH59?sbb_=jkzG1t=~O$7Dm>4zOAeps zERChJm&FZrj%Rt|g^7a2k*j&V6FnC1Bi=vJYw-c%gA;ugUqgIoqTk|ci4RWWzuP44?VvWT&65lj2Wbw_!w@j?H_*UXY+*qp(w-LX6Vx6~*_#MRW4$UkEq zYfhf6pcRueRxw~3XQ2>0jY^2f{P6+94`d4KPfde`SI&U(hY25`Uk`F&sy@FE^n`J| zx7O!p8|iKLAA0+-y$78hb8gjwgHU2UG=n(Q~|hXs>=yw5fQc_8fcr+YjB(Q3xND z@D1>JiPr*e=zJ?|~TQJbqbdfT}>4`D(W^xC7-MD2NyMvWf1 z2(-%V{M@|dUHbzTRd{wSRPpBFEliJ(js%&K=^vF!L9V2O_CmEnyih91b$~RJN+n+$ zBtOaO4f;x@4hOV6>er(bkNmWPi~`IXBBn&+!#F_-uAoOjuYx`W{R##Y3@VWNh`(0B zItA+$Y*WoEflfC+&V%n{o7BVw%tZ>ZemG#ZUqQ%L{aQ|c%--=qo4e62Dq0a_scx7< zLV`e7q&I>JjwN7L|k*WH_uyi(3lT#Ck zIYHAeDx

!7$kUB-Ci=n+fwDQ)>SSbu2-7F-QE{X(1@Rytv3bZIy#bfZqypTv)uC zFFWg6R$My?m_5il@l3=@mayLGp3(_4rB3oj7Z%E=tv*w+93$~T-a0kNS=Tm$U^t8c z0O1A!iEum~s$GiJp78MjAGsbaUXMO9?n~jNt$caT|0(KQH3p`9`LpV<@D3D&P#eWP z+3hP5S4O2$?cf%-G%jJ|Kdt(nP+-ITqGCTzfR5KTzP##FLEoxMe_Edv+AY)Z93s%n z4gYT5yY>)y*i=b>uZvr+!D*vMMH1m!fioiY;TehQ#732v2{hY=IEdzsF$AO38}y$> zJ?|V?8=}RFCX5Sn=I2Uh7Z=XL&tU-71Shw~tKg;jMdU+4LLt>yx2&$DjmzpTu2Z=M zDAHA3{W*12`Zi`hz)9qO6JJ3?bGN>iAm~~MT3A|xhmB_t_SF^^YIv97ywnj_x&*xn5j*0RN&uSx zC!I?nCC@BK^zPSC2Y^3_+AnyG?XTqcN&f;%)v^AJNs}}23k*IKM4W=aMHT2) z48RFFw3Dd^m@p7v@*ur+`|Fa*KT3CQsUolm&2?s>gX>x=onCa~Tzf-`Zxfy1LlNAF0fgpS=JpYp&MF+7qE<*91!&z2!?>EKtWb1?RAhN|-iiW$C_4 zATj60Ix-elXH+Wyd2dV=+bdbCQqts;jhz0CnNRvguf@#o@)3x|beKWsM zYv;gwj$>@z;^G0^EUe&* zT2|u}W@n0!K!6f^D-n;qiR#W*63N)#33Vo?` zwjx~idZA_RE0x_K&PwqIXjvD-S%;)BEnC|nLnYRRxhzFXNEieH9Fee)k+?JyDB_xO zWgXVDV67Oy(8^J>#x;jbGjNYkHZ);DN#r0Ma z$yv_3gn#a7cWt-9T=!Fn5R0_IW9*&5&ifMRsgu`avU)OaPSJU zyH7+@54@F$^B$~+~e#rK4p5|b9sOe*y2FMg$jMwuBX`etT1jgq%GQ5mGhX8mov z85t1fSRZ**J@XPXc)hk7bF31K7U9-V8dJnRo+4<_vpR~V%RTixUx%9m4~KWJiZj$R zm$KY7#(VAYhMFn3wBP5gb^q(V_H;;3Ew)xMTHRak`zlZRtZg{zo{GWA>Y0kw`)3A# z+rVrf_6Bd`73_y_jMndMS{j(sc%V#5YRl39Z(^Np{H7k9hHurI_Vxy@h#6@P=su+@ zd3sZr!G7J_hB|KEy8|7)!R8?SgI^^c^3q^^ZF3Nj$(`(~Yt)(J*|5uVwc|*9kb2j7 zcX9s?_c76`d(h5jD=}`?8RKu4L_WA6-{9@m{wB83+v8$0o47gWKJR|n zH1O_b?!4|fTBM6(m)4Ph zb_~9JrgT7Y2Z|jJjw$zpyyn#4eo%;n0DZ{M{_=swhU14Hf8uy?1*s~w*;L4kkm7LH zN|0?IO8ze1z9wxQ|79wpNVMjkAumWZoIV#Ga>M(8JoLBLN-R@at=ku@=`7S@j1R5O z6p#FRw0TWZh>-FAg#(TCPllKBn3U?l2kw-g2Z^bwCT&7$^AM!bmIIMh)Pjje~;R|g$@WOsPeQEL!D^BV(4#QH3{ zdx=;ejmjJE8~U-w{;=^=2mE*P5ts@;qdr`8AA}`(@{RM4Dxo2Qh>rAJlkh?nZBS zDpEo}{C`WQFym!O$^O{zWdFGOn32-}9WqW`Pl|M0IztA>?WBY?zSD(q;u$AG?@)bt zRF#Xpgb}mMpAd%nA5$He48!=<2QTZ>33*gD78EXKQF%4Y)rewDs!3n6F>&!l7f(Ab zA}*fXa&~?!T#pZe{W?z9Yg#31{O8mTBQBm-V-kWHLB_B&^sW-74%{a#ikf1nmHsLCo=_=oarFBgecdWMA+Vzo}(4$*fz6n8UYis-vV~s=K zgl8;=%aM=&2aaxc$&t8BLUy2~rEV{$emFxztc&6j2jo>LlPB=PpnF@On6HwD64b-I zk>XDxT@Yf>mn=w#Q)LFh=Kht3ytHX-6u4`o$IHm`X){mVW>EubfC`<|=kUB5!*9b~ z^}ZQyGBMo_gvoNR}vB-Maa&)+n?dVT3amBRo~U|3HY= z{Oke$L)6%kf%uD!8NM zyY+g!hG({n#QXw={wpf?ssgF9``@b~>t!1eZf$pfVmWL_r!(D`*Czboz*^vpc#?m6 zkNE#gZ+{~J{x*I@@gG%SB%bCqC>XvGDs*w*7~f{VO(nn>eg2=RnZ^LNHS$aJM3>R& z_|{2epAXXEM6Kp3Y%;*U0;xLhiLCv5LG=dpB@(!|7WI3PSgy#=yA4NW zXWPAxJHvY)TUV(Rr?{xIW-cpqk>u`weW)Qj@=lJ`+rEAyI67(EoeuuL(nUrrJTO7S4`7WTAv8> zLak3&bflc34oo^+;Yw8}C6^pxBfpFog!PUp3^%+XJb5|xa^z*KT3g-2(!l(ptj-N`lcUf-33ly^)) zA6E4etZkSsvOUh#U3vh9PahZq3OT&hM|%fh*SOjS*>P*wacksyko`314cd-dW0u83 zQKmPYuManSU{garxE(df)AR;NK1+{L`y?TbDGnt66X6RzZd z3T*$D_{**dh3U@G3(nvgqBe$C?6R}$(%}$dTx2!7JB8L%#iBqD&K(z= z1@r#YmH+e8%U}G=PwY2*OFDP88ieEg84m6}kj3)B1-fG(k;P7wW>jnO6|O<|%4 zzi~eteVKV~#?}GdUxpza$~l;33#yrj=>-&JIxfn#7tj)TSa7%;{UpW-%;?J3jVW#-&OGcDUgyTPIqvy3I6<14Z`|Wu;1w#_n4a# z%}HtU%oIEn-F|4%e8;k*CY=Sl#VTw8`Q7phymAVZcCQS9Ee5PE+v5QC39_rD)G-qS zHx0Z0-&I0G@Bbb_5Pc_?Ow~&>HcjqOc|#Ai;DQZbeR=KCI+WzLP(ZiMFh6iI{QnLk zVRCp{?#=!bnl-UR!dvNvMg~nd9Y8}IN5bb?hfoUIeVAYQ%y9zd@J&*w zzxLmvdDrvOX9;QGP=c^x037Cb8%zqJlgQ<;jxkA=@aBfn#X@r0d47+WmH9rltc*_p zs!oGdvIMCqfYL>gG9hm=DMAdt!*^p!iAtRaKW?u{p*5YtZ6Tz!oLeCOM}}k^@w6kC zu@sK^2Rq+U+_}Cuj}!_Jn=GGVO&L6MDALW`zGBEJ;F399^Z*bI=!yYmE)(1p;fPr} z;Oh$R-pOSsq7=c+>G(Cq?W*%UhNsrk4w|x{w>z_!?M9Z>#*7cKeROnmtk}ApvfUE- zMmCI&=9FvZk>_FSgg1Xtg3wBNw)19lT|Ey^$laOAqL9`yDlfG1Ps(&Ll*X|k^wT5qH;r`5mALW;OnrqC^V4yj; zFLM4}=pp41`jBRb7L{L?&@Jk6n?q`SvsFH^+=m!UwnXuzfkqC`FA<(0IQII% zGuMU+JOhgu4C@!SkB_4ily9!p9mkB@NfR6WFa3+_{|!wvyHkJD4`FGg`4?^EHZ4so-v8=nHBHH0r2 zm+G+;eYCM-oM+ZH8Fg>Xr=%D%gorb~w0>#B(nc_Q+Iuy+ybg?UeMfvlM|@*P9G+Uc zeD6veY&x;L77m)Tt7@f&HL~qzW@Qzodb>*YdYvVEojn%d*?Md^u&xzD)t{Ljd9_b? zYI%JNh#s9q#{k9GYEjas78bbEaHO^{D&$O=tf~b}4ARIKseByQ7eiorv^;s1S^zi= z2fSJk5OT6O7lI$z^%pBDX}INzu@)BqqZ9*jdr%{12fCa!*xHr+4lA$e2X_NvI-X&R)w;g(+FqcAapN40D5Noj zum4I<%!FBTuh+TZ26)ne-kP~;8*Ky;{zHLHF|PeTB9L_t#CdIXcApbx!)&FOQbqC5i(i>z6wPUI%Jxmo@6l_$1VOst4jRY6$htvCn(TOw zcItgbeM<`WJsrdkDmYX+?VMwDPvGjKaob9ZZo)#XI>z)S;Ec8o-I|fOVxG@AP zQbHX44#s4`tdmm)OFK02 zf(LGF~}kfFHl@a?6sA`pv$F<_tB5^X8>P{@KU@C+>|^@GP2L^ z_i%nQqw5l|8D{Wabiq&{(PEf&MWp!fWL z4f|naJrMKW$P!om*cjb78}ZhEUvz0b=iU%XBiZF)5doalH*CDJ-UgIk+?Z_$b2hef zHikKy+Buu*>zW%lrPeh!UUm*EV*dXkwdt~4l-`ijin-0MjJFwVg|94+@VCCX;i`R= zK5S@iXl`t7vR-WQwxW`Ao8)BP?X-FX<&<&pn=eE^5%tI0Z$qqHZRKoit|fm{b8X0} z;!2g}&EB@=W{fyDUrj7+A+5i3bBpE|jB6_?jpH57twZoF!tpG2=7OQPQNGPuIMggQ zZwp&^XRC!g*@~)9(YxCubou3?x83lU{PJzyXc)WQyT`?5wn_cR+v&>A+<_9%?cOeL zx3|aJ>(s1Fo|Z_UVtzRug_B~cIWKle)a0A}zQ+Hw6Y)6G5F;0W8R8#7d(0TMv$N*z zd}(0?MmOkq?qFbr=l*-&Yuo$1@728)=*h@;0v$4~3vc=bZus_!rx%NIc4_JtW5s9X zk_d{18;mY;SrEr_$fn|?M+wzhb`l!KVIu7B_8~1@kFSdc*!!y_hT)m3KwW^D^2 z6FYTq8uxycLbsxw{1oAhC#+gukc)66pC21@BXE3%!~*oSMmB zwiHi3ci`F%)hYxbRGW%C+SvRIIy9g!+;z7?;DM7pWFzCBkBvr`0^^{+NW;|FB{oQSh>Y3ITfS$21%!9sO42{v!n+ z)1y0A$wgPuz`C4j4&DA)!(O;7(S$1gKh^uDJ2Rk|a3Ii@@*IWTf7w7!Rpw{S|Cqm? z$#Lk+CskZVTvD)c9`MW?Qjidf!s*(?L23=#32+)*1t+VZBgK{HZ8lwys$d9Cn{L;4 zo&Jqu)V`wXyT!03Qbm_&yV_uqno-5>QE;yU@hp8qk9I0KreK$X6{)M-51TSusf?HI zR{XR%dN7RC-^~jM>4mX78L!X;jfVc=PusAu+f@j>LDz%JEd{z=VVg zNK$nE!n6vfF4;5e2+dEga0YmSaxCKVxo#R3BNp$fPp0l+Oor~^JNBjU9mQAUVWzGN z^Afq3SeO{QL~5U5)8y8)q`mz~Ax1)yoQzjoKs1Dfg0+#Y)tWXVz0NyH0hzEXrP`!- zBTX9T<=DKwMejf!j2&&;Rf7Afq{e96#JM0vHQbO-An_81X@lbyw*IHS+?x|xw_0aM zT10n!hkrujzqUp9F-V*Tpy(LPo$2y%vz0E(i1Dy+t|kn4DIjAihg=xSGYaKxuEv|l zzgGOn;j<L z#i$t+#l=M-_wjnp!=a#;2a_qhhAJ_*R0^jG#w|g@oRVofXSn?2%0X4=dOZAaTT6b7`rwbWThcqH(r6apbGOGZ58hbn#qX{w z{~5|7rD~V&v+{$j^57(JZQvK>ka}rc@Yjnxs|WkKVYoE?;G^r%MH=A$h^Z^5;Ewdm zH?Q?i(muTK+M79VBe?1&utazj)Rgl!BLnPAZ*kAPtuEgyhJI+);>semJghljb=(%> zwL_M>F3i0hyml?6M!-(jdfOPKVbq_I5`h7)t8<-~wR*U@4jej#qIG&{q`nckruqU$ z3g-ZFM9_!a@6Ic^rS+X(T-wn2C3ra+`+~4ZS~3WZyeZ7t*v{FE0TBw>81}%Kl3bD4 z2r*Q-iESs=?~RI!Md}_)-5aLXyU``f{Vm|bTS-5nFIsGWJ9fZg37$UMjvZ{r z4z**4+p#0<*ips2W3R&Vgp0SOxn52`!>df(tZ`Kr@K3rmka zRhDU*p@EZ;HAxlV?bV^g&tpk9;hkST5(?z8 zlkaDZY@O7qt|f<{oN?_WRk(RHy66G5+986m@o~Q@G~p>u1z#ps4q=ccN@LG0wdq&< zYIEhwPNz!jxQ4Yc%f0*HGviN;KlTt;QqFQ=-mR4K!&gjpv+t4O?lH{(RKrCO zBS;QoAY#cYQA!7w{}Fd{OpmGLgsjCe|D4dS&f{;zA&!%gTu>Kh72MBB>CY)9WWBLJ zw{z#t;&JEqNHQ5i5pF&Cor?<`qR5+mUe&dke}-{2_NvOSiTfYCr-Li1dK!s?pvL&Q zdK-SO4iRl++pDUP?dH+W-1oGJNs;U*n1>&cpSB2YYtjYixhU zsXlZprp~?L)DuK*1nW$Aj{ROs2wZXkF#pdrRqITqev7Anp^0gW^qqkBFO~OK3O=Xc z^9m#%_m8Uc#|Vy{^#2bTSk14{Jz_>{F=0~BMK9jgQ3hUR2V=dpryvG<9VW{i>kifJ zN<#n)av#A0_m6>SwaeU;N4YhS!8Pg-Gh9~b6y6i8wM_g`hnZ-4?k+uOJqU(v<2rwM z^Z-R_E-}M@biXD+SaU09+#go&*D8>6GQUr;eg$&W<%^%@uTd~WK$kYIaN|a!u4$P@6wWwhulH!iCsCJHq^wW(+vSDrZK{)ZGZruGrV%((=GE;g=X;Tku%B;mG( zQF}v=T4~>C&Nh81YufrM5;Vx!E^@{jiend`x#?i!iz<~lSsmhRYGq^l6Y9cqdYIQk z=lScVdr)rIGi?slik&VL?R@QA#lvU7^rQ)ts~-}>UI znL6OysDSxKKZYGfe%;`%=(k4V7=x0QkN1S}<_F^|&l7JAVQkzgtrYLLw8+4(WN?Y3 zt8T>E-@w&iA{gCsnR(cIA*m-~kPo#}E3wqy_4YU?n5Lp)t>R|QHw#!l7qEWLVEvrI z`niDha{=q;0@lw3te*>5KNqlmE_7HwXRv-w4a&EKBn^Di!pxuKHbZL4qy)!J3f(PM zj^sp%`3v-ZYvKzrb6<;Vr&xk(5;8shC7T<~Xb4p%Sw8I_X?V4@RfJrK+-!OKH_v@> zxxh5jh9Ux1^iU9q3SE?oR7^OQv44ilBF{YAqlWAy@#BI;tX&QZG> zOVolB5WD{KDzZ(%DFr{-PFmg_N>VQJPP7Qz)s?$$);f?*16-o|Eidjvlj6^{dkrf5 zg<_XgsR0Djzc)x>TQlK&ge@|5PO0K&RNNZUZcVW9YkdDOtu)UqIOvLPLYtyD<~OJf z_R*yUvrcxL4r{D<0G*Q}`d=37hs6W9pD#NVB0N{X2au~p@OWJ2$mJr`CdRO`tbv-w z6nP&LD4(ckDKwmlky9^{uT;4lHtlfAVLpsa33MqcMAWw6W=ax@e+XnbTN182T}+{|P# z>TG)$agt;6Nhi%L!%BT&o-5X8-~gc)$F@P5%aTLKZY$N$hvK`KRzuYt`YOy7H`Y@B zarJ+Tp1)tg2NYQG6>L|SM6KaD)Nc0hxBez zc+R4^p>Y~U50krMcj#Q@nTC;$KU(%awW6tDn-;a{gr~JCrBJl2dgZmn-^W9Y{j{0% z1hd(88y!6&l->q!^}ZN(cP)emnEVo(W#?1c5;N3~%^Vte(#n$;x_-dBl#T@EaxpIC z6MyeUfver&*BO{bXAh9q83LfqdmXxcig9<$O3Z=G#Aj1Ca<}^*p_^S>i$g|=vuep=okw0UqAZ6R9()*2TbkfZH_HXkY2Z5=5t zo9akqNM%{FIW2YY5`GY z`-i@)@Mi}}nd1;aT~s;cOGwK@l3BSAk-Z~v4oGYsaf)*#W4stw78^)8aImu2BzBa| zN_yQ@7Rj}^sFcmhZBAGg4ytStuTVc69#^YedAgEtQNL!%p>|2^7refXI{Un}a^-Nv zagJQQBkj^^0oR>c&uBrtK7cWscNyEd#qN(G$8I@#k$3z^ypWt`8-)9mRul zi>FS>JgVya4%NC27is!yG!F}JTP9+#0XH3>$HK^N<4}>CYYq`%V*D3I#^nyWSPR&Z zlEcE+K&5rKTo$6u1rIJR%@;tv9LgZ;arS&jT{6qbIa=y8)Ggd+f~JU(CP50cMUx_%0!7d^ zX@dUIqOI#nmek*O?(EKxlAX39bob7k`@VN(&OPTl=VYnl`mzE8P_}vhv9YmZM{j(; zlDcQ)@WEJvA~P!4{euTSe$U9U2Zj%#+eeRW@sxm8`M4T)e?P zRm&fe4)YVL0%2Bt&)QHWDD0vtzD@`yMJ92&CC1n3T1%q|9^Lm%&qx^KNx7VFO#Zvf z0H;<4U#4SA_Nb0qJbiADet`(z4O9~Nek7(7>`&AZND=Oeuh(??qT4IDvG%Jmv_$y| zC5u^PJSa9y$wYX)Xd^xG`*D3_VN6 zunx)iR*k{mrf@3h_g;ubtr7M`BP`d0-7A|u?ajsY{w;hBESly1POgxbhF2fqwC>Vx199(Sw=)*ia!0imr86F--MJ*q#)Jtim;a}h+;jQw zT9W4S@3NXcR=dM&(XYGd7m`z(yFJ8_m@57%MpgZNTq7+F`?sH~1J$dk73k}_`8;@o z`E2<=%qOe)at0HtS&UTHzvqRFf3KKi3#~6E%dPBy zO~_#rFD0EE)<5bG`}e-m@KWN_3GyFP{&K7GPbGvcSX}*La@t$$x@hhB_u1;Z@P3Pr z#6s0mXa4s_qk7P1*^FGe)gJ>}@BoMf%&RSS7q|Ei`VaXJ7kd`B@@#Dr#o#~UQl@)} zf5PHT(|yDr^&jKvxc|v8_svFcobZ3zKVIya?hR9K5fYT#Wi{uO|%1cmREyu=wP%T@n>Wexwd|1$~Qrgoh`2= zoPD){lILpA9cLcw#@rNN^_P?9cNK3Fa}*@RIvIXMzW>Kz>*sIK?z#_cH-HlS2`$1F zaIO3>Y~}jByc3T`A2zze2ei^OeLLu!JBoMs#S46+i+B2SqEADF%ANj`RjG!tS*jvU(tu_9R*({;21f3F8nU%oDs0_yZ!t65~s)K5TN2G3t*styTQ+@`u@Q= zjD7 zh5-v4v;c$;xjeod=(h*GG=6r zPiUN(QpXrX^aQg4mw7vLV@&W;+jYCGl#O6Qr5srek*Wo=D`(J&tLvsCO^I>TbOWjwi;8 zBQjO^TOyoqiAwtD;e!W{sQL!SL-?_3RArS0 z;|l~|RDrh?{GNih3Brk8zb|}ytl1ZQQpF@)soHvsKM^1U4(fi*8gVHLdt=DVMSG?7 zyTs=3j$EgL3A@h=9@OjWqBTzyrH{GH_XV?$hqljnfx-kj{rH?%jn3rI4Ut@hfpLMe zT=^jZn!Wf}qV|r}YYn;hfHyjS2J@T+uW$-24tsWDSh&#KNgN_mFeH|T9fs3Qemipk<6ZenqpST0!9Zu+_|r6OrerQna0Zy7&n1$Tut;q~d<%80k*+|RtB z*q2qsHpNV?XS^!2cfCewUn9U<&)n%5-yPP4WI-R(L@9SEhkzx~5#E7KM5Dp#b{}U6 zUmXAHTXatbh9Mg&O~KsR6K0FS@L+*pU%+5MaadxdrH0}GBQ{g0WW$dI!an<`SDLHU z+l@*~6;cw~Yj3SxJEQFtjF-r+|B16kid33wRy|2~Zs(j%#BJ(GTg^I8=%2=qyXg`@ zrJEu0GWecmyrST{3XC<9i+`&g(#PyOOd7}x@jtF+;D5-Nhz!Na!bsY9BjS%3zl}WU z|4IStkdQKE4z%E#poq|s&evDGtXnx?7mcWhRF%>3Q8j)si-V|~hCn3USC zw&a|9FXM|}?#NINsFCpRa(iQyUnm3&cKPY0a#9vuRXB`G_=!n3H#Vx|`$;oEs@(j2 zUt`K{#o3soZE%-5lB>k01cg$aALV2AaO48;hI_Oq?*Laor!?jf1It8wjl zuywH9I6gzewz>}STFH~fx4^Z`icsR0NMC75<8b^vmi;7qq4uwrSj~3nex7tJ(}AJD zh_)rnD<{%SdFp(>6WLF$9sNnBhi48LHtLB`tC&dvlpux2cp$`DKjzF2GU zC_F9fWV!Q}tcRYitCxpkOv%Gj%1pNw7LW(D8TLc&kYt;*d9nQ>uEdLS^M@`hTn6lN z$;=FPDW>7dx!gBf?iS0{63?~J3hsXe)Z**lF2tkY;?_8pX471p(e(=P&ePUSMq-^h%EbSZbhNx+7h z>I%}7=~VyD|2=vS=yH&66`B$2v5(PN+<@qR6)s2+z4MD{d!bP7)_dR^NV?<;Sx-%9=M6N&fCM^e*Koa zY|Y-m`p3YpB^EYg&!nqN@6ozNLI*_mbVXrM#uw{-iXG zf8h6&yPvF6tlTYbD0(XR<$F1C{+q=f^i{f8-=1<0?IsHwid)M)$QRnN0GM3Zy0ERh z6~0U(qu6Rss_A|Bw_$6K(b#I=9Pq!}Ln()Wpm= z$BpP(&ZO}CSs4(4M=^oTjN&fOcn`s|yjP7g5tNCQwhYh47zj|nf&sxWnZs+2cj##> zMxLD-a3HB)QvMBw@*ngi$%dSbF%JxUpB_?4uhh*Ux=3viZ-R@OL0=DXOJM*v)pO_i zi?jVSDeg$kKQgy^JbdYRQ2QL)nTV(w$lO^Gwn9`7l-4_)19vH%gIW^>NYs>dX^V#e zDi7Web~DKl^Oq1CYVe$<d?Krl-|y?_`5G9`bu6lyZb}` zG08)RhX;>%2RO{_9~m4wd~Dbo9UB}wWzo0q2lr^X^QSDE5~Mex%@E73@&(nu6a}5aIC*=ww)Op`9v?Ae&@9 z?t^W9J+Ja_YQg6fJgdOG)qYAh%nNsdn#>KdJ75Sy8tmq-*qK0GB{?*do8f$o++oaV z8w=)QZ8SEJpHq#WQ^9sli~+OYfU+mFTOB5_Nb@@JFvUWsw$ecOL;!u4bGJ#jUZ-5K z0@6&&7&0XyLu4)iA%s>|3KsdJvy zN4K#f9eSDgr?uVjQyho05>K%hm=4_9-IQd$;; zELv7{vfz&S#72IxkPF>Tu50;vytYF}%Z3b8E>v*_Wn-O2B3eE3el{%}02DIv{;CUK zYt4X=xbii8y2*0!^e#%;{betMbEj=#$E%oe7voIuZWg<5GW4mF3;nx>>jE$anY8TB zsFeGIvJd;m=$_3AK^OjIuB44}=VWqq)iTU{=2GfEQY(U$(Ask0s@Xen0{)QraKRb{ zVo=&c!ZCmA^7ySJz|0qB!^E_p3Sr50=VClxa-lecQI-&JfK0pp5#m=6BDn(i0`)>3po75s(n z+X#(;5RM>9wY=X}y(#zcwY<+(YXwRf<$YhNR}@$&^-y!K$H+=Wsqs=fbl*l}_m%U* zY7KkXr)HvBUs_)3sA^T-KzbfdCuY4erGhafJfOfDv`I@jG_n7uLRKol&Vz}{EG5gmv2-m9m~0_%de=``2E*alTA)7C-?Pgw=T!_Zgf9; z!~TCoWmLcW5=PVXXQZ)YdI1NDyYcBhse3d&0gRR`WVhbfQ4O|aiJ+HMIv3sSQ?_Wz zEZcRvAf-|%2U3q3vgXsyjg75cRBN<$S;s>V%FzxDrfd5CCl)gcDBIOQRBgX5EzJy~Sj?=VM8nV$GRdYuXTyx9p*CPP z$-HWcHdOeumjn0Q5X=1g7jMMSTS4NhQ2EmE;)B9EvF8 z9#z{GB@CRKm+BsrrTsepR}C<0abYV3Dv_MXTMI@FaKeHwwP1z>XxB731A6cj8wdVE%Czpmg71$H|5 zQ^o8!^52TtnWSDTX0jENNtlS@R;BG!V7_baR%}p#2`fx&VL}awEzALgQ=(03uWwT@ zsz4fxrYz?UHd4WGG98n%7$#9{U85+2wW_LH!;gY3Tvjr7W2bm-WVf21)D8b>y=7=G+R%{)SvLX)VcI{t8|&68S4!=W-h}o#gGytk3L*Ex!iSj}7@I z($l$4T#T*9y3zZ|uKbm{XOgS(S9AB}yxjdcFF7S*JO3`J8Qsh8rPV{3TXJKW-rUBz SXX>8G+*bEYb|=uyYySfe^%~Cr literal 0 HcmV?d00001 diff --git a/bin/calibrateHeatSource b/bin/calibrateHeatSource new file mode 100755 index 00000000..17a61437 --- /dev/null +++ b/bin/calibrateHeatSource @@ -0,0 +1,2617 @@ +#!/usr/bin/env python3 +""" +Calibrate the projected depth distribution used by AdditiveFOAM heat sources. + +For each experimental condition, the routine builds a response curve from CFD +results obtained with A = 0 and trial values of B. It infers a local n from the +measured melt-pool depths, then fits the production relation + + n = clip(A * log2(depth / spot_size) + B, 0, 9) + +Local posterior uncertainty is used in both the global fit and its uncertainty +band. Simulation results and calibration state are cached by process condition. +""" + +from __future__ import annotations + +import argparse +import ast +import hashlib +import math +import os +import re +import shutil +import subprocess +import textwrap +from dataclasses import dataclass +from datetime import datetime +from io import BytesIO +from pathlib import Path +from typing import Any + +import arviz as az +import matplotlib as mpl +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +import pymc as pm +import pytensor.tensor as pt +import seaborn as sns +import yaml + +from scipy.interpolate import PchipInterpolator +from scipy.optimize import least_squares + +from reportlab.lib import colors +from reportlab.lib.enums import TA_CENTER, TA_LEFT +from reportlab.lib.pagesizes import letter +from reportlab.lib.styles import ParagraphStyle, getSampleStyleSheet +from reportlab.lib.units import inch +from reportlab.platypus import ( + Image, + PageBreak, + Paragraph, + SimpleDocTemplate, + Spacer, + Table, + TableStyle, +) + + +# Input data, configuration, and persistent simulation cache +def _create_data_fingerprint(data_list: list) -> str: + sorted_data_str = str(sorted(data_list)) + return hashlib.sha256(sorted_data_str.encode()).hexdigest() + + +def load_yaml_file(filepath: str): + if not os.path.exists(filepath): + return [] + with open(filepath, "r") as f: + return yaml.safe_load(f) or [] + + +def save_yaml_file(data, filepath: str | Path): + Path(filepath).parent.mkdir(parents=True, exist_ok=True) + with open(filepath, "w") as f: + yaml.dump(data, f, default_flow_style=False, sort_keys=False, indent=2) + + +def _resolve_config_path(value: str, config_dir: Path, key: str) -> str: + expanded = os.path.expanduser(os.path.expandvars(str(value))) + if "$" in expanded: + raise ValueError(f"Unresolved environment variable in {key}: {value}") + path = Path(expanded) + if not path.is_absolute(): + path = config_dir / path + return str(path.resolve()) + + +def load_config(filepath: str) -> dict: + with open(filepath, "r") as f: + raw = yaml.safe_load(f) or {} + + required_sections = ("paths", "case", "calibration", "bayesian", "report") + missing_sections = [ + key for key in required_sections if not isinstance(raw.get(key), dict) + ] + if missing_sections: + raise ValueError( + "Configuration requires mapping sections: " + + ", ".join(missing_sections) + ) + + config_dir = Path(filepath).resolve().parent + paths = raw["paths"] + for key in ("experiments", "template", "campaign"): + if key not in paths: + raise ValueError( + f"Missing required configuration entry: paths.{key}" + ) + + experiments_file = _resolve_config_path( + paths["experiments"], config_dir, "paths.experiments" + ) + template_case = _resolve_config_path( + paths["template"], config_dir, "paths.template" + ) + campaign_dir = Path( + _resolve_config_path(paths["campaign"], config_dir, "paths.campaign") + ) + + case = raw["case"] + calibration = raw["calibration"] + bounds = calibration.get("bounds", [0.0, 9.0]) + if not isinstance(bounds, list) or len(bounds) != 2: + raise ValueError("calibration.bounds must contain [minimum, maximum]") + if float(bounds[0]) >= float(bounds[1]): + raise ValueError("calibration.bounds minimum must be less than maximum") + + initial_values = [ + float(value) + for value in calibration.get("initial_values", [0.0, 4.5, 9.0]) + ] + if not initial_values: + raise ValueError("calibration.initial_values must not be empty") + if len(set(initial_values)) != len(initial_values): + raise ValueError( + "calibration.initial_values must not contain duplicates" + ) + if any( + value < float(bounds[0]) or value > float(bounds[1]) + for value in initial_values + ): + raise ValueError( + "calibration.initial_values must lie within calibration.bounds" + ) + + max_simulations = int(calibration.get("max_simulations_per_experiment", 7)) + if max_simulations < len(initial_values): + raise ValueError( + "calibration.max_simulations_per_experiment cannot be smaller than " + "the number of initial_values" + ) + + token = str(calibration.get("token", "B")) + if not re.fullmatch(r"[A-Za-z][A-Za-z0-9_]*", token): + raise ValueError("calibration.token must be an alphanumeric identifier") + + adaptive_learning = { + "enabled": True, + "initial_B_values": initial_values, + "B_min": float(bounds[0]), + "B_max": float(bounds[1]), + "max_simulations_per_experiment": max_simulations, + "depth_tolerance_microns": float( + calibration.get("depth_tolerance_microns", 10.0) + ), + "posterior_std_tolerance": float( + calibration.get("posterior_std_tolerance", 0.25) + ), + "pchip_grid_points": int(calibration.get("pchip_grid_points", 1000)), + } + + report = dict(raw["report"]) + report["output_pdf"] = str( + campaign_dir / "reports" / "calibration_report.pdf" + ) + + return { + "experiments_file": experiments_file, + "template_case": template_case, + "campaign_dir": str(campaign_dir), + "simulations_file": str(campaign_dir / "simulations.yml"), + "calibration_state_file": str(campaign_dir / "calibration_state.yml"), + "run_command": str(case.get("command", "./Allrun")), + "keep_successful_cases": bool(case.get("keep_successful", False)), + "melt_pool_dimensions_isovalue": case.get( + "melt_pool_isovalue", "liquidus" + ), + "calibration_token": token, + "adaptive_learning": adaptive_learning, + "bayesian": dict(raw["bayesian"]), + "final_fit": dict(raw.get("final_fit", {})), + "report": report, + } + + +def parse_experiments(raw_experiments: list) -> pd.DataFrame: + if not raw_experiments: + return pd.DataFrame( + columns=[ + "parameters", + "depths_list", + "normalized_depths_list", + "fingerprint", + ] + ) + + records = [] + for exp in raw_experiments: + params = exp["parameters"] + spot_size = params.get("Spot_Size_microns") + depths = exp["Measured_Depth_microns"] + # The projected-source closure uses depth relative to lateral spot size. + normalized_depths = [float(d) / float(spot_size) for d in depths] + + records.append( + { + "parameters": params, + "depths_list": depths, + "normalized_depths_list": normalized_depths, + "fingerprint": _create_data_fingerprint(depths), + } + ) + return pd.DataFrame(records) + + +def parse_simulations(raw_simulations: list) -> pd.DataFrame: + if not raw_simulations: + return pd.DataFrame() + + records = [] + for sim_run in raw_simulations: + params = sim_run["parameters"] + n_values = sim_run["n"] + depths = sim_run["Simulated_Depth_microns"] + spot_size = params.get("Spot_Size_microns") + + for n, depth in zip(n_values, depths): + records.append( + { + **params, + "n": float(n), + "Simulated_Depth_microns": float(depth), + "Normalized_Simulated_Depth": float(depth) + / float(spot_size), + } + ) + return pd.DataFrame(records) + + +def params_equal(a: dict, b: dict) -> bool: + return str(sorted(a.items())) == str(sorted(b.items())) + + +def append_simulation_result( + raw_simulations: list, + params: dict, + B: float, + depth_microns: float, + simulations_path: str, +): + target = None + for sim_run in raw_simulations: + if params_equal(sim_run["parameters"], params): + target = sim_run + break + + if target is None: + target = { + "parameters": dict(params), + "n": [], + "Simulated_Depth_microns": [], + } + raw_simulations.append(target) + + existing = [float(v) for v in target["n"]] + match_index = None + for i, value in enumerate(existing): + if abs(value - B) < 1e-8: + match_index = i + break + + if match_index is None: + target["n"].append(float(B)) + target["Simulated_Depth_microns"].append(float(depth_microns)) + else: + target["Simulated_Depth_microns"][match_index] = float(depth_microns) + + pairs = sorted( + zip(target["n"], target["Simulated_Depth_microns"]), + key=lambda x: float(x[0]), + ) + target["n"] = [float(x[0]) for x in pairs] + target["Simulated_Depth_microns"] = [float(x[1]) for x in pairs] + save_yaml_file(raw_simulations, simulations_path) + + +# AdditiveFOAM template rendering and case execution +def format_value(value) -> str: + if isinstance(value, float): + return f"{value:.10g}" + return str(value) + + +def render_text(text: str, values: dict) -> str: + for key, value in values.items(): + text = text.replace(f"<<{key}>>", format_value(value)) + return text + + +def render_case(case_dir: Path, values: dict): + for path in case_dir.rglob("*"): + if path.is_file(): + try: + text = path.read_text() + except UnicodeDecodeError: + continue + path.write_text(render_text(text, values)) + + +def find_template_placeholders(case_dir: Path) -> set[str]: + placeholders = set() + for path in case_dir.rglob("*"): + if not path.is_file(): + continue + try: + text = path.read_text() + except UnicodeDecodeError: + continue + placeholders.update(re.findall(r"<<([A-Za-z][A-Za-z0-9_]*)>>", text)) + return placeholders + + +def validate_template_case(template_case: Path, calibration_token: str): + if not template_case.is_dir(): + raise FileNotFoundError( + f"Template case does not exist: {template_case}" + ) + + required_paths = [ + template_case / "0", + template_case / "constant" / "heatSourceDict", + template_case / "constant" / "scanPath", + template_case / "constant" / "transportProperties", + template_case / "system" / "controlDict", + template_case / "Allrun", + ] + missing = [ + str(path.relative_to(template_case)) + for path in required_paths + if not path.exists() + ] + if missing: + raise ValueError("Template case is missing: " + ", ".join(missing)) + + placeholder_locations = { + "constant/heatSourceDict": [calibration_token], + "constant/scanPath": ["power", "velocity"], + "system/controlDict": ["endTime", "writeInterval"], + } + for relative_path, required_tokens in placeholder_locations.items(): + text = (template_case / relative_path).read_text() + for token in required_tokens: + count = text.count(f"<<{token}>>") + if count != 1: + raise ValueError( + f"Template {relative_path} must contain <<{token}>> " + "exactly once; " + f"found {count}" + ) + + +def validate_additivefoam_environment(): + if os.environ.get("WM_PROJECT_VERSION") != "14": + raise RuntimeError( + "OpenFOAM-14 is required. Source OpenFOAM-14/etc/bashrc " + "before running." + ) + for variable in ("ADDITIVEFOAM_PROJECT_DIR", "ADDITIVEFOAM_ETC"): + if not os.environ.get(variable): + raise RuntimeError( + f"{variable} is not set. Source AdditiveFOAM/etc/bashrc " + "before running." + ) + if shutil.which("additiveFoam") is None: + raise RuntimeError("additiveFoam is not available on PATH") + + +def resolve_material_isovalue( + case_dir: Path, selector: str | float | int +) -> float: + if isinstance(selector, (float, int)): + return float(selector) + + selector_text = str(selector).strip().lower() + try: + return float(selector_text) + except ValueError: + pass + + if selector_text not in {"solidus", "liquidus"}: + raise ValueError( + "case.melt_pool_isovalue must be 'solidus', 'liquidus', or " + "a temperature" + ) + + foam_dictionary = shutil.which("foamDictionary") + if foam_dictionary is None: + raise RuntimeError( + "foamDictionary is required to resolve material isovalues" + ) + + transport_properties = case_dir / "constant" / "transportProperties" + # Expanding the dictionary resolves material files included through + # $ADDITIVEFOAM_ETC before thermoPath is inspected. + completed = subprocess.run( + [ + foam_dictionary, + "-expand", + "-entry", + "thermoPath", + "-value", + str(transport_properties), + ], + cwd=case_dir, + text=True, + stdout=subprocess.PIPE, + stderr=subprocess.STDOUT, + ) + if completed.returncode != 0: + raise RuntimeError( + "Could not read thermoPath from transportProperties:\n" + + completed.stdout + ) + + number = r"[-+]?(?:\d+(?:\.\d*)?|\.\d+)(?:[eE][-+]?\d+)?" + points = [ + (float(temperature), float(fraction)) + for temperature, fraction in re.findall( + rf"\(\s*({number})\s+({number})\s*\)", completed.stdout + ) + ] + target_fraction = 0.0 if selector_text == "liquidus" else 1.0 + for temperature, fraction in points: + if abs(fraction - target_fraction) <= 1.0e-12: + return temperature + + raise RuntimeError( + f"thermoPath does not contain an exact {selector_text} point " + f"with alpha.solid={target_fraction:g}" + ) + + +def compute_end_time_from_scan_path(scan_path_file: Path) -> float: + rows = [] + for line in scan_path_file.read_text().splitlines(): + line = line.strip() + if not line or line.startswith("//") or line.lower().startswith("mode"): + continue + fields = line.split() + if len(fields) < 6: + continue + rows.append([float(v) for v in fields[:6]]) + + end_time = 0.0 + for i in range(1, len(rows)): + mode = int(rows[i][0]) + if mode != 0: + continue + x0, y0, z0 = rows[i - 1][1:4] + x1, y1, z1 = rows[i][1:4] + velocity = rows[i][5] + distance = math.sqrt((x1 - x0) ** 2 + (y1 - y0) ** 2 + (z1 - z0) ** 2) + end_time += distance / velocity + return end_time + + +def find_melt_pool_dimensions_file(case_dir: Path, isovalue: float) -> Path: + output_dir = case_dir / "postProcessing" / "meltPoolDimensions" + candidates = [] + for path in output_dir.glob("*.csv"): + try: + candidates.append((float(path.stem), path)) + except ValueError: + continue + + tolerance = max(1.0e-9, abs(isovalue) * 1.0e-9) + matches = [ + path for value, path in candidates if abs(value - isovalue) <= tolerance + ] + if len(matches) == 1: + return matches[0] + + available = ", ".join(path.name for _, path in candidates) or "none" + raise FileNotFoundError( + f"Expected meltPoolDimensions output for {isovalue:g} K; " + f"available: {available}" + ) + + +def parse_melt_pool_dimensions( + case_dir: Path, selector: str | float | int +) -> dict: + isovalue = resolve_material_isovalue(case_dir, selector) + csv_path = find_melt_pool_dimensions_file(case_dir, isovalue) + df = pd.read_csv(csv_path) + if df.empty: + raise ValueError(f"Melt-pool dimensions output is empty: {csv_path}") + # Calibration uses the maximum penetration recorded during the scan. + row = df.loc[df["depth(m)"].idxmax()] + return { + "isovalue_K": isovalue, + "melt_pool_dimensions_file_name": csv_path.name, + "time_s": float(row["time(s)"]), + "length_microns": float(row["length(m)"]) * 1e6, + "width_microns": float(row["width(m)"]) * 1e6, + "depth_microns": float(row["depth(m)"]) * 1e6, + } + + +def _looks_like_time_dir(path: Path) -> bool: + if not path.is_dir() or path.name == "0": + return False + try: + float(path.name) + return True + except ValueError: + return False + + +def finalize_run(case_dir: Path, metrics: dict, keep_case: bool): + save_yaml_file(metrics, case_dir / "metrics.yml") + + if keep_case: + return + + for path in case_dir.iterdir(): + if path.is_dir() and ( + path.name.startswith("processor") or _looks_like_time_dir(path) + ): + shutil.rmtree(path, ignore_errors=True) + + +def format_case_value(value, decimals: int | None = None) -> str: + value = float(value) + if decimals is not None: + text = f"{value:.{decimals}f}" + else: + text = f"{value:.10g}" + text = text.rstrip("0").rstrip(".") if "." in text else text + text = text.replace("-", "m").replace(".", "p") + return text + + +def make_condition_key(params: dict) -> str: + power = format_case_value(params["Power_W"]) + speed = format_case_value(params["Speed_mm_s"]) + spot = format_case_value(params["Spot_Size_microns"]) + return f"P{power}_V{speed}_D{spot}" + + +def make_run_key(B: float) -> str: + B_key = format_case_value(B) + return f"B{B_key}" + + +def run_additivefoam_case(case_dir: Path, run_command: str) -> int: + completed = subprocess.run( + run_command, + cwd=case_dir, + shell=True, + text=True, + stdout=subprocess.PIPE, + stderr=subprocess.STDOUT, + ) + (case_dir / "log.run").write_text(completed.stdout) + return completed.returncode + + +def run_one_additivefoam_simulation( + params: dict, B: float, config: dict +) -> float: + validate_additivefoam_environment() + campaign_dir = Path(config["campaign_dir"]) + template_case = Path(config["template_case"]) + condition_key = make_condition_key(params) + run_key = make_run_key(B) + + case_dir = campaign_dir / "cases" / condition_key / run_key + if case_dir.exists(): + shutil.rmtree(case_dir) + shutil.copytree(template_case, case_dir) + + velocity_m_s = float(params["Speed_mm_s"]) * 1e-3 + values = { + config["calibration_token"]: float(B), + "power": float(params["Power_W"]), + "velocity": velocity_m_s, + } + + # Render the scan first because its geometry and velocity determine endTime. + render_case(case_dir, values) + end_time = compute_end_time_from_scan_path( + case_dir / "constant" / "scanPath" + ) + values["endTime"] = end_time + values["writeInterval"] = end_time + render_case(case_dir, values) + unresolved = sorted(find_template_placeholders(case_dir)) + if unresolved: + raise RuntimeError( + "Rendered case contains unresolved placeholders: " + + ", ".join(unresolved) + ) + + print( + f" Running AdditiveFOAM: B={B:.4f}, " + f"P={params['Power_W']}, v={velocity_m_s:.4g} m/s" + ) + return_code = run_additivefoam_case(case_dir, config["run_command"]) + if return_code != 0: + raise RuntimeError(f"AdditiveFOAM run failed in {case_dir}") + + metrics = parse_melt_pool_dimensions( + case_dir, + config["melt_pool_dimensions_isovalue"], + ) + metrics.update( + { + "parameters": params, + "B": float(B), + "power": float(params["Power_W"]), + "velocity": velocity_m_s, + "endTime": end_time, + "writeInterval": end_time, + } + ) + finalize_run( + case_dir, metrics, bool(config.get("keep_successful_cases", False)) + ) + return float(metrics["depth_microns"]) + + +# Local response surrogate and Bayesian calibration +def find_matching_simulation_subset( + sim_df: pd.DataFrame, params: dict +) -> pd.DataFrame: + if sim_df.empty: + return pd.DataFrame() + + mask = pd.Series(True, index=sim_df.index) + for key, value in params.items(): + if key not in sim_df.columns: + return pd.DataFrame() + mask &= sim_df[key] == value + return sim_df.loc[mask] + + +def get_sorted_simulation_curve(model_subset: pd.DataFrame): + curve = model_subset.sort_values("n") + n_coords = curve["n"].values.astype(float) + model_normalized_depths = curve["Normalized_Simulated_Depth"].values.astype( + float + ) + return n_coords, model_normalized_depths + + +def build_pchip_surrogate_grid( + n_coords: np.ndarray, + model_values: np.ndarray, + n_grid: int = 1000, +) -> tuple[np.ndarray, np.ndarray]: + order = np.argsort(n_coords) + n_sorted = np.asarray(n_coords, dtype=float)[order] + y_sorted = np.asarray(model_values, dtype=float)[order] + + unique_n, unique_idx = np.unique(n_sorted, return_index=True) + unique_y = y_sorted[unique_idx] + + if len(unique_n) < 3: + return unique_n, unique_y + + # PCHIP preserves the sampled curve shape without extrapolating beyond + # CFD data. + n_dense = np.linspace(unique_n.min(), unique_n.max(), int(n_grid)) + y_dense = PchipInterpolator(unique_n, unique_y, extrapolate=False)(n_dense) + return n_dense.astype(float), np.asarray(y_dense, dtype=float) + + +def pchip_grid_from_subset( + model_subset: pd.DataFrame, + value_column: str = "Normalized_Simulated_Depth", + n_grid: int = 1000, +) -> tuple[np.ndarray, np.ndarray]: + curve = model_subset.sort_values("n") + return build_pchip_surrogate_grid( + curve["n"].values.astype(float), + curve[value_column].values.astype(float), + n_grid=n_grid, + ) + + +def proposed_B_values_for_initial_curve( + model_subset: pd.DataFrame, config: dict +) -> list[float]: + existing = set() + if not model_subset.empty: + existing = {round(float(v), 8) for v in model_subset["n"].values} + + adaptive = config["adaptive_learning"] + B_values = adaptive.get("initial_B_values", [0.0, 4.5, 9.0]) + return [float(v) for v in B_values if round(float(v), 8) not in existing] + + +def propose_next_B( + model_subset: pd.DataFrame, target_depth: float, config: dict +): + adaptive = config["adaptive_learning"] + B_min = float(adaptive.get("B_min", 0.0)) + B_max = float(adaptive.get("B_max", 9.0)) + + B_values, depths = pchip_grid_from_subset( + model_subset, + value_column="Simulated_Depth_microns", + n_grid=int(adaptive.get("pchip_grid_points", 1000)), + ) + + for i in range(len(B_values) - 1): + d0, d1 = depths[i], depths[i + 1] + if (d0 <= target_depth <= d1) or (d1 <= target_depth <= d0): + if abs(d1 - d0) < 1e-12: + B = 0.5 * (B_values[i] + B_values[i + 1]) + else: + B = B_values[i] + (target_depth - d0) * ( + B_values[i + 1] - B_values[i] + ) / (d1 - d0) + B = min(max(float(B), B_min), B_max) + existing = model_subset["n"].values.astype(float) + if np.min(np.abs(existing - B)) > 1e-5: + return B + return 0.5 * (B_values[i] + B_values[i + 1]) + return None + + +def ensure_simulation_curve( + job_row, sim_df: pd.DataFrame, raw_simulations: list, config: dict +): + if not config["adaptive_learning"].get("enabled", True): + return sim_df, raw_simulations + + params = job_row["parameters"] + model_subset = find_matching_simulation_subset(sim_df, params) + # A missing condition starts with the configured B grid. Existing cache + # entries are skipped, so interrupted campaigns resume without rerunning + # them. + if len(model_subset) < 2: + for B in proposed_B_values_for_initial_curve(model_subset, config): + depth = run_one_additivefoam_simulation(params, B, config) + append_simulation_result( + raw_simulations, params, B, depth, config["simulations_file"] + ) + raw_simulations = load_yaml_file(config["simulations_file"]) + sim_df = parse_simulations(raw_simulations) + model_subset = find_matching_simulation_subset(sim_df, params) + return sim_df, raw_simulations + + +def linear_interp_pt(n_val, n_data_pt, column_data_pt): + n_data_pt, column_data_pt, n_val = ( + pt.cast(n_data_pt, "float64"), + pt.cast(column_data_pt, "float64"), + pt.cast(n_val, "float64"), + ) + idx, data_len = pt.searchsorted(n_data_pt, n_val), pt.shape(n_data_pt)[0] + idx = pt.clip(idx, 1, data_len - 1) + idx_lower, idx_upper = idx - 1, idx + n_lower, n_upper = n_data_pt[idx_lower], n_data_pt[idx_upper] + val_lower, val_upper = column_data_pt[idx_lower], column_data_pt[idx_upper] + return val_lower + (val_upper - val_lower) * (n_val - n_lower) / ( + n_upper - n_lower + ) + + +def perform_bayesian_calibration( + n_coords: np.ndarray, + model_values: np.ndarray, + observed_values: list[float], + config: dict, +) -> az.InferenceData: + + # Measurement scatter sets the likelihood width, with a 5% floor for + # repeated measurements that are identical or nearly identical. + sigma_est = max(0.05 * np.mean(observed_values), np.std(observed_values)) + + bayes = config.get("bayesian", {}) + with pm.Model() as model: + n = pm.Uniform("n", lower=n_coords.min(), upper=n_coords.max()) + n_data_pt = pt.constant(n_coords) + model_values_pt = pt.constant(model_values) + + predicted_value = pm.Deterministic( + "predicted_value", linear_interp_pt(n, n_data_pt, model_values_pt) + ) + pm.Normal( + "likelihood", + mu=predicted_value, + sigma=sigma_est, + observed=observed_values, + ) + print( + f" Sampling posterior with {len(observed_values)} observations " + f"(sigma_est={sigma_est:.3f})..." + ) + trace = pm.sample( + draws=int(bayes.get("draws", 2000)), + tune=int(bayes.get("tune", 1000)), + cores=int(bayes.get("cores", 4)), + progressbar=bool(bayes.get("progressbar", True)), + target_accept=float(bayes.get("target_accept", 0.9)), + random_seed=int(bayes.get("random_seed", 42)), + ) + return trace + + +def summarize_local_n_posterior(trace: az.InferenceData) -> dict: + n_samples = trace.posterior["n"].values.flatten() + + # Use the sampled point with highest log density as the reported local + # value. + if "lp" in trace.sample_stats: + logp = trace.sample_stats["lp"].values.flatten() + n_mode = float(n_samples[np.argmax(logp)]) + else: + counts, edges = np.histogram(n_samples, bins=100) + i = int(np.argmax(counts)) + n_mode = float(0.5 * (edges[i] + edges[i + 1])) + + q025, q975 = np.quantile(n_samples, [0.025, 0.975]) + n_variance = float(np.var(n_samples)) + + return { + "n_mode": n_mode, + "n_mean": float(np.mean(n_samples)), + "n_median": float(np.median(n_samples)), + "n_variance": n_variance, + "n_std": float(np.sqrt(n_variance)), + "n_95_low": float(q025), + "n_95_high": float(q975), + "samples": n_samples, + } + + +def maybe_refine_curve(job_row, trace, sim_df, raw_simulations, config): + adaptive = config["adaptive_learning"] + if not adaptive.get("enabled", True): + return trace, sim_df, raw_simulations + + params = job_row["parameters"] + model_subset = find_matching_simulation_subset(sim_df, params) + max_sims = int(adaptive.get("max_simulations_per_experiment", 7)) + posterior_tol = float(adaptive.get("posterior_std_tolerance", 0.25)) + depth_tol = float(adaptive.get("depth_tolerance_microns", 10.0)) + + # Add CFD points only while both the posterior and depth match remain broad. + while len(model_subset) < max_sims: + n_samples = trace.posterior["n"].values.flatten() + if np.std(n_samples) <= posterior_tol: + break + + target_depth = float(np.mean(job_row["depths_list"])) + best_error = np.min( + np.abs( + model_subset["Simulated_Depth_microns"].values - target_depth + ) + ) + if best_error <= depth_tol: + break + + B_next = propose_next_B(model_subset, target_depth, config) + if B_next is None: + break + + depth = run_one_additivefoam_simulation(params, B_next, config) + append_simulation_result( + raw_simulations, + params, + B_next, + depth, + config["simulations_file"], + ) + raw_simulations = load_yaml_file(config["simulations_file"]) + sim_df = parse_simulations(raw_simulations) + model_subset = find_matching_simulation_subset(sim_df, params) + + n_coords, model_normalized_depths = get_sorted_simulation_curve( + model_subset + ) + n_surrogate, depth_surrogate = build_pchip_surrogate_grid( + n_coords, + model_normalized_depths, + n_grid=int(adaptive.get("pchip_grid_points", 1000)), + ) + trace = perform_bayesian_calibration( + n_surrogate, + depth_surrogate, + job_row["normalized_depths_list"], + config, + ) + + return trace, sim_df, raw_simulations + + +# Global response fit, uncertainty propagation, and reporting +@dataclass(frozen=True) +class ReportTheme: + page_size: tuple = letter + + margin_left: float = 0.65 * inch + margin_right: float = 0.65 * inch + margin_top: float = 0.65 * inch + margin_bottom: float = 0.65 * inch + + ornl_green: str = "#4B7F2A" + dark_green: str = "#2F5D1E" + gray: str = "#555555" + light_gray: str = "#E8E8E8" + very_light_gray: str = "#F7F7F7" + black: str = "#111111" + red: str = "#B22222" + blue: str = "#1F77B4" + orange: str = "#D55E00" + + body_font: str = "Helvetica" + body_bold_font: str = "Helvetica-Bold" + mono_font: str = "Courier" + + plot_dpi: int = 220 + + +REPORT_THEME = ReportTheme() + + +def _escape_reportlab_text(value: Any) -> str: + text = str(value) + return text.replace("&", "&").replace("<", "<").replace(">", ">") + + +def condition_label(params_dict: dict) -> str: + if not isinstance(params_dict, dict): + return str(params_dict) + + return ", ".join( + f"{str(k).replace('_', ' ')}={v}" for k, v in params_dict.items() + ) + + +def wrapped_condition_label(params_dict: dict, width: int = 58) -> str: + label = condition_label(params_dict) + return "\n".join( + textwrap.wrap( + label, + width=width, + break_long_words=False, + break_on_hyphens=False, + ) + ) + + +def _parameter_key(params: dict) -> tuple: + if not isinstance(params, dict): + return tuple() + return tuple(sorted(params.items())) + + +def _parse_param_key(param_key: Any) -> dict: + if isinstance(param_key, dict): + return param_key + + if isinstance(param_key, (list, tuple)): + try: + return dict(param_key) + except Exception: + return {"condition": str(param_key)} + + if isinstance(param_key, str): + try: + parsed = ast.literal_eval(param_key) + if isinstance(parsed, dict): + return parsed + if isinstance(parsed, (list, tuple)): + return dict(parsed) + except Exception: + pass + + return {"condition": str(param_key)} + + +def filter_simulations_by_params( + simulations_df: pd.DataFrame, + params_dict: dict, +) -> pd.DataFrame: + if simulations_df.empty: + return simulations_df.copy() + + mask = pd.Series(True, index=simulations_df.index) + for key, value in params_dict.items(): + if key not in simulations_df.columns: + return simulations_df.iloc[0:0].copy() + mask &= simulations_df[key] == value + + return simulations_df.loc[mask].sort_values("n") + + +def report_plot_context(): + return mpl.rc_context( + { + "figure.dpi": REPORT_THEME.plot_dpi, + "savefig.dpi": REPORT_THEME.plot_dpi, + "font.family": "DejaVu Sans", + "font.size": 10, + "axes.titlesize": 10, + "axes.labelsize": 10, + "xtick.labelsize": 10, + "ytick.labelsize": 10, + "legend.fontsize": 10, + "axes.grid": True, + "grid.color": REPORT_THEME.black, + "grid.linestyle": "-", + "grid.linewidth": 0.6, + "axes.edgecolor": REPORT_THEME.black, + "axes.linewidth": 0.8, + "figure.constrained_layout.use": True, + "pdf.fonttype": 42, + "ps.fonttype": 42, + } + ) + + +def fig_to_reportlab_image( + fig: mpl.figure.Figure, + max_width: float, + max_height: float | None = None, + dpi: int = REPORT_THEME.plot_dpi, +) -> Image: + buffer = BytesIO() + fig.savefig( + buffer, + format="png", + dpi=dpi, + bbox_inches="tight", + pad_inches=0.12, + ) + buffer.seek(0) + plt.close(fig) + + img = Image(buffer) + aspect = img.imageHeight / float(img.imageWidth) + img.drawWidth = max_width + img.drawHeight = max_width * aspect + + if max_height is not None and img.drawHeight > max_height: + scale = max_height / img.drawHeight + img.drawWidth *= scale + img.drawHeight *= scale + + return img + + +def plot_calibration_overview_pages( + results_df: pd.DataFrame, + simulations_df: pd.DataFrame, + plots_per_page: int = 1, +) -> list[mpl.figure.Figure]: + if results_df.empty: + return [] + + figs: list[mpl.figure.Figure] = [] + + with report_plot_context(): + for plot_idx, (_, row) in enumerate(results_df.iterrows(), start=1): + params_dict = row["parameters"] + sim_curve = filter_simulations_by_params( + simulations_df, params_dict + ) + + fig, ax = plt.subplots( + figsize=(10.5, 6.4), + layout="constrained", + ) + + if not sim_curve.empty: + sim_curve = sim_curve.sort_values("n") + + n_plot, y_plot = pchip_grid_from_subset( + sim_curve, + value_column="Normalized_Simulated_Depth", + n_grid=1000, + ) + + ax.plot( + n_plot, + y_plot, + color="black", + linewidth=1.6, + label="PCHIP simulation surrogate", + ) + + ax.scatter( + sim_curve["n"], + sim_curve["Normalized_Simulated_Depth"], + s=26, + color="black", + zorder=3, + ) + + calibrated_depth = float( + np.interp( + float(row["calibrated_n"]), + n_plot, + y_plot, + ) + ) + else: + calibrated_depth = float(row["mean_normalized_depth"]) + + if "calibrated_n_95_low" in row and "calibrated_n_95_high" in row: + interval_low = min( + float(row["calibrated_n_95_low"]), + float(row["calibrated_n_95_high"]), + ) + interval_high = max( + float(row["calibrated_n_95_low"]), + float(row["calibrated_n_95_high"]), + ) + ax.plot( + [interval_low, interval_high], + [calibrated_depth, calibrated_depth], + color=REPORT_THEME.red, + marker="|", + markersize=9, + linewidth=1.2, + label="Local n with 95% interval", + zorder=4, + ) + xerr = None + calibration_label = "_nolegend_" + else: + xerr = row["calibrated_n_std"] + calibration_label = "Local n" + + ax.errorbar( + x=row["calibrated_n"], + y=calibrated_depth, + xerr=xerr, + fmt="o", + color=REPORT_THEME.red, + ecolor=REPORT_THEME.red, + capsize=4, + markersize=7, + linewidth=1.2, + label=calibration_label, + zorder=5, + ) + + depths = np.asarray(row["normalized_depths_list"], dtype=float) + depths = depths[np.isfinite(depths)] + + if len(depths) > 0: + ax.scatter( + np.full(len(depths), row["calibrated_n"]), + depths, + edgecolor=REPORT_THEME.blue, + facecolor="none", + linewidth=1.0, + s=42, + label="Experiment", + zorder=4, + ) + + ax.set_title( + wrapped_condition_label(params_dict, width=90), + fontsize=10, + fontweight="bold", + pad=10, + ) + + ax.set_xlabel("B used in A=0 local calibration") + ax.set_ylabel("Depth / spot size") + ax.margins(x=0.08, y=0.12) + + ax.legend( + loc="best", + frameon=True, + framealpha=0.92, + borderpad=0.45, + handlelength=1.6, + labelspacing=0.35, + ) + + ax.text( + 0.99, + 0.02, + f"Local calibration {plot_idx} of {len(results_df)}", + transform=ax.transAxes, + ha="right", + va="bottom", + fontsize=10, + color=REPORT_THEME.gray, + ) + + figs.append(fig) + + return figs + + +def plot_posterior_distributions( + posterior_samples_dict: dict, + max_legend_items: int = 10, +) -> mpl.figure.Figure | None: + if not posterior_samples_dict: + return None + + with report_plot_context(): + fig, ax = plt.subplots( + figsize=(10.5, 6.4), + layout="constrained", + ) + + color_values = plt.cm.viridis( + np.linspace(0.08, 0.92, len(posterior_samples_dict)) + ) + + plotted = 0 + + for i, (param_key, samples) in enumerate( + posterior_samples_dict.items() + ): + params_dict = _parse_param_key(param_key) + label = condition_label(params_dict) + + if i >= max_legend_items: + label = "_nolegend_" + + values = np.asarray(samples, dtype=float) + values = values[np.isfinite(values)] + + if len(values) == 0: + continue + + sns.histplot( + values, + label=label, + color=color_values[i], + kde=True, + alpha=0.38, + stat="probability", + ax=ax, + ) + + plotted += 1 + + ax.set_title("Local Calibration Distributions") + ax.set_xlabel("Calibrated n") + ax.set_ylabel("Probability") + ax.margins(x=0.04, y=0.10) + + handles, labels = ax.get_legend_handles_labels() + + if handles: + if len(posterior_samples_dict) > max_legend_items: + omitted = len(posterior_samples_dict) - max_legend_items + labels.append( + f"... {omitted} additional conditions omitted from legend" + ) + handles.append(mpl.lines.Line2D([], [], color="none")) + + ax.legend( + handles, + labels, + title="Process parameters", + loc="upper center", + bbox_to_anchor=(0.5, -0.16), + ncol=2, + frameon=False, + fontsize=10, + title_fontsize=10, + ) + + if plotted == 0: + plt.close(fig) + return None + + return fig + + +def _weighted_linear_initial_guess( + z_obs: np.ndarray, + n_obs: np.ndarray, + n_sigma_obs: np.ndarray, +) -> tuple[np.ndarray, np.ndarray, float]: + """ + Weighted least-squares initial guess for n = A*z + B. + + Returns theta = [A, B], approximate covariance, and reduced chi-square. + """ + X = np.column_stack([z_obs, np.ones_like(z_obs)]) + weights = 1.0 / np.maximum(n_sigma_obs, 1.0e-12) ** 2 + + sqrt_w = np.sqrt(weights) + Xw = X * sqrt_w[:, None] + yw = n_obs * sqrt_w + + theta, *_ = np.linalg.lstsq(Xw, yw, rcond=None) + + residual = (n_obs - X @ theta) / np.maximum(n_sigma_obs, 1.0e-12) + dof = max(len(n_obs) - 2, 1) + reduced_chi2 = float(np.sum(residual**2) / dof) + + xtwx = Xw.T @ Xw + try: + covariance = np.linalg.inv(xtwx) * max(reduced_chi2, 1.0) + except np.linalg.LinAlgError: + covariance = np.linalg.pinv(xtwx) * max(reduced_chi2, 1.0) + + return theta.astype(float), covariance.astype(float), reduced_chi2 + + +def _robust_weighted_curve_fit( + z_obs: np.ndarray, + n_obs: np.ndarray, + n_sigma_obs: np.ndarray, + config: dict, +) -> dict: + """ + Fit the local n estimates with uncertainty weighting and a robust loss. + + The robust loss limits the influence of isolated conditions that are + inconsistent with the global response. + """ + final_fit_cfg = config.get("final_fit", {}) + bayes_cfg = config.get("bayesian", {}) + + sigma_floor = float( + final_fit_cfg.get( + "local_n_sigma_floor", + bayes_cfg.get("local_n_sigma_floor", 1.0e-6), + ) + ) + n_sigma_obs = np.maximum(np.asarray(n_sigma_obs, dtype=float), sigma_floor) + + theta0, covariance0, reduced_chi2_0 = _weighted_linear_initial_guess( + z_obs, + n_obs, + n_sigma_obs, + ) + + loss = str(final_fit_cfg.get("loss", "soft_l1")) + f_scale = float(final_fit_cfg.get("f_scale", 1.0)) + + def residual(theta): + A, B = theta + return (n_obs - (A * z_obs + B)) / n_sigma_obs + + used_scipy = False + if least_squares is not None and len(n_obs) >= 3: + sol = least_squares( + residual, + x0=theta0, + loss=loss, + f_scale=f_scale, + ) + theta = sol.x.astype(float) + used_scipy = True + else: + theta = theta0.astype(float) + loss = "linear_weighted_lstsq" + + A, B = theta + weighted_residual = residual(theta) + dof = max(len(n_obs) - 2, 1) + reduced_chi2 = float(np.sum(weighted_residual**2) / dof) + + X = np.column_stack([z_obs, np.ones_like(z_obs)]) + J = X / n_sigma_obs[:, None] + jt_j = J.T @ J + try: + covariance = np.linalg.inv(jt_j) * max(reduced_chi2, 1.0) + except np.linalg.LinAlgError: + covariance = np.linalg.pinv(jt_j) * max(reduced_chi2, 1.0) + + return { + "A": float(A), + "B": float(B), + "A_std": float(np.sqrt(max(covariance[0, 0], 0.0))), + "B_std": float(np.sqrt(max(covariance[1, 1], 0.0))), + "covariance": covariance.tolist(), + "weighted_residuals": weighted_residual.astype(float).tolist(), + "reduced_chi2": reduced_chi2, + "initial_A": float(theta0[0]), + "initial_B": float(theta0[1]), + "initial_reduced_chi2": reduced_chi2_0, + "loss": loss, + "f_scale": f_scale, + "used_scipy_least_squares": bool(used_scipy), + "local_n_sigma_floor": sigma_floor, + } + + +def _posterior_samples_for_results( + results_df: pd.DataFrame, + posterior_samples_dict: dict | None, +) -> list[np.ndarray | None]: + """ + Return local calibration samples aligned with rows in results_df. + + The posterior_samples_dict keys are expected to be + str(sorted(params.items())). Missing samples are returned as None; + downstream uncertainty propagation then falls back to a Gaussian + approximation using the stored calibrated value and standard uncertainty. + """ + samples_by_row: list[np.ndarray | None] = [] + posterior_samples_dict = posterior_samples_dict or {} + + for _, row in results_df.iterrows(): + key = str(sorted(row["parameters"].items())) + samples = posterior_samples_dict.get(key) + if samples is None: + samples_by_row.append(None) + continue + + values = np.asarray(samples, dtype=float) + values = values[np.isfinite(values)] + if len(values) == 0: + samples_by_row.append(None) + else: + samples_by_row.append(values) + + return samples_by_row + + +def _bootstrap_robust_weighted_fit_band( + z_obs: np.ndarray, + n_mode_obs: np.ndarray, + n_sigma_obs: np.ndarray, + x_range: np.ndarray, + posterior_samples_by_row: list[np.ndarray | None], + config: dict, + x_log_floor: float, +) -> dict | None: + """ + Empirical 95% response-curve interval from local calibration uncertainty. + + Each realization samples one n value from each local posterior, refits A and + B, and evaluates n = clip(A * log2(x) + B, 0, 9). Pointwise quantiles give a + generally asymmetric interval. + """ + final_fit_cfg = config.get("final_fit", {}) + n_bootstrap = int(final_fit_cfg.get("bootstrap_samples", 1000)) + seed = int(final_fit_cfg.get("bootstrap_random_seed", 12345)) + min_success = int( + final_fit_cfg.get("bootstrap_min_success", max(50, n_bootstrap // 10)) + ) + + if n_bootstrap <= 0 or len(n_mode_obs) < 2: + return None + + rng = np.random.default_rng(seed) + curves = [] + A_samples = [] + B_samples = [] + + z_range = np.log2(np.maximum(x_range, x_log_floor)) + + # Each realization propagates the local posteriors through a complete + # refit, retaining covariance between A and B in the response band. + for _ in range(n_bootstrap): + sampled_n = np.empty_like(n_mode_obs, dtype=float) + + for i, samples in enumerate(posterior_samples_by_row): + if samples is not None and len(samples) > 0: + sampled_n[i] = float(rng.choice(samples)) + else: + sigma = max(float(n_sigma_obs[i]), 1.0e-8) + sampled_n[i] = float(rng.normal(float(n_mode_obs[i]), sigma)) + + sampled_n = np.clip(sampled_n, 0.0, 9.0) + + try: + fit_b = _robust_weighted_curve_fit( + z_obs=z_obs, + n_obs=sampled_n, + n_sigma_obs=n_sigma_obs, + config=config, + ) + except Exception: + continue + + A_b = float(fit_b["A"]) + B_b = float(fit_b["B"]) + if not np.isfinite(A_b) or not np.isfinite(B_b): + continue + + curve = np.clip(A_b * z_range + B_b, 0.0, 9.0) + curves.append(curve) + A_samples.append(A_b) + B_samples.append(B_b) + + if len(curves) < min_success: + return None + + curves = np.asarray(curves, dtype=float) + A_samples = np.asarray(A_samples, dtype=float) + B_samples = np.asarray(B_samples, dtype=float) + + lower, median, upper = np.quantile(curves, [0.025, 0.5, 0.975], axis=0) + + return { + "x": x_range.astype(float).tolist(), + "lower": lower.astype(float).tolist(), + "median": median.astype(float).tolist(), + "upper": upper.astype(float).tolist(), + "A_samples": A_samples.astype(float).tolist(), + "B_samples": B_samples.astype(float).tolist(), + "A_95_low": float(np.quantile(A_samples, 0.025)), + "A_50": float(np.quantile(A_samples, 0.5)), + "A_95_high": float(np.quantile(A_samples, 0.975)), + "B_95_low": float(np.quantile(B_samples, 0.025)), + "B_50": float(np.quantile(B_samples, 0.5)), + "B_95_high": float(np.quantile(B_samples, 0.975)), + "n_successful_bootstrap": int(len(curves)), + "n_requested_bootstrap": int(n_bootstrap), + "random_seed": int(seed), + } + + +def plot_final_fit( + x_obs: np.ndarray, + n_obs: np.ndarray, + n_stds_obs: np.ndarray, + A_m: float, + B_m: float, + covariance: np.ndarray | None, + x_log_floor: float, + n_95_low_obs: np.ndarray | None = None, + n_95_high_obs: np.ndarray | None = None, + bootstrap_band: dict | None = None, +) -> mpl.figure.Figure: + """ + Plot the fitted projected-source response and its uncertainty. + + ``x_obs`` is measured depth divided by spot size. Values below one are + valid. + The curve and interval use n = clip(A * log2(x) + B, 0, 9). + """ + with report_plot_context(): + fig, ax = plt.subplots( + figsize=(10.5, 6.4), + layout="constrained", + ) + + if n_95_low_obs is not None and n_95_high_obs is not None: + interval_a = np.asarray(n_95_low_obs, dtype=float) + interval_b = np.asarray(n_95_high_obs, dtype=float) + lower = np.minimum(interval_a, interval_b) + upper = np.maximum(interval_a, interval_b) + yerr = np.vstack( + [ + np.maximum(n_obs - lower, 0.0), + np.maximum(upper - n_obs, 0.0), + ] + ) + point_label = "Local calibration with 95% interval" + else: + yerr = n_stds_obs + point_label = "Local calibration with standard uncertainty" + + ax.errorbar( + x_obs, + n_obs, + yerr=yerr, + fmt="o", + color=REPORT_THEME.blue, + ecolor=REPORT_THEME.blue, + capsize=3, + markersize=5, + linewidth=1.0, + label=point_label, + zorder=4, + ) + + x_min = max(float(np.nanmin(x_obs)) * 0.95, x_log_floor) + x_max = float(np.nanmax(x_obs)) * 1.05 + x_range = np.geomspace(x_min, x_max, 300) + z_range = np.log2(np.maximum(x_range, x_log_floor)) + + n_raw_pred = A_m * z_range + B_m + mean_pred = np.clip(n_raw_pred, 0.0, 9.0) + + ax.plot( + x_range, + mean_pred, + color=REPORT_THEME.orange, + linewidth=2.0, + label=f"Recommended calibration: n = {A_m:.4f} log2(x) + {B_m:.4f}", + zorder=3, + ) + + plotted_bootstrap_band = False + if bootstrap_band is not None: + required = {"x", "lower", "upper"} + if required.issubset(set(bootstrap_band.keys())): + x_boot = np.asarray(bootstrap_band["x"], dtype=float) + lower_boot = np.asarray(bootstrap_band["lower"], dtype=float) + upper_boot = np.asarray(bootstrap_band["upper"], dtype=float) + if ( + len(x_boot) == len(lower_boot) == len(upper_boot) + and len(x_boot) > 1 + and np.all(np.isfinite(x_boot)) + ): + ax.fill_between( + x_boot, + lower_boot, + upper_boot, + color=REPORT_THEME.orange, + alpha=0.22, + label="Calibrated 95% fit interval", + zorder=2, + ) + plotted_bootstrap_band = True + + if ( + not plotted_bootstrap_band + and covariance is not None + and np.all(np.isfinite(covariance)) + ): + design = np.column_stack([z_range, np.ones_like(z_range)]) + pred_var = np.einsum("ij,jk,ik->i", design, covariance, design) + pred_std = np.sqrt(np.maximum(pred_var, 0.0)) + lower_pred = np.clip(n_raw_pred - 1.96 * pred_std, 0.0, 9.0) + upper_pred = np.clip(n_raw_pred + 1.96 * pred_std, 0.0, 9.0) + ax.fill_between( + x_range, + lower_pred, + upper_pred, + color=REPORT_THEME.orange, + alpha=0.22, + label="Approximate 95% fit interval", + zorder=2, + ) + + ax.set_xscale("log", base=2) + ax.set_title("AdditiveFOAM Calibration Fit") + ax.set_xlabel("x = measured depth / spot size") + ax.set_ylabel("Shape factor n") + ax.margins(x=0.04, y=0.12) + + ax.legend( + loc="best", + frameon=True, + framealpha=0.92, + borderpad=0.45, + handlelength=1.6, + ) + + return fig + + +def fit_and_plot_heteroskedastic_model( + calibrated_results_df: pd.DataFrame, + config: dict, + posterior_samples_dict: dict | None = None, +): + """ + Fit and plot the production projected-source relation. + + Each local calibration contributes an n estimate and its uncertainty. The + robust, uncertainty-weighted relation is + + n = clip(A * log2(x) + B, 0, 9) + + where x is measured depth divided by spot size. The reported 95% response + interval is propagated from the local calibration uncertainty. + """ + if calibrated_results_df.empty or len(calibrated_results_df) < 2: + return None, None + + x_raw = calibrated_results_df["mean_normalized_depth"].values.astype(float) + + bayes = config.get("bayesian", {}) + report_cfg = config.get("report", {}) + x_log_floor = float( + report_cfg.get("x_log_floor", bayes.get("x_log_floor", 1.0e-12)) + ) + + if np.any(~np.isfinite(x_raw)): + raise ValueError("mean_normalized_depth contains non-finite values.") + + if np.any(x_raw <= 0.0): + print( + "\nWARNING: Some mean_normalized_depth values are <= 0. " + f"These will be clipped to x_log_floor={x_log_floor:g} for log2(x)." + ) + + x_obs = np.maximum(x_raw, x_log_floor) + z_obs = np.log2(x_obs) + + n_obs = calibrated_results_df["calibrated_n"].values.astype(float) + + if "calibrated_n_variance" in calibrated_results_df.columns: + n_var_obs = calibrated_results_df[ + "calibrated_n_variance" + ].values.astype(float) + n_var_obs = np.where(np.isfinite(n_var_obs), n_var_obs, np.nan) + else: + n_var_obs = np.full_like(n_obs, np.nan, dtype=float) + + n_std_fallback = calibrated_results_df["calibrated_n_std"].values.astype( + float + ) + n_var_fallback = n_std_fallback**2 + n_var_obs = np.where( + np.isfinite(n_var_obs) & (n_var_obs >= 0.0), n_var_obs, n_var_fallback + ) + n_stds_obs = np.sqrt(np.maximum(n_var_obs, 0.0)) + + if "calibrated_n_95_low" in calibrated_results_df.columns: + n_95_low_obs = calibrated_results_df[ + "calibrated_n_95_low" + ].values.astype(float) + else: + n_95_low_obs = np.maximum(n_obs - 1.96 * n_stds_obs, 0.0) + + if "calibrated_n_95_high" in calibrated_results_df.columns: + n_95_high_obs = calibrated_results_df[ + "calibrated_n_95_high" + ].values.astype(float) + else: + n_95_high_obs = np.minimum(n_obs + 1.96 * n_stds_obs, 9.0) + + print("\n--- Performing final uncertainty-weighted calibration fit ---") + + fit = _robust_weighted_curve_fit( + z_obs=z_obs, + n_obs=n_obs, + n_sigma_obs=n_stds_obs, + config=config, + ) + + A_m = fit["A"] + B_m = fit["B"] + covariance = np.asarray(fit["covariance"], dtype=float) + + print("--- Recommended calibration: n = A*log2(x) + B ---") + print(f"A = {A_m:.6f}") + print(f"B = {B_m:.6f}") + print(f"n = clip({A_m:.6f}*log2(x) + {B_m:.6f}, 0, 9)") + + x_min_plot = max(float(np.nanmin(x_obs)) * 0.95, x_log_floor) + x_max_plot = float(np.nanmax(x_obs)) * 1.05 + x_range = np.geomspace(x_min_plot, x_max_plot, 300) + + posterior_samples_by_row = _posterior_samples_for_results( + calibrated_results_df, + posterior_samples_dict, + ) + bootstrap_band = _bootstrap_robust_weighted_fit_band( + z_obs=z_obs, + n_mode_obs=n_obs, + n_sigma_obs=n_stds_obs, + x_range=x_range, + posterior_samples_by_row=posterior_samples_by_row, + config=config, + x_log_floor=x_log_floor, + ) + + fit_fig = plot_final_fit( + x_obs=x_obs, + n_obs=n_obs, + n_stds_obs=n_stds_obs, + A_m=A_m, + B_m=B_m, + covariance=covariance, + x_log_floor=x_log_floor, + n_95_low_obs=n_95_low_obs, + n_95_high_obs=n_95_high_obs, + bootstrap_band=bootstrap_band, + ) + + x_min = float(np.nanmin(x_obs)) + x_max = float(np.nanmax(x_obs)) + n_raw_at_x_min = A_m * np.log2(x_min) + B_m + n_raw_at_x_max = A_m * np.log2(x_max) + B_m + + fit_summary = { + "fit_method": "Robust uncertainty-weighted response fit", + "A": float(A_m), + "B": float(B_m), + "A_std_covariance_approx": float(fit["A_std"]), + "B_std_covariance_approx": float(fit["B_std"]), + "covariance": fit["covariance"], + "loss": fit["loss"], + "f_scale": float(fit["f_scale"]), + "reduced_chi2": float(fit["reduced_chi2"]), + "initial_A": float(fit["initial_A"]), + "initial_B": float(fit["initial_B"]), + "initial_reduced_chi2": float(fit["initial_reduced_chi2"]), + "used_scipy_least_squares": bool(fit["used_scipy_least_squares"]), + "local_n_sigma_floor": float(fit["local_n_sigma_floor"]), + "x_min": x_min, + "x_max": x_max, + "x_log_floor": float(x_log_floor), + "n_raw_at_x_min": float(n_raw_at_x_min), + "n_raw_at_x_max": float(n_raw_at_x_max), + "n_applied_at_x_min": float(np.clip(n_raw_at_x_min, 0.0, 9.0)), + "n_applied_at_x_max": float(np.clip(n_raw_at_x_max, 0.0, 9.0)), + "uses_x_clipped_at_1": False, + "uses_n_clipped_to_0_9": True, + "weighted_residuals": fit["weighted_residuals"], + "local_interval_display": "calibrated value with 95% interval", + "fit_band_display": "calibrated empirical 95% interval", + } + + if bootstrap_band is not None: + fit_summary.update( + { + "bootstrap_samples_requested": int( + bootstrap_band["n_requested_bootstrap"] + ), + "bootstrap_samples_successful": int( + bootstrap_band["n_successful_bootstrap"] + ), + "bootstrap_random_seed": int(bootstrap_band["random_seed"]), + "A_bootstrap_95_low": float(bootstrap_band["A_95_low"]), + "A_bootstrap_median": float(bootstrap_band["A_50"]), + "A_bootstrap_95_high": float(bootstrap_band["A_95_high"]), + "B_bootstrap_95_low": float(bootstrap_band["B_95_low"]), + "B_bootstrap_median": float(bootstrap_band["B_50"]), + "B_bootstrap_95_high": float(bootstrap_band["B_95_high"]), + } + ) + else: + fit_summary.update( + { + "bootstrap_samples_requested": int( + config.get("final_fit", {}).get("bootstrap_samples", 1000) + ), + "bootstrap_samples_successful": 0, + "bootstrap_note": ( + "Empirical fit interval unavailable; plot fell back to " + "covariance approximation if possible." + ), + } + ) + + return fit_summary, fit_fig + + +def posterior_samples_for_state(trace, max_samples: int = 1000) -> list[float]: + samples = trace.posterior["n"].values.flatten() + + if len(samples) > max_samples: + idx = np.linspace(0, len(samples) - 1, max_samples).astype(int) + samples = samples[idx] + + return [float(v) for v in samples] + + +def build_report_dataframe( + exp_df: pd.DataFrame, + final_state_df: pd.DataFrame, +) -> pd.DataFrame: + rows = [] + + if final_state_df.empty: + return pd.DataFrame() + + final_state_df = final_state_df.copy() + final_state_df["_param_key"] = final_state_df["parameters"].apply( + _parameter_key + ) + + for _, exp_row in exp_df.iterrows(): + params = exp_row["parameters"] + param_key = _parameter_key(params) + + matches = final_state_df[final_state_df["_param_key"] == param_key] + if matches.empty: + continue + + state_row = matches.iloc[0] + + mean_depth = float(np.mean(exp_row["depths_list"])) + mean_normalized_depth = float( + np.mean(exp_row["normalized_depths_list"]) + ) + + calibrated_n = float(state_row["calibrated_n"]) + calibrated_n_std = float(state_row["calibrated_n_std"]) + + rows.append( + { + "parameters": params, + "calibrated_n": calibrated_n, + "calibrated_n_mode": float( + state_row.get("calibrated_n_mode", calibrated_n) + ), + "calibrated_n_mean": float( + state_row.get("calibrated_n_mean", calibrated_n) + ), + "calibrated_n_median": float( + state_row.get("calibrated_n_median", calibrated_n) + ), + "calibrated_n_variance": float( + state_row.get("calibrated_n_variance", calibrated_n_std**2) + ), + "calibrated_n_std": calibrated_n_std, + "calibrated_n_95_low": float( + state_row.get( + "calibrated_n_95_low", calibrated_n - calibrated_n_std + ) + ), + "calibrated_n_95_high": float( + state_row.get( + "calibrated_n_95_high", calibrated_n + calibrated_n_std + ) + ), + "mean_depth_microns": mean_depth, + "mean_normalized_depth": mean_normalized_depth, + "mean_x": mean_normalized_depth, + "normalized_depths_list": exp_row["normalized_depths_list"], + } + ) + + return pd.DataFrame(rows) + + +def save_summary_table(results_df: pd.DataFrame, filepath: Path): + if results_df.empty: + return + + rows = [] + + for _, row in results_df.iterrows(): + params = row["parameters"] + + rows.append( + { + **params, + "mean_depth_microns": row["mean_depth_microns"], + "x_depth_over_spot_size": row["mean_x"], + "calibrated_n": row["calibrated_n"], + "calibrated_n_mode": row.get( + "calibrated_n_mode", row["calibrated_n"] + ), + "calibrated_n_mean": row.get( + "calibrated_n_mean", row["calibrated_n"] + ), + "calibrated_n_median": row.get( + "calibrated_n_median", row["calibrated_n"] + ), + "calibrated_n_variance": row.get( + "calibrated_n_variance", + row["calibrated_n_std"] ** 2, + ), + "calibrated_n_std": row["calibrated_n_std"], + "calibrated_n_95_low": row.get( + "calibrated_n_95_low", + row["calibrated_n"] - row["calibrated_n_std"], + ), + "calibrated_n_95_high": row.get( + "calibrated_n_95_high", + row["calibrated_n"] + row["calibrated_n_std"], + ), + } + ) + + pd.DataFrame(rows).to_csv(filepath, index=False) + + +def build_report_styles(): + styles = getSampleStyleSheet() + + styles.add( + ParagraphStyle( + name="ReportTitle", + parent=styles["Title"], + fontName=REPORT_THEME.body_bold_font, + fontSize=18, + leading=22, + textColor=colors.HexColor(REPORT_THEME.dark_green), + spaceAfter=14, + alignment=TA_LEFT, + ) + ) + + styles.add( + ParagraphStyle( + name="SectionHeading", + parent=styles["Heading2"], + fontName=REPORT_THEME.body_bold_font, + fontSize=12, + leading=15, + textColor=colors.HexColor(REPORT_THEME.dark_green), + spaceBefore=12, + spaceAfter=6, + ) + ) + + styles.add( + ParagraphStyle( + name="BodySmall", + parent=styles["BodyText"], + fontName=REPORT_THEME.body_font, + fontSize=10, + leading=11, + textColor=colors.HexColor(REPORT_THEME.black), + spaceAfter=5, + ) + ) + + styles.add( + ParagraphStyle( + name="BodySmallBold", + parent=styles["BodyText"], + fontName=REPORT_THEME.body_bold_font, + fontSize=10, + leading=11, + textColor=colors.HexColor(REPORT_THEME.black), + spaceAfter=5, + ) + ) + + styles.add( + ParagraphStyle( + name="MonoBlock", + parent=styles["BodyText"], + fontName=REPORT_THEME.mono_font, + fontSize=10, + leading=10, + textColor=colors.HexColor(REPORT_THEME.black), + leftIndent=10, + spaceBefore=4, + spaceAfter=8, + ) + ) + + styles.add( + ParagraphStyle( + name="Caption", + parent=styles["BodyText"], + fontName=REPORT_THEME.body_font, + fontSize=10, + leading=9, + textColor=colors.HexColor(REPORT_THEME.black), + alignment=TA_CENTER, + spaceBefore=4, + spaceAfter=10, + ) + ) + + return styles + + +def make_key_value_table( + items: list[tuple[str, Any]], + col_widths: list[float], +): + styles = build_report_styles() + + data = [] + for key, value in items: + key_text = _escape_reportlab_text(key) + value_text = _escape_reportlab_text(value) + data.append( + [ + Paragraph(f"{key_text}", styles["BodySmall"]), + Paragraph(value_text, styles["BodySmall"]), + ] + ) + + table = Table(data, colWidths=col_widths, hAlign="LEFT") + table.setStyle( + TableStyle( + [ + ("VALIGN", (0, 0), (-1, -1), "TOP"), + ("LEFTPADDING", (0, 0), (-1, -1), 4), + ("RIGHTPADDING", (0, 0), (-1, -1), 4), + ("TOPPADDING", (0, 0), (-1, -1), 3), + ("BOTTOMPADDING", (0, 0), (-1, -1), 3), + ( + "GRID", + (0, 0), + (-1, -1), + 0.25, + colors.HexColor(REPORT_THEME.light_gray), + ), + ("BACKGROUND", (0, 0), (0, -1), colors.HexColor("#F5F5F5")), + ] + ) + ) + return table + + +def report_footer(canvas, doc): + canvas.saveState() + + width, _ = doc.pagesize + + canvas.setStrokeColor(colors.HexColor(REPORT_THEME.light_gray)) + canvas.setLineWidth(0.5) + canvas.line( + doc.leftMargin, + 0.45 * inch, + width - doc.rightMargin, + 0.45 * inch, + ) + + canvas.setFont(REPORT_THEME.body_font, 7) + canvas.setFillColor(colors.HexColor(REPORT_THEME.gray)) + canvas.drawString( + doc.leftMargin, + 0.28 * inch, + "Oak Ridge National Laboratory", + ) + canvas.drawRightString( + width - doc.rightMargin, + 0.28 * inch, + f"Page {doc.page}", + ) + + canvas.restoreState() + + +def write_report( + results_df: pd.DataFrame, + simulations_df: pd.DataFrame, + posterior_samples_dict: dict, + fit_summary: dict | None, + fit_fig: mpl.figure.Figure | None, + config: dict, +): + output_pdf = config.get("report", {}).get( + "output_pdf", + "reports/calibration_report.pdf", + ) + + output_pdf = Path(output_pdf) + if not output_pdf.is_absolute(): + output_pdf = Path(config["campaign_dir"]) / output_pdf + + output_pdf.parent.mkdir(parents=True, exist_ok=True) + + summary_csv = output_pdf.parent / "calibration_summary.csv" + save_summary_table(results_df, summary_csv) + + styles = build_report_styles() + + doc = SimpleDocTemplate( + str(output_pdf), + pagesize=REPORT_THEME.page_size, + rightMargin=REPORT_THEME.margin_right, + leftMargin=REPORT_THEME.margin_left, + topMargin=REPORT_THEME.margin_top, + bottomMargin=REPORT_THEME.margin_bottom, + title="AdditiveFOAM Calibration Report", + author="Oak Ridge National Laboratory", + subject="AdditiveFOAM calibration", + ) + + story = [] + + story.append( + Paragraph("AdditiveFOAM Calibration Report", styles["ReportTitle"]) + ) + story.append( + Paragraph( + "Automated calibration report for AdditiveFOAM", + styles["BodySmall"], + ) + ) + story.append(Spacer(1, 0.15 * inch)) + + run_items = [ + ("Generated", datetime.now().strftime("%Y-%m-%d %H:%M")), + ("Experiments file", config.get("experiments_file", "not specified")), + ("Simulation cache", config.get("simulations_file", "not specified")), + ("Template case", config.get("template_case", "not specified")), + ("Campaign directory", config.get("campaign_dir", "not specified")), + ( + "Melt pool dimensions isovalue", + config.get("melt_pool_dimensions_isovalue", "not specified"), + ), + ("Number of calibrated experiments", len(results_df)), + ] + + if not results_df.empty: + run_items.extend( + [ + ( + "Range x = depth / spot size", + f"[{results_df['mean_x'].min():.4f}, " + f"{results_df['mean_x'].max():.4f}]", + ), + ] + ) + + story.append(Paragraph("Run Summary", styles["SectionHeading"])) + story.append( + make_key_value_table(run_items, col_widths=[2.15 * inch, 4.95 * inch]) + ) + story.append(Spacer(1, 0.2 * inch)) + + if fit_summary: + A_std = fit_summary.get("A_std_covariance_approx", float("nan")) + B_std = fit_summary.get("B_std_covariance_approx", float("nan")) + fit_items = [ + ("Calibration relation", f"n = A * log2(x) + B"), + ("A", f"{fit_summary['A']:.6f}"), + ("B", f"{fit_summary['B']:.6f}"), + ( + "A standard error", + f"{A_std:.6f}", + ), + ( + "B standard error", + f"{B_std:.6f}", + ), + ( + "A 95% interval", + f"[{fit_summary.get('A_bootstrap_95_low', float('nan')):.6f}, " + f"{fit_summary.get('A_bootstrap_95_high', float('nan')):.6f}]", + ), + ( + "B 95% interval", + f"[{fit_summary.get('B_bootstrap_95_low', float('nan')):.6f}, " + f"{fit_summary.get('B_bootstrap_95_high', float('nan')):.6f}]", + ), + ( + "Weighted residual statistic", + f"{fit_summary.get('reduced_chi2', float('nan')):.4f}", + ), + ( + "Calibration x range", + f"{fit_summary['x_min']:.4f} to {fit_summary['x_max']:.4f}", + ), + ("n at minimum x", f"{fit_summary['n_applied_at_x_min']:.6f}"), + ("n at maximum x", f"{fit_summary['n_applied_at_x_max']:.6f}"), + ] + + story.append(Paragraph("Calibration Fit", styles["SectionHeading"])) + story.append( + make_key_value_table( + fit_items, col_widths=[2.15 * inch, 4.95 * inch] + ) + ) + story.append(Spacer(1, 0.2 * inch)) + + if not results_df.empty: + story.append(PageBreak()) + + story.append( + Paragraph("Calibration Summary Table", styles["SectionHeading"]) + ) + + table_df = results_df.copy() + available_cols = [ + col + for col in [ + "mean_depth_microns", + "mean_x", + "calibrated_n", + "calibrated_n_95_low", + "calibrated_n_95_high", + ] + if col in table_df.columns + ] + + table_df = table_df[available_cols].round(6) + + header_map = { + "mean_depth_microns": "Target mean, microns", + "mean_x": "depth / spot size", + "calibrated_n": "n estimate", + "calibrated_n_95_low": "n, lower 95%", + "calibrated_n_95_high": "n, upper 95%", + } + + table_data = [[header_map.get(col, col) for col in available_cols]] + table_data.extend(table_df.values.tolist()) + + col_width_lookup = { + "mean_depth_microns": 1.65 * inch, + "mean_x": 1.65 * inch, + "calibrated_n": 1.35 * inch, + "calibrated_n_95_low": 1.20 * inch, + "calibrated_n_95_high": 1.20 * inch, + } + col_widths = [col_width_lookup[col] for col in available_cols] + + summary_table = Table( + table_data, + repeatRows=1, + hAlign="LEFT", + colWidths=col_widths, + ) + + summary_table.setStyle( + TableStyle( + [ + ( + "BACKGROUND", + (0, 0), + (-1, 0), + colors.HexColor(REPORT_THEME.dark_green), + ), + ("TEXTCOLOR", (0, 0), (-1, 0), colors.white), + ("FONTNAME", (0, 0), (-1, 0), REPORT_THEME.body_bold_font), + ("FONTSIZE", (0, 0), (-1, -1), 10), + ( + "GRID", + (0, 0), + (-1, -1), + 0.25, + colors.HexColor(REPORT_THEME.light_gray), + ), + ( + "ROWBACKGROUNDS", + (0, 1), + (-1, -1), + [ + colors.white, + colors.HexColor(REPORT_THEME.very_light_gray), + ], + ), + ("VALIGN", (0, 0), (-1, -1), "MIDDLE"), + ("ALIGN", (1, 1), (-1, -1), "RIGHT"), + ("LEFTPADDING", (0, 0), (-1, -1), 4), + ("RIGHTPADDING", (0, 0), (-1, -1), 4), + ("TOPPADDING", (0, 0), (-1, -1), 3), + ("BOTTOMPADDING", (0, 0), (-1, -1), 3), + ] + ) + ) + + story.append(summary_table) + + story.append(PageBreak()) + + max_plot_width = 7.1 * inch + max_plot_height = 6.3 * inch + figure_number = 1 + + overview_figs = plot_calibration_overview_pages(results_df, simulations_df) + + for i, fig in enumerate(overview_figs, start=1): + condition = condition_label(results_df.iloc[i - 1]["parameters"]) + condition_escaped = _escape_reportlab_text(condition) + + story.append( + Paragraph(f"Local Calibration {i}", styles["SectionHeading"]) + ) + story.append( + fig_to_reportlab_image( + fig, + max_width=max_plot_width, + max_height=max_plot_height, + dpi=REPORT_THEME.plot_dpi, + ) + ) + story.append( + Paragraph( + f"Figure {figure_number}: " + f"Local calibration for {condition_escaped}.", + styles["Caption"], + ) + ) + figure_number += 1 + story.append(PageBreak()) + + posterior_fig = plot_posterior_distributions(posterior_samples_dict) + + if posterior_fig is not None: + story.append( + Paragraph( + "Local Calibration Distributions", styles["SectionHeading"] + ) + ) + story.append( + fig_to_reportlab_image( + posterior_fig, + max_width=max_plot_width, + max_height=max_plot_height, + dpi=REPORT_THEME.plot_dpi, + ) + ) + story.append( + Paragraph( + f"Figure {figure_number}: " + "Probability distributions for the local calibrations.", + styles["Caption"], + ) + ) + figure_number += 1 + story.append(PageBreak()) + + if fit_fig is not None: + story.append(Paragraph("Calibration Fit", styles["SectionHeading"])) + story.append( + fig_to_reportlab_image( + fit_fig, + max_width=max_plot_width, + max_height=max_plot_height, + dpi=REPORT_THEME.plot_dpi, + ) + ) + story.append( + Paragraph( + f"Figure {figure_number}: " + "Final AdditiveFOAM calibration fit. Symbols show local n " + "estimates with 95% intervals. The solid curve is n = A " + "log2(x) + B, where x is measured depth divided by spot " + "size. The shaded region is the 95% interval propagated " + "from local calibration uncertainty.", + styles["Caption"], + ) + ) + figure_number += 1 + + doc.build( + story, + onFirstPage=report_footer, + onLaterPages=report_footer, + ) + + print(f"\nWrote report: {output_pdf}") + print(f"Wrote summary table: {summary_csv}") + + +def main(): + parser = argparse.ArgumentParser( + description=( + "Run an AdditiveFOAM heat-source calibration campaign and generate " + "the fitted projected-source closure." + ) + ) + parser.add_argument( + "--config", + default="config.yml", + help="Calibration configuration file (default: config.yml)", + ) + args = parser.parse_args() + config = load_config(args.config) + if not Path(config["experiments_file"]).is_file(): + raise FileNotFoundError( + f"Experiments file does not exist: {config['experiments_file']}" + ) + validate_template_case( + Path(config["template_case"]), + config["calibration_token"], + ) + + # Fingerprints allow unchanged conditions to reuse their saved posterior + # state. + raw_experiments = load_yaml_file(config["experiments_file"]) + raw_simulations = load_yaml_file(config["simulations_file"]) + current_state_data = load_yaml_file(config["calibration_state_file"]) + + exp_df = parse_experiments(raw_experiments) + sim_df = parse_simulations(raw_simulations) + state_df = pd.DataFrame(current_state_data) + + if exp_df.empty: + print("No experimental data found. Exiting.") + return + + to_process_list, fresh_states, state_lookup = [], [], {} + if not state_df.empty: + for _, row in state_df.iterrows(): + state_lookup[str(sorted(row["parameters"].items()))] = row.to_dict() + + for _, exp_row in exp_df.iterrows(): + param_key = str(sorted(exp_row["parameters"].items())) + if ( + param_key in state_lookup + and state_lookup[param_key].get("fingerprint") + == exp_row["fingerprint"] + ): + fresh_states.append(state_lookup[param_key]) + else: + to_process_list.append(exp_row) + + to_process_df = pd.DataFrame(to_process_list) + + print("\n--- Adaptive Run Summary ---") + print(f"Found {len(exp_df)} total experimental parameter sets.") + print(f"Found {len(fresh_states)} up-to-date calibrations.") + print(f"Found {len(to_process_df)} new or stale experiments to process.") + + newly_calibrated_states = [] + posterior_samples = {} + + # Stage 1: infer a local n independently for each process condition. + if not to_process_df.empty: + for _, job_row in to_process_df.iterrows(): + params = job_row["parameters"] + print(f"\n--- Processing parameters: {params} ---") + + sim_df, raw_simulations = ensure_simulation_curve( + job_row, + sim_df, + raw_simulations, + config, + ) + model_subset = find_matching_simulation_subset(sim_df, params) + if model_subset.empty: + print(" No matching simulation data found.") + continue + + n_coords, model_normalized_depths = get_sorted_simulation_curve( + model_subset + ) + n_surrogate, depth_surrogate = build_pchip_surrogate_grid( + n_coords, + model_normalized_depths, + n_grid=int( + config.get("adaptive_learning", {}).get( + "pchip_grid_points", 1000 + ) + ), + ) + + trace = perform_bayesian_calibration( + n_surrogate, + depth_surrogate, + job_row["normalized_depths_list"], + config, + ) + trace, sim_df, raw_simulations = maybe_refine_curve( + job_row, + trace, + sim_df, + raw_simulations, + config, + ) + posterior_summary = summarize_local_n_posterior(trace) + + calibrated_n = posterior_summary["n_mode"] + calibrated_n_std = posterior_summary["n_std"] + + new_state = { + "parameters": params, + "calibrated_n": float(calibrated_n), + "calibrated_n_mode": float(posterior_summary["n_mode"]), + "calibrated_n_mean": float(posterior_summary["n_mean"]), + "calibrated_n_median": float(posterior_summary["n_median"]), + "calibrated_n_variance": float(posterior_summary["n_variance"]), + "calibrated_n_std": float(calibrated_n_std), + "calibrated_n_95_low": float(posterior_summary["n_95_low"]), + "calibrated_n_95_high": float(posterior_summary["n_95_high"]), + "fingerprint": job_row["fingerprint"], + } + newly_calibrated_states.append(new_state) + samples = posterior_samples_for_state(trace) + new_state["posterior_n_samples"] = samples + posterior_samples[str(sorted(params.items()))] = samples + + print( + f" -> Local n estimate: {calibrated_n:.3f} " + f"(mean={posterior_summary['n_mean']:.3f}, " + f"std={calibrated_n_std:.3f})" + ) + + final_state = fresh_states + newly_calibrated_states + save_yaml_file(final_state, config["calibration_state_file"]) + final_state_df = pd.DataFrame(final_state) + + print("\n\n" + "=" * 30) + print(" FINAL CALIBRATION STATE ") + print("=" * 30) + if not final_state_df.empty: + columns = ["parameters", "calibrated_n", "calibrated_n_std"] + print(final_state_df[columns].to_string(index=False)) + else: + print("No calibrated results exist in the state file.") + + for state in final_state: + key = str(sorted(state["parameters"].items())) + if key not in posterior_samples and "posterior_n_samples" in state: + posterior_samples[key] = state["posterior_n_samples"] + + # Stage 2: fit A and B across the completed local calibrations. + report_df = build_report_dataframe(exp_df, final_state_df) + if not report_df.empty: + fit_summary, fit_fig = fit_and_plot_heteroskedastic_model( + report_df, + config, + posterior_samples, + ) + if fit_summary: + save_yaml_file( + fit_summary, + Path(config["campaign_dir"]) / "calibration_fit.yml", + ) + if config.get("report", {}).get("enabled", True): + write_report( + report_df, + sim_df, + posterior_samples, + fit_summary, + fit_fig, + config, + ) + + +if __name__ == "__main__": + main() diff --git a/etc/bashrc b/etc/bashrc index 106e370e..18c73084 100644 --- a/etc/bashrc +++ b/etc/bashrc @@ -77,6 +77,26 @@ if [ -d "$ADDITIVEFOAM_PROJECT_DIR/bin" ]; then fi +# AdditiveFOAM Python environment +additivefoamVenv="$ADDITIVEFOAM_PROJECT_DIR/.venv" + +if [ -d "$additivefoamVenv" ] +then + if [ ! -x "$additivefoamVenv/bin/python" ] || \ + [ ! -r "$additivefoamVenv/bin/activate" ] + then + echo "WARNING: The AdditiveFOAM Python environment is incomplete." >&2 + echo " Run $ADDITIVEFOAM_PROJECT_DIR/Allwmake to rebuild it." >&2 + elif [ "${VIRTUAL_ENV:-}" != "$additivefoamVenv" ] + then + . "$additivefoamVenv/bin/activate" || \ + echo "WARNING: Unable to activate $additivefoamVenv." >&2 + fi +fi + +unset additivefoamVenv + + # AdditiveFOAM directories export ADDITIVEFOAM_TUTORIALS="$ADDITIVEFOAM_PROJECT_DIR/tutorials" export ADDITIVEFOAM_APPLICATIONS="$ADDITIVEFOAM_PROJECT_DIR/applications" diff --git a/tutorials/nLightAFX/constant/nLightAFX.cfg b/etc/heatSources/nLightAFX-1000.cfg similarity index 84% rename from tutorials/nLightAFX/constant/nLightAFX.cfg rename to etc/heatSources/nLightAFX-1000.cfg index 2b61e957..1d493239 100644 --- a/tutorials/nLightAFX/constant/nLightAFX.cfg +++ b/etc/heatSources/nLightAFX-1000.cfg @@ -1,5 +1,5 @@ /*---------------------------------------------------------------------------*\ - nLight AFX-1000 heat source mode include file + ORNL-characterized nLight AFX-1000 heat source modes Values are from the ORNL-characterized AFX beam profiles. @@ -8,35 +8,20 @@ radius [m] sigma [m] - Usage pattern: + The including heatSourceDict must define depth, innerA, innerB, outerA, + and outerB before including this file. For example: - #include "heatSourceDict.cfg" + depth 5.0e-5; + innerA 0.0; + innerB 1.0; + outerA 0.0; + outerB 1.0; - - { - ... + #include "$ADDITIVEFOAM_ETC/heatSources/nLightAFX-1000.cfg" - heatSourceModel nLightAFX; - - nLightAFXCoeffs - { - $Index0 - - transient true; - nPoints (10 10 10); - } - } + Select a mode inside nLightAFXCoeffs with $Index0 through $Index6. \*---------------------------------------------------------------------------*/ -// Initial heat-source depth for inner and outer components. -depth 5.0e-5; - -// Axial projected-Gaussian closure coefficients: n = A*log2(x) + B -innerA 0.0; -innerB 1.0; -outerA 0.0; -outerB 1.0; - // 0: D4 = 109.69 um Index0 { diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 00000000..5a803561 --- /dev/null +++ b/requirements.txt @@ -0,0 +1,10 @@ +arviz>=0.20,<0.23 +matplotlib +numpy +pandas +pymc +pytensor +PyYAML +reportlab +scipy +seaborn diff --git a/tutorials/heatSourceCalibration/README.md b/tutorials/heatSourceCalibration/README.md new file mode 100644 index 00000000..28b13ffc --- /dev/null +++ b/tutorials/heatSourceCalibration/README.md @@ -0,0 +1,106 @@ +# Heat-source calibration tutorial + +This tutorial calibrates the projected depth-distribution closure used by +AdditiveFOAM heat sources. The supplied worked example uses a +`projectedGaussian` source with SS316L. The laser spot size (D4sigma) was +measured to be 109.69 microns. + +## Installation and setup + +Build AdditiveFOAM against OpenFOAM 14, then source both environments: + +```bash +source /path/to/OpenFOAM-14/etc/bashrc +source /path/to/AdditiveFOAM/etc/bashrc +``` + +Copy the tutorial before running it: + +```bash +mkdir -p "$FOAM_RUN/AdditiveFOAM" +cp -r "$ADDITIVEFOAM_TUTORIALS/heatSourceCalibration" \ + "$FOAM_RUN/AdditiveFOAM/heatSourceCalibration" +cd "$FOAM_RUN/AdditiveFOAM/heatSourceCalibration" +``` + +The repository installation creates the Python environment at +`$ADDITIVEFOAM_PROJECT_DIR/.venv`. It is activated by the AdditiveFOAM +environment: + +```bash +calibrateHeatSource --help +``` + +Do not run a calibration inside `$ADDITIVEFOAM_TUTORIALS`; the campaign writes +cases and reports beneath the tutorial directory. + +## Run the calibration + +Review `system/decomposeParDict` in the template before starting. The supplied +research configuration requests 6 MPI ranks for each CFD case, evaluates ten +trial values for each of five experiments, and uses 2,000 posterior draws. + +```bash +calibrateHeatSource --config config.yml +``` + +The configuration resolves relative paths from the location of `config.yml`. +Environment variables and `~` are supported in the three `paths` entries. + +Generated output is contained in `campaign/`: + +```text +campaign/ +├── cases/ +│ └── P187p5_V500_D109p69/ +│ ├── B0/ +│ ├── B4p5/ +│ └── B9/ +├── simulations.yml +├── calibration_state.yml +├── calibration_fit.yml +└── reports/ + ├── calibration_report.pdf + └── calibration_summary.csv +``` + +Successful cases retain their rendered inputs, solver log, post-processing +output, and `metrics.yml`. With `keep_successful: false`, processor and numeric +time directories are removed after their melt-pool dimensions are recorded. + +## Material-derived liquidus + +The case includes: + +```foam +#include "$ADDITIVEFOAM_ETC/materials/SS316L.cfg" +``` + +The configuration selects `melt_pool_isovalue: liquidus`. The calibration +command expands `constant/transportProperties` with `foamDictionary`, obtains +the temperature paired with `alpha.solid = 0`, and reads the matching file from +`postProcessing/meltPoolDimensions`. The SS316L liquidus is therefore not +duplicated in `config.yml`, `heatSourceDict`, or `controlDict`. + +## Shared projected-source parameter + +The template exposes the trial value once in `constant/heatSourceDict`: + +```foam +calibrationA 0.0; +calibrationB <>; +``` + +The template references these aliases directly from the projected Gaussian +coefficients: + +```foam +projectedGaussianCoeffs +{ + dimensions (54.845e-6 54.845e-6 5.0e-5); + A $calibrationA; + B $calibrationB; +} +``` + +Power, speed, end time, and write interval are rendered independently. diff --git a/tutorials/heatSourceCalibration/config.yml b/tutorials/heatSourceCalibration/config.yml new file mode 100644 index 00000000..615435af --- /dev/null +++ b/tutorials/heatSourceCalibration/config.yml @@ -0,0 +1,33 @@ +paths: + experiments: experiments.yml + template: template + campaign: campaign + +case: + command: ./Allrun + keep_successful: false + melt_pool_isovalue: liquidus + +calibration: + token: B + initial_values: [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0] + bounds: [0.0, 9.0] + max_simulations_per_experiment: 10 + depth_tolerance_microns: 1.0 + posterior_std_tolerance: 0.15 + pchip_grid_points: 1000 + +bayesian: + draws: 2000 + tune: 1000 + cores: 8 + target_accept: 0.9 + random_seed: 42 + progressbar: true + +final_fit: + bootstrap_samples: 1000 + bootstrap_random_seed: 12345 + +report: + enabled: true diff --git a/tutorials/heatSourceCalibration/experiments.yml b/tutorials/heatSourceCalibration/experiments.yml new file mode 100644 index 00000000..479003fe --- /dev/null +++ b/tutorials/heatSourceCalibration/experiments.yml @@ -0,0 +1,25 @@ +- parameters: + Power_W: 187.5 + Speed_mm_s: 500 + Spot_Size_microns: 109.69 + Measured_Depth_microns: [87.05,98.1] +- parameters: + Power_W: 300.0 + Speed_mm_s: 500 + Spot_Size_microns: 109.69 + Measured_Depth_microns: [168.57,190.67] +- parameters: + Power_W: 412.5 + Speed_mm_s: 500 + Spot_Size_microns: 109.69 + Measured_Depth_microns: [272.19,264.07] +- parameters: + Power_W: 525.0 + Speed_mm_s: 500 + Spot_Size_microns: 109.69 + Measured_Depth_microns: [359.24,353.71] +- parameters: + Power_W: 637.5 + Speed_mm_s: 500 + Spot_Size_microns: 109.69 + Measured_Depth_microns: [464.25,467.01] diff --git a/tutorials/heatSourceCalibration/template/0/T b/tutorials/heatSourceCalibration/template/0/T new file mode 100644 index 00000000..66221b3f --- /dev/null +++ b/tutorials/heatSourceCalibration/template/0/T @@ -0,0 +1,50 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 14 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + object T; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 0 0 1 0 0 0]; + +internalField uniform 300; + +boundaryField +{ + bottom + { + type zeroGradient; + value uniform 300; + } + + top + { + type mixedTemperature; + h 10.0; + Tinf uniform 300; + value uniform 300; + } + + sides + { + type zeroGradient; + value uniform 300; + } + + internalFaces + { + type internal; + } +} + + +// ************************************************************************* // diff --git a/tutorials/heatSourceCalibration/template/0/U b/tutorials/heatSourceCalibration/template/0/U new file mode 100644 index 00000000..626f1188 --- /dev/null +++ b/tutorials/heatSourceCalibration/template/0/U @@ -0,0 +1,45 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 14 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volVectorField; + object U; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 1 -1 0 0 0 0]; + +internalField uniform (0 0 0); + +boundaryField +{ + bottom + { + type noSlip; + } + + top + { + type marangoni; + value uniform (0 0 0); + } + + sides + { + type noSlip; + } + + internalFaces + { + type internal; + } +} + +// ************************************************************************* // diff --git a/tutorials/heatSourceCalibration/template/0/p_rgh b/tutorials/heatSourceCalibration/template/0/p_rgh new file mode 100644 index 00000000..89657bfd --- /dev/null +++ b/tutorials/heatSourceCalibration/template/0/p_rgh @@ -0,0 +1,50 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 14 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + object p_rgh; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 2 -2 0 0 0 0]; + +internalField uniform 0; + +boundaryField +{ + bottom + { + type fixedFluxPressure; + rho rhok; + value uniform 0; + } + + top + { + type fixedFluxPressure; + rho rhok; + value uniform 0; + } + + sides + { + type fixedFluxPressure; + rho rhok; + value uniform 0; + } + + internalFaces + { + type internal; + } +} + +// ************************************************************************* // diff --git a/tutorials/heatSourceCalibration/template/Allclean b/tutorials/heatSourceCalibration/template/Allclean new file mode 100755 index 00000000..0dc3f1c6 --- /dev/null +++ b/tutorials/heatSourceCalibration/template/Allclean @@ -0,0 +1,7 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # Run from this directory + +# Source tutorial clean functions +. $WM_PROJECT_DIR/bin/tools/CleanFunctions + +cleanCase diff --git a/tutorials/heatSourceCalibration/template/Allrun b/tutorials/heatSourceCalibration/template/Allrun new file mode 100755 index 00000000..d4e57173 --- /dev/null +++ b/tutorials/heatSourceCalibration/template/Allrun @@ -0,0 +1,13 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # Run from this directory + +# Source tutorial run functions +. $WM_PROJECT_DIR/bin/tools/RunFunctions +application=`getApplication` + +runApplication blockMesh + +runApplication decomposePar + +runParallel $application + diff --git a/tutorials/heatSourceCalibration/template/README.md b/tutorials/heatSourceCalibration/template/README.md new file mode 100644 index 00000000..50194b84 --- /dev/null +++ b/tutorials/heatSourceCalibration/template/README.md @@ -0,0 +1,29 @@ +# Calibration case template + +This is the OpenFOAM-14 AdditiveFOAM case copied for every trial value in the +heat-source calibration campaign. It uses a `projectedGaussian` heat source +with SS316L. The laser spot size (D4sigma) was measured to be 109.69 microns. + +The calibration command expects this full case structure: + +```text +0/ +constant/ +system/ +Allrun +Allclean +``` + +The required renderer placeholders are: + +- `constant/heatSourceDict`: `<>`, exposed once as `calibrationB` +- `constant/scanPath`: `<>`, `<>` +- `system/controlDict`: `<>`, `<>` + +OpenFOAM dictionary aliases propagate `$calibrationB` into the projected +Gaussian closure. A template for another projected source can use the same +alias for every applicable `B` coefficient. + +The SS316L material configuration supplies `thermoPath`, emissivity, and the +Marangoni coefficient. Transient source depth and `meltPoolDimensions` obtain +their default isovalues from that material configuration. diff --git a/tutorials/heatSourceCalibration/template/constant/dynamicMeshDict b/tutorials/heatSourceCalibration/template/constant/dynamicMeshDict new file mode 100644 index 00000000..6f2c46cc --- /dev/null +++ b/tutorials/heatSourceCalibration/template/constant/dynamicMeshDict @@ -0,0 +1,39 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 14 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + format ascii; + class dictionary; + object dynamicMeshDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +topoChanger +{ + type refiner; + + libs ("libfvMeshTopoChangers.so"); + + refineInterval 1; + + field refinementField; + + lowerRefineLevel 0.1; + + upperRefineLevel 10.0; + + nBufferLayers 5; + + maxRefinement 1; + + maxCells 20000000; + + dumpLevel true; +} + +// ************************************************************************* // diff --git a/tutorials/heatSourceCalibration/template/constant/g b/tutorials/heatSourceCalibration/template/constant/g new file mode 100644 index 00000000..082543df --- /dev/null +++ b/tutorials/heatSourceCalibration/template/constant/g @@ -0,0 +1,22 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 14 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class uniformDimensionedVectorField; + location "constant"; + object g; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 1 -2 0 0 0 0]; +value (0 0 -9.81); + + +// ************************************************************************* // diff --git a/tutorials/heatSourceCalibration/template/constant/heatSourceDict b/tutorials/heatSourceCalibration/template/constant/heatSourceDict new file mode 100644 index 00000000..067d0733 --- /dev/null +++ b/tutorials/heatSourceCalibration/template/constant/heatSourceDict @@ -0,0 +1,66 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 14 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object heatSourceDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +sources (beam); + +beam +{ + pathName scanPath; + + absorptionModel Kelly; + + KellyCoeffs + { + eta0 0.27; + etaMin 0.27; + aspectRatioSwitch 0.0; + geometry cylinder; + } + + heatSourceModel projectedGaussian; + + projectedGaussianCoeffs + { + dimensions (54.845e-6 54.845e-6 30.0e-6); + A 0.0; + B <>; + + transient true; + nPoints (10 10 10); + } +} + +refinementModel +{ + refinementModel targetCellLoad; + + refinementTemperature 1000; + + buffers + { + beam (100.0e-6 100.0e-6 600e-6); + } + + targetCellLoadCoeffs + { + targetCellsPerProc 5000; + nBufferVolumes 4; + maxSearchIter 10; + timeTolerance 1e-4; + } +} + +// ************************************************************************* // diff --git a/tutorials/heatSourceCalibration/template/constant/scanPath b/tutorials/heatSourceCalibration/template/constant/scanPath new file mode 100644 index 00000000..a3da99e9 --- /dev/null +++ b/tutorials/heatSourceCalibration/template/constant/scanPath @@ -0,0 +1,3 @@ +Mode X Y Z Power Param +1 0.000 0.000 0 0 0 +0 0.002 0.000 0 <> <> diff --git a/tutorials/heatSourceCalibration/template/constant/transportProperties b/tutorials/heatSourceCalibration/template/constant/transportProperties new file mode 100644 index 00000000..88c76a6c --- /dev/null +++ b/tutorials/heatSourceCalibration/template/constant/transportProperties @@ -0,0 +1,18 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 14 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object transportProperties; +} + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +#include "$ADDITIVEFOAM_ETC/materials/SS316L.cfg" +// ************************************************************************* // diff --git a/tutorials/heatSourceCalibration/template/system/blockMeshDict b/tutorials/heatSourceCalibration/template/system/blockMeshDict new file mode 100644 index 00000000..19d4bb13 --- /dev/null +++ b/tutorials/heatSourceCalibration/template/system/blockMeshDict @@ -0,0 +1,85 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 14 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object blockMeshDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +xmin -0.0005; +xmax 0.0025; +ymin -0.0005; +ymax 0.0005; +zmin -0.0008; +zmax 0.0; + + +vertices +( + ($xmin $ymin $zmin) //0 + ($xmax $ymin $zmin) //1 + ($xmax $ymax $zmin) //2 + ($xmin $ymax $zmin) //3 + ($xmin $ymin $zmax) //4 + ($xmax $ymin $zmax) //5 + ($xmax $ymax $zmax) //6 + ($xmin $ymax $zmax) //7 +); + +blocks +( + hex (0 1 2 3 4 5 6 7) (75 25 40) simpleGrading (1 1 1) +); + +edges +( +); + +boundary +( + bottom + { + type wall; + faces + ( + (0 3 2 1) + ); + } + top + { + type wall; + faces + ( + (4 5 6 7) + ); + } + sides + { + type wall; + faces + ( + (0 4 7 3) + (2 6 5 1) + (1 5 4 0) + (3 7 6 2) + ); + } + internalFaces + { + type internal; + faces (); + } +); + +mergePatchPairs +( +); + +// ************************************************************************* // diff --git a/tutorials/heatSourceCalibration/template/system/controlDict b/tutorials/heatSourceCalibration/template/system/controlDict new file mode 100644 index 00000000..4a7607db --- /dev/null +++ b/tutorials/heatSourceCalibration/template/system/controlDict @@ -0,0 +1,69 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 14 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2; + format ascii; + class dictionary; + location "system"; + object controlDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +application additiveFoam; + +startFrom startTime; + +startTime 0; + +stopAt endTime; + +endTime <>; + +deltaT 1e-07; + +writeControl adjustableRunTime; + +writeInterval <>; + +purgeWrite 0; + +writeFormat binary; + +writePrecision 8; + +writeCompression off; + +timeFormat general; + +timePrecision 8; + +runTimeModifiable yes; + +adjustTimeStep yes; + +maxCo 0.5; + +maxDi 1; + +maxAlphaCo 1; + +functions +{ + meltPoolDimensions + { + libs ("libadditiveFoamFunctionObjects.so"); + + type meltPoolDimensions; + enabled true; + + scanPathAngle 0.0; + } +} + +// ************************************************************************* // diff --git a/tutorials/heatSourceCalibration/template/system/decomposeParDict b/tutorials/heatSourceCalibration/template/system/decomposeParDict new file mode 100644 index 00000000..50c8da78 --- /dev/null +++ b/tutorials/heatSourceCalibration/template/system/decomposeParDict @@ -0,0 +1,27 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 14 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object decomposeParDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +numberOfSubdomains 32; + +method simple; + +simple +{ + n (1 8 4); +} + +// ************************************************************************* // diff --git a/tutorials/heatSourceCalibration/template/system/fvSchemes b/tutorials/heatSourceCalibration/template/system/fvSchemes new file mode 100644 index 00000000..24f14a00 --- /dev/null +++ b/tutorials/heatSourceCalibration/template/system/fvSchemes @@ -0,0 +1,50 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 14 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object fvSchemes; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +ddtSchemes +{ + default Euler; +} + +gradSchemes +{ + default Gauss linear; +} + +divSchemes +{ + default Gauss upwind; +} + +laplacianSchemes +{ + default Gauss linear corrected; + laplacian(kappa,T) Gauss harmonic corrected; +} + +interpolationSchemes +{ + default linear; +} + +snGradSchemes +{ + default corrected; +} + + +// ************************************************************************* // diff --git a/tutorials/heatSourceCalibration/template/system/fvSolution b/tutorials/heatSourceCalibration/template/system/fvSolution new file mode 100644 index 00000000..d1f0ea8b --- /dev/null +++ b/tutorials/heatSourceCalibration/template/system/fvSolution @@ -0,0 +1,76 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 14 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object fvSolution; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +solvers +{ + "pcorr.*" + { + solver PCG; + preconditioner DIC; + tolerance 1e-2; + relTol 0; + } + + p_rgh + { + solver GAMG; + tolerance 1e-06; + relTol 0.01; + smoother DIC; + } + + p_rghFinal + { + $p_rgh; + relTol 0; + } + + "T.*" + { + solver PBiCGStab; + preconditioner DILU; + tolerance 1e-15; + relTol 0; + minIter 1; + maxIter 20; + } +} + +PIMPLE +{ + momentumPredictor no; + nOuterCorrectors 0; + nCorrectors 1; + nNonOrthogonalCorrectors 0; + pRefCell 0; + pRefValue 0; + + + nThermoCorrectors 20; + thermoTolerance 1e-8; + explicitSolve true; +} + +relaxationFactors +{ + equations + { + ".*" 1; + } +} + +// ************************************************************************* // diff --git a/tutorials/nLightAFX/README.md b/tutorials/nLightAFX/README.md index 3ac05dab..f272c294 100644 --- a/tutorials/nLightAFX/README.md +++ b/tutorials/nLightAFX/README.md @@ -27,7 +27,7 @@ Use `./Allclean` to remove generated mesh, decomposition, and result files. The important files for this tutorial are: ```text -constant/nLightAFX.cfg +$ADDITIVEFOAM_ETC/heatSources/nLightAFX-1000.cfg ``` Defines the ORNL-characterized AFX mode parameters. @@ -55,7 +55,13 @@ heatSourceModel nLightAFX; The corresponding coefficient dictionary is: ```foam -#include "nLightAFX.cfg" +depth 5.0e-5; +innerA 0.0; +innerB 1.0; +outerA 0.0; +outerB 1.0; + +#include "$ADDITIVEFOAM_ETC/heatSources/nLightAFX-1000.cfg" nLightAFXCoeffs { @@ -86,7 +92,8 @@ $Index6; ## Characterized AFX modes -The included `nLightAFX.cfg` file contains the ORNL-characterized AFX beam parameters for modes 0 through 6. Each mode defines: +The shared `nLightAFX-1000.cfg` file contains the ORNL-characterized AFX beam +parameters for modes 0 through 6. Each mode defines: ```foam dimensions @@ -133,7 +140,7 @@ where `x` is the ratio between the current heat source depth and lateral heat so ## Example mode -A typical mode block from `nLightAFX.cfg` looks like: +A typical mode block from `nLightAFX-1000.cfg` looks like: ```foam Index3 diff --git a/tutorials/nLightAFX/constant/heatSourceDict b/tutorials/nLightAFX/constant/heatSourceDict index 03250667..9b5c50f8 100644 --- a/tutorials/nLightAFX/constant/heatSourceDict +++ b/tutorials/nLightAFX/constant/heatSourceDict @@ -17,7 +17,13 @@ FoamFile sources (beam); -#include "nLightAFX.cfg" +depth 5.0e-5; +innerA 0.0; +innerB 1.0; +outerA 0.0; +outerB 1.0; + +#include "$ADDITIVEFOAM_ETC/heatSources/nLightAFX-1000.cfg" beam { From 5d521eb5340305c01552e4fe3a0330b556e4aed2 Mon Sep 17 00:00:00 2001 From: John Date: Mon, 20 Jul 2026 19:19:42 -0400 Subject: [PATCH 2/5] update build environment and decomposition --- Allwmake | 54 ------------------- README.md | 25 ++++++--- etc/bashrc | 20 ------- requirements.txt | 20 +++---- .../AMB2018-02-B/system/decomposeParDict | 2 +- tutorials/heatSourceCalibration/README.md | 12 +++-- .../template/system/decomposeParDict | 9 +--- tutorials/multiBeam/system/decomposeParDict | 2 +- .../multiLayerPBF/system/decomposeParDict | 2 +- tutorials/nLightAFX/system/decomposeParDict | 2 +- tutorials/tabulated/system/decomposeParDict | 2 +- 11 files changed, 41 insertions(+), 109 deletions(-) diff --git a/Allwmake b/Allwmake index 591905b9..830f8594 100755 --- a/Allwmake +++ b/Allwmake @@ -27,57 +27,3 @@ $ADDITIVEFOAM_BUILD_FLAGS \ applications/Allwmake $targetType $* #------------------------------------------------------------------------------ -# Build the Python environment - -pythonCommand= - -if command -v python3 >/dev/null 2>&1 -then - pythonCommand=python3 -elif command -v python >/dev/null 2>&1 && \ - python -c 'import sys; raise SystemExit(sys.version_info.major != 3)' -then - pythonCommand=python -fi - -if [ -z "$pythonCommand" ] -then - echo "WARNING: Python 3 was not found; skipping the Python environment." >&2 -elif ! "$pythonCommand" -c \ - 'import sys; raise SystemExit(sys.version_info < (3, 10))' -then - echo "WARNING: Python 3.10 or newer is required; skipping the Python environment." >&2 -elif [ ! -f "$ADDITIVEFOAM_PROJECT_DIR/requirements.txt" ] -then - echo "WARNING: requirements.txt was not found; skipping the Python environment." >&2 -else - additivefoamVenv="$ADDITIVEFOAM_PROJECT_DIR/.venv" - - if [ ! -x "$additivefoamVenv/bin/python" ] || \ - [ ! -r "$additivefoamVenv/bin/activate" ] - then - echo "Creating the AdditiveFOAM Python environment" - "$pythonCommand" -m venv "$additivefoamVenv" || exit 1 - fi - - if ! "$additivefoamVenv/bin/python" -c \ - 'import sys; raise SystemExit(sys.version_info < (3, 10))' - then - echo "ERROR: $additivefoamVenv uses Python older than 3.10." >&2 - echo " Remove .venv and run Allwmake again." >&2 - exit 1 - fi - - echo "Installing AdditiveFOAM Python dependencies" - "$additivefoamVenv/bin/python" -m pip install \ - -r "$ADDITIVEFOAM_PROJECT_DIR/requirements.txt" || exit 1 - "$additivefoamVenv/bin/python" -m pip check || exit 1 - - echo "Source $ADDITIVEFOAM_PROJECT_DIR/etc/bashrc to activate the environment" - - unset additivefoamVenv -fi - -unset pythonCommand - -#------------------------------------------------------------------------------ diff --git a/README.md b/README.md index af14a54d..0b972080 100644 --- a/README.md +++ b/README.md @@ -42,21 +42,30 @@ script from the repository root: ```sh ./Allwmake -source etc/bashrc ``` -The second `source` activates the newly created Python environment. The master -build script uses Python 3.10 or newer to create `.venv` and install -the dependencies in `requirements.txt`. If Python is unavailable, the compiled -applications are still built and `Allwmake` prints a warning. +The Python utilities require Python 3.10 or newer. To use a dedicated virtual +environment, create it and install the pinned dependencies from the repository +root: ```sh +python3 -m venv .venv +. .venv/bin/activate +python -m pip install -r requirements.txt +python -m pip check calibrateHeatSource --help ``` -For regular use, source both environments in each new shell or add them to your -shell startup file. The AdditiveFOAM environment activates `.venv` when it is -available: +To use an existing Python environment instead, activate that environment and +install the same requirements into it: + +```sh +python -m pip install -r /path/to/AdditiveFOAM/requirements.txt +python -m pip check +``` + +The AdditiveFOAM environment does not create or activate a Python environment. +Source the OpenFOAM and AdditiveFOAM environments in each new shell: ```sh source /path/to/OpenFOAM-14/etc/bashrc diff --git a/etc/bashrc b/etc/bashrc index 18c73084..106e370e 100644 --- a/etc/bashrc +++ b/etc/bashrc @@ -77,26 +77,6 @@ if [ -d "$ADDITIVEFOAM_PROJECT_DIR/bin" ]; then fi -# AdditiveFOAM Python environment -additivefoamVenv="$ADDITIVEFOAM_PROJECT_DIR/.venv" - -if [ -d "$additivefoamVenv" ] -then - if [ ! -x "$additivefoamVenv/bin/python" ] || \ - [ ! -r "$additivefoamVenv/bin/activate" ] - then - echo "WARNING: The AdditiveFOAM Python environment is incomplete." >&2 - echo " Run $ADDITIVEFOAM_PROJECT_DIR/Allwmake to rebuild it." >&2 - elif [ "${VIRTUAL_ENV:-}" != "$additivefoamVenv" ] - then - . "$additivefoamVenv/bin/activate" || \ - echo "WARNING: Unable to activate $additivefoamVenv." >&2 - fi -fi - -unset additivefoamVenv - - # AdditiveFOAM directories export ADDITIVEFOAM_TUTORIALS="$ADDITIVEFOAM_PROJECT_DIR/tutorials" export ADDITIVEFOAM_APPLICATIONS="$ADDITIVEFOAM_PROJECT_DIR/applications" diff --git a/requirements.txt b/requirements.txt index 5a803561..65e7397a 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,10 +1,10 @@ -arviz>=0.20,<0.23 -matplotlib -numpy -pandas -pymc -pytensor -PyYAML -reportlab -scipy -seaborn +arviz==0.22.0 +matplotlib==3.10.8 +numpy==1.26.4 +pandas==2.3.3 +pymc==5.25.1 +pytensor==2.31.7 +PyYAML==6.0.2 +reportlab==5.0.0 +scipy==1.15.3 +seaborn==0.13.2 diff --git a/tutorials/AMB2018-02-B/system/decomposeParDict b/tutorials/AMB2018-02-B/system/decomposeParDict index 8d6b9f48..da6a8474 100644 --- a/tutorials/AMB2018-02-B/system/decomposeParDict +++ b/tutorials/AMB2018-02-B/system/decomposeParDict @@ -15,7 +15,7 @@ FoamFile } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -numberOfSubdomains 6; +numberOfSubdomains 8; method scotch; diff --git a/tutorials/heatSourceCalibration/README.md b/tutorials/heatSourceCalibration/README.md index 28b13ffc..0e0ab753 100644 --- a/tutorials/heatSourceCalibration/README.md +++ b/tutorials/heatSourceCalibration/README.md @@ -23,11 +23,13 @@ cp -r "$ADDITIVEFOAM_TUTORIALS/heatSourceCalibration" \ cd "$FOAM_RUN/AdditiveFOAM/heatSourceCalibration" ``` -The repository installation creates the Python environment at -`$ADDITIVEFOAM_PROJECT_DIR/.venv`. It is activated by the AdditiveFOAM -environment: +Install the Python dependencies into the environment you intend to use before +running the calibration: ```bash +python -m pip install \ + -r "$ADDITIVEFOAM_PROJECT_DIR/requirements.txt" +python -m pip check calibrateHeatSource --help ``` @@ -37,8 +39,8 @@ cases and reports beneath the tutorial directory. ## Run the calibration Review `system/decomposeParDict` in the template before starting. The supplied -research configuration requests 6 MPI ranks for each CFD case, evaluates ten -trial values for each of five experiments, and uses 2,000 posterior draws. +configuration uses 8 MPI ranks for each AdditiveFOAM case, evaluates ten trial +values for each of five experiments, and uses 2,000 posterior draws. ```bash calibrateHeatSource --config config.yml diff --git a/tutorials/heatSourceCalibration/template/system/decomposeParDict b/tutorials/heatSourceCalibration/template/system/decomposeParDict index 50c8da78..da6a8474 100644 --- a/tutorials/heatSourceCalibration/template/system/decomposeParDict +++ b/tutorials/heatSourceCalibration/template/system/decomposeParDict @@ -15,13 +15,8 @@ FoamFile } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -numberOfSubdomains 32; +numberOfSubdomains 8; -method simple; - -simple -{ - n (1 8 4); -} +method scotch; // ************************************************************************* // diff --git a/tutorials/multiBeam/system/decomposeParDict b/tutorials/multiBeam/system/decomposeParDict index 8d6b9f48..da6a8474 100644 --- a/tutorials/multiBeam/system/decomposeParDict +++ b/tutorials/multiBeam/system/decomposeParDict @@ -15,7 +15,7 @@ FoamFile } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -numberOfSubdomains 6; +numberOfSubdomains 8; method scotch; diff --git a/tutorials/multiLayerPBF/system/decomposeParDict b/tutorials/multiLayerPBF/system/decomposeParDict index 8d6b9f48..da6a8474 100644 --- a/tutorials/multiLayerPBF/system/decomposeParDict +++ b/tutorials/multiLayerPBF/system/decomposeParDict @@ -15,7 +15,7 @@ FoamFile } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -numberOfSubdomains 6; +numberOfSubdomains 8; method scotch; diff --git a/tutorials/nLightAFX/system/decomposeParDict b/tutorials/nLightAFX/system/decomposeParDict index 8d6b9f48..da6a8474 100644 --- a/tutorials/nLightAFX/system/decomposeParDict +++ b/tutorials/nLightAFX/system/decomposeParDict @@ -15,7 +15,7 @@ FoamFile } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -numberOfSubdomains 6; +numberOfSubdomains 8; method scotch; diff --git a/tutorials/tabulated/system/decomposeParDict b/tutorials/tabulated/system/decomposeParDict index 8d6b9f48..da6a8474 100644 --- a/tutorials/tabulated/system/decomposeParDict +++ b/tutorials/tabulated/system/decomposeParDict @@ -15,7 +15,7 @@ FoamFile } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -numberOfSubdomains 6; +numberOfSubdomains 8; method scotch; From bd15c959bab31a7c1fc3b7f4a169e99fbec00408 Mon Sep 17 00:00:00 2001 From: John Date: Mon, 20 Jul 2026 19:37:07 -0400 Subject: [PATCH 3/5] fix up CI --- .github/workflows/CI.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 85a97343..cf98addc 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -42,4 +42,4 @@ jobs: cd userCase blockMesh decomposePar - mpirun -n 6 --oversubscribe --allow-run-as-root additiveFoam -parallel + mpirun -n 8 --oversubscribe --allow-run-as-root additiveFoam -parallel From 2455d9fbb348ac6a85881fe5bfb2f4ae08238c3f Mon Sep 17 00:00:00 2001 From: John Date: Mon, 20 Jul 2026 19:50:36 -0400 Subject: [PATCH 4/5] fixup --- bin/calibrateHeatSource | 27 +++++++++++++++++-- tutorials/heatSourceCalibration/README.md | 22 ++++++++------- .../heatSourceCalibration/template/README.md | 12 ++++++--- .../template/constant/heatSourceDict | 4 ++- 4 files changed, 48 insertions(+), 17 deletions(-) diff --git a/bin/calibrateHeatSource b/bin/calibrateHeatSource index 17a61437..00004d30 100755 --- a/bin/calibrateHeatSource +++ b/bin/calibrateHeatSource @@ -64,6 +64,22 @@ def _create_data_fingerprint(data_list: list) -> str: return hashlib.sha256(sorted_data_str.encode()).hexdigest() +def spot_size_2sigma_metres(params: dict) -> float: + try: + spot_size_d4sigma_microns = float(params["Spot_Size_microns"]) + except (KeyError, TypeError, ValueError) as error: + raise ValueError( + "Each experiment requires a numeric Spot_Size_microns value" + ) from error + + if not math.isfinite(spot_size_d4sigma_microns): + raise ValueError("Spot_Size_microns must be finite") + if spot_size_d4sigma_microns <= 0.0: + raise ValueError("Spot_Size_microns must be greater than zero") + + return 0.5e-6 * spot_size_d4sigma_microns + + def load_yaml_file(filepath: str): if not os.path.exists(filepath): return [] @@ -209,7 +225,8 @@ def parse_experiments(raw_experiments: list) -> pd.DataFrame: records = [] for exp in raw_experiments: params = exp["parameters"] - spot_size = params.get("Spot_Size_microns") + spot_size_2sigma_metres(params) + spot_size = float(params["Spot_Size_microns"]) depths = exp["Measured_Depth_microns"] # The projected-source closure uses depth relative to lateral spot size. normalized_depths = [float(d) / float(spot_size) for d in depths] @@ -355,7 +372,10 @@ def validate_template_case(template_case: Path, calibration_token: str): raise ValueError("Template case is missing: " + ", ".join(missing)) placeholder_locations = { - "constant/heatSourceDict": [calibration_token], + "constant/heatSourceDict": [ + calibration_token, + "spotSize2Sigma", + ], "constant/scanPath": ["power", "velocity"], "system/controlDict": ["endTime", "writeInterval"], } @@ -593,10 +613,12 @@ def run_one_additivefoam_simulation( shutil.copytree(template_case, case_dir) velocity_m_s = float(params["Speed_mm_s"]) * 1e-3 + spot_size_2sigma_m = spot_size_2sigma_metres(params) values = { config["calibration_token"]: float(B), "power": float(params["Power_W"]), "velocity": velocity_m_s, + "spotSize2Sigma": spot_size_2sigma_m, } # Render the scan first because its geometry and velocity determine endTime. @@ -632,6 +654,7 @@ def run_one_additivefoam_simulation( "B": float(B), "power": float(params["Power_W"]), "velocity": velocity_m_s, + "spot_size_2sigma_m": spot_size_2sigma_m, "endTime": end_time, "writeInterval": end_time, } diff --git a/tutorials/heatSourceCalibration/README.md b/tutorials/heatSourceCalibration/README.md index 0e0ab753..607310c2 100644 --- a/tutorials/heatSourceCalibration/README.md +++ b/tutorials/heatSourceCalibration/README.md @@ -84,25 +84,27 @@ the temperature paired with `alpha.solid = 0`, and reads the matching file from `postProcessing/meltPoolDimensions`. The SS316L liquidus is therefore not duplicated in `config.yml`, `heatSourceDict`, or `controlDict`. -## Shared projected-source parameter +## Rendered projected-source parameters -The template exposes the trial value once in `constant/heatSourceDict`: +The template exposes the lateral source dimension once in +`constant/heatSourceDict`: ```foam -calibrationA 0.0; -calibrationB <>; +spotSize2Sigma <>; ``` -The template references these aliases directly from the projected Gaussian -coefficients: +The projected Gaussian coefficients use that value for both lateral dimensions +and render the trial value directly: ```foam projectedGaussianCoeffs { - dimensions (54.845e-6 54.845e-6 5.0e-5); - A $calibrationA; - B $calibrationB; + dimensions ($spotSize2Sigma $spotSize2Sigma 30.0e-6); + A 0.0; + B <>; } ``` -Power, speed, end time, and write interval are rendered independently. +For each experiment, `spotSize2Sigma` is half of `Spot_Size_microns` and is +converted from microns to metres. Power, speed, end time, and write interval +are also rendered for each case. diff --git a/tutorials/heatSourceCalibration/template/README.md b/tutorials/heatSourceCalibration/template/README.md index 50194b84..f209c6df 100644 --- a/tutorials/heatSourceCalibration/template/README.md +++ b/tutorials/heatSourceCalibration/template/README.md @@ -16,13 +16,17 @@ Allclean The required renderer placeholders are: -- `constant/heatSourceDict`: `<>`, exposed once as `calibrationB` +- `constant/heatSourceDict`: `<>` and `<>` - `constant/scanPath`: `<>`, `<>` - `system/controlDict`: `<>`, `<>` -OpenFOAM dictionary aliases propagate `$calibrationB` into the projected -Gaussian closure. A template for another projected source can use the same -alias for every applicable `B` coefficient. +The projected Gaussian closure reads `<>` directly. The +`$spotSize2Sigma` alias supplies both lateral dimensions from one rendered +value. + +The calibration command converts `Spot_Size_microns` from the measured D4sigma +diameter to a 2sigma radius in metres and assigns it to both lateral source +dimensions. The SS316L material configuration supplies `thermoPath`, emissivity, and the Marangoni coefficient. Transient source depth and `meltPoolDimensions` obtain diff --git a/tutorials/heatSourceCalibration/template/constant/heatSourceDict b/tutorials/heatSourceCalibration/template/constant/heatSourceDict index 067d0733..2fa3a379 100644 --- a/tutorials/heatSourceCalibration/template/constant/heatSourceDict +++ b/tutorials/heatSourceCalibration/template/constant/heatSourceDict @@ -16,6 +16,8 @@ FoamFile sources (beam); +spotSize2Sigma <>; + beam { pathName scanPath; @@ -34,7 +36,7 @@ beam projectedGaussianCoeffs { - dimensions (54.845e-6 54.845e-6 30.0e-6); + dimensions ($spotSize2Sigma $spotSize2Sigma 30.0e-6); A 0.0; B <>; From 6da10d05d0a739e68085daf49794f3e220f69808 Mon Sep 17 00:00:00 2001 From: John Date: Mon, 20 Jul 2026 20:07:31 -0400 Subject: [PATCH 5/5] fixed up spot size --- bin/calibrateHeatSource | 48 +++++++++---------- tutorials/heatSourceCalibration/README.md | 18 ++++--- .../heatSourceCalibration/template/README.md | 13 ++--- .../template/constant/heatSourceDict | 4 +- 4 files changed, 44 insertions(+), 39 deletions(-) diff --git a/bin/calibrateHeatSource b/bin/calibrateHeatSource index 00004d30..9949c256 100755 --- a/bin/calibrateHeatSource +++ b/bin/calibrateHeatSource @@ -6,7 +6,7 @@ For each experimental condition, the routine builds a response curve from CFD results obtained with A = 0 and trial values of B. It infers a local n from the measured melt-pool depths, then fits the production relation - n = clip(A * log2(depth / spot_size) + B, 0, 9) + n = clip(A * log2(depth / heat_source_radius) + B, 0, 9) Local posterior uncertainty is used in both the global fit and its uncertainty band. Simulation results and calibration state are cached by process condition. @@ -64,7 +64,7 @@ def _create_data_fingerprint(data_list: list) -> str: return hashlib.sha256(sorted_data_str.encode()).hexdigest() -def spot_size_2sigma_metres(params: dict) -> float: +def heat_source_radius_microns(params: dict) -> float: try: spot_size_d4sigma_microns = float(params["Spot_Size_microns"]) except (KeyError, TypeError, ValueError) as error: @@ -77,7 +77,7 @@ def spot_size_2sigma_metres(params: dict) -> float: if spot_size_d4sigma_microns <= 0.0: raise ValueError("Spot_Size_microns must be greater than zero") - return 0.5e-6 * spot_size_d4sigma_microns + return 0.5 * spot_size_d4sigma_microns def load_yaml_file(filepath: str): @@ -225,11 +225,9 @@ def parse_experiments(raw_experiments: list) -> pd.DataFrame: records = [] for exp in raw_experiments: params = exp["parameters"] - spot_size_2sigma_metres(params) - spot_size = float(params["Spot_Size_microns"]) + heat_source_radius = heat_source_radius_microns(params) depths = exp["Measured_Depth_microns"] - # The projected-source closure uses depth relative to lateral spot size. - normalized_depths = [float(d) / float(spot_size) for d in depths] + normalized_depths = [float(d) / heat_source_radius for d in depths] records.append( { @@ -251,7 +249,7 @@ def parse_simulations(raw_simulations: list) -> pd.DataFrame: params = sim_run["parameters"] n_values = sim_run["n"] depths = sim_run["Simulated_Depth_microns"] - spot_size = params.get("Spot_Size_microns") + heat_source_radius = heat_source_radius_microns(params) for n, depth in zip(n_values, depths): records.append( @@ -260,7 +258,7 @@ def parse_simulations(raw_simulations: list) -> pd.DataFrame: "n": float(n), "Simulated_Depth_microns": float(depth), "Normalized_Simulated_Depth": float(depth) - / float(spot_size), + / heat_source_radius, } ) return pd.DataFrame(records) @@ -374,7 +372,7 @@ def validate_template_case(template_case: Path, calibration_token: str): placeholder_locations = { "constant/heatSourceDict": [ calibration_token, - "spotSize2Sigma", + "heatSourceRadius", ], "constant/scanPath": ["power", "velocity"], "system/controlDict": ["endTime", "writeInterval"], @@ -613,12 +611,12 @@ def run_one_additivefoam_simulation( shutil.copytree(template_case, case_dir) velocity_m_s = float(params["Speed_mm_s"]) * 1e-3 - spot_size_2sigma_m = spot_size_2sigma_metres(params) + heat_source_radius_m = heat_source_radius_microns(params) * 1.0e-6 values = { config["calibration_token"]: float(B), "power": float(params["Power_W"]), "velocity": velocity_m_s, - "spotSize2Sigma": spot_size_2sigma_m, + "heatSourceRadius": heat_source_radius_m, } # Render the scan first because its geometry and velocity determine endTime. @@ -654,7 +652,7 @@ def run_one_additivefoam_simulation( "B": float(B), "power": float(params["Power_W"]), "velocity": velocity_m_s, - "spot_size_2sigma_m": spot_size_2sigma_m, + "heat_source_radius_m": heat_source_radius_m, "endTime": end_time, "writeInterval": end_time, } @@ -1209,7 +1207,7 @@ def plot_calibration_overview_pages( ) ax.set_xlabel("B used in A=0 local calibration") - ax.set_ylabel("Depth / spot size") + ax.set_ylabel("d / (D4sigma / 2)") ax.margins(x=0.08, y=0.12) ax.legend( @@ -1572,8 +1570,8 @@ def plot_final_fit( """ Plot the fitted projected-source response and its uncertainty. - ``x_obs`` is measured depth divided by spot size. Values below one are - valid. + ``x_obs`` is measured depth divided by half the D4sigma spot size. Values + below one are valid. The curve and interval use n = clip(A * log2(x) + B, 0, 9). """ with report_plot_context(): @@ -1674,7 +1672,7 @@ def plot_final_fit( ax.set_xscale("log", base=2) ax.set_title("AdditiveFOAM Calibration Fit") - ax.set_xlabel("x = measured depth / spot size") + ax.set_xlabel("x = d / (D4sigma / 2)") ax.set_ylabel("Shape factor n") ax.margins(x=0.04, y=0.12) @@ -1702,8 +1700,9 @@ def fit_and_plot_heteroskedastic_model( n = clip(A * log2(x) + B, 0, 9) - where x is measured depth divided by spot size. The reported 95% response - interval is propagated from the local calibration uncertainty. + where x is measured depth divided by half the D4sigma spot size. The + reported 95% response interval is propagated from the local calibration + uncertainty. """ if calibrated_results_df.empty or len(calibrated_results_df) < 2: return None, None @@ -1833,6 +1832,7 @@ def fit_and_plot_heteroskedastic_model( "x_min": x_min, "x_max": x_max, "x_log_floor": float(x_log_floor), + "x_definition": "d / (D4sigma / 2)", "n_raw_at_x_min": float(n_raw_at_x_min), "n_raw_at_x_max": float(n_raw_at_x_max), "n_applied_at_x_min": float(np.clip(n_raw_at_x_min, 0.0, 9.0)), @@ -1971,7 +1971,7 @@ def save_summary_table(results_df: pd.DataFrame, filepath: Path): { **params, "mean_depth_microns": row["mean_depth_microns"], - "x_depth_over_spot_size": row["mean_x"], + "x_depth_over_half_d4sigma": row["mean_x"], "calibrated_n": row["calibrated_n"], "calibrated_n_mode": row.get( "calibrated_n_mode", row["calibrated_n"] @@ -2221,7 +2221,7 @@ def write_report( run_items.extend( [ ( - "Range x = depth / spot size", + "Range x = d / (D4sigma / 2)", f"[{results_df['mean_x'].min():.4f}, " f"{results_df['mean_x'].max():.4f}]", ), @@ -2303,7 +2303,7 @@ def write_report( header_map = { "mean_depth_microns": "Target mean, microns", - "mean_x": "depth / spot size", + "mean_x": "d / (D4sigma / 2)", "calibrated_n": "n estimate", "calibrated_n_95_low": "n, lower 95%", "calibrated_n_95_high": "n, upper 95%", @@ -2442,8 +2442,8 @@ def write_report( f"Figure {figure_number}: " "Final AdditiveFOAM calibration fit. Symbols show local n " "estimates with 95% intervals. The solid curve is n = A " - "log2(x) + B, where x is measured depth divided by spot " - "size. The shaded region is the 95% interval propagated " + "log2(x) + B, where x = d / (D4sigma / 2). The shaded " + "region is the 95% interval propagated " "from local calibration uncertainty.", styles["Caption"], ) diff --git a/tutorials/heatSourceCalibration/README.md b/tutorials/heatSourceCalibration/README.md index 607310c2..48eaaef3 100644 --- a/tutorials/heatSourceCalibration/README.md +++ b/tutorials/heatSourceCalibration/README.md @@ -2,8 +2,8 @@ This tutorial calibrates the projected depth-distribution closure used by AdditiveFOAM heat sources. The supplied worked example uses a -`projectedGaussian` source with SS316L. The laser spot size (D4sigma) was -measured to be 109.69 microns. +`projectedGaussian` source with SS316L. The laser D4sigma diameter was measured +to be 109.69 microns, giving a 2sigma heat-source radius of 54.845 microns. ## Installation and setup @@ -90,7 +90,7 @@ The template exposes the lateral source dimension once in `constant/heatSourceDict`: ```foam -spotSize2Sigma <>; +heatSourceRadius <>; ``` The projected Gaussian coefficients use that value for both lateral dimensions @@ -99,12 +99,16 @@ and render the trial value directly: ```foam projectedGaussianCoeffs { - dimensions ($spotSize2Sigma $spotSize2Sigma 30.0e-6); + dimensions ($heatSourceRadius $heatSourceRadius 30.0e-6); A 0.0; B <>; } ``` -For each experiment, `spotSize2Sigma` is half of `Spot_Size_microns` and is -converted from microns to metres. Power, speed, end time, and write interval -are also rendered for each case. +For each experiment, `Spot_Size_microns` is the measured D4sigma diameter. The +calibration command divides it by two to obtain the heat-source radius, then +uses that radius for both lateral source dimensions and to normalize measured +and simulated depth. The calibration coordinate is therefore +`d / (D4sigma / 2)`, matching `d / min(dimensions.x(), dimensions.y())` in the +heat-source models. Power, speed, end time, and write interval are also +rendered for each case. diff --git a/tutorials/heatSourceCalibration/template/README.md b/tutorials/heatSourceCalibration/template/README.md index f209c6df..26a907c4 100644 --- a/tutorials/heatSourceCalibration/template/README.md +++ b/tutorials/heatSourceCalibration/template/README.md @@ -2,7 +2,8 @@ This is the OpenFOAM-14 AdditiveFOAM case copied for every trial value in the heat-source calibration campaign. It uses a `projectedGaussian` heat source -with SS316L. The laser spot size (D4sigma) was measured to be 109.69 microns. +with SS316L. The laser D4sigma diameter was measured to be 109.69 microns, +giving a 2sigma heat-source radius of 54.845 microns. The calibration command expects this full case structure: @@ -16,17 +17,17 @@ Allclean The required renderer placeholders are: -- `constant/heatSourceDict`: `<>` and `<>` +- `constant/heatSourceDict`: `<>` and `<>` - `constant/scanPath`: `<>`, `<>` - `system/controlDict`: `<>`, `<>` The projected Gaussian closure reads `<>` directly. The -`$spotSize2Sigma` alias supplies both lateral dimensions from one rendered +`$heatSourceRadius` alias supplies both lateral dimensions from one rendered value. -The calibration command converts `Spot_Size_microns` from the measured D4sigma -diameter to a 2sigma radius in metres and assigns it to both lateral source -dimensions. +The calibration command divides `Spot_Size_microns` by two, converts the +result to metres, assigns it to both lateral source dimensions, and uses the +same radius to normalize measured and simulated depth. The SS316L material configuration supplies `thermoPath`, emissivity, and the Marangoni coefficient. Transient source depth and `meltPoolDimensions` obtain diff --git a/tutorials/heatSourceCalibration/template/constant/heatSourceDict b/tutorials/heatSourceCalibration/template/constant/heatSourceDict index 2fa3a379..8aa56f5c 100644 --- a/tutorials/heatSourceCalibration/template/constant/heatSourceDict +++ b/tutorials/heatSourceCalibration/template/constant/heatSourceDict @@ -16,7 +16,7 @@ FoamFile sources (beam); -spotSize2Sigma <>; +heatSourceRadius <>; beam { @@ -36,7 +36,7 @@ beam projectedGaussianCoeffs { - dimensions ($spotSize2Sigma $spotSize2Sigma 30.0e-6); + dimensions ($heatSourceRadius $heatSourceRadius 30.0e-6); A 0.0; B <>;