diff --git a/modules/optimization/examples/bimaterial/comparison.png b/modules/optimization/examples/bimaterial/comparison.png new file mode 100644 index 000000000000..67d37a9c8630 Binary files /dev/null and b/modules/optimization/examples/bimaterial/comparison.png differ diff --git a/modules/optimization/examples/bimaterial/gold/main_out_optInfo_0001_lmvm.csv b/modules/optimization/examples/bimaterial/gold/main_out_optInfo_0001_lmvm.csv new file mode 100644 index 000000000000..40e8a40ef31d --- /dev/null +++ b/modules/optimization/examples/bimaterial/gold/main_out_optInfo_0001_lmvm.csv @@ -0,0 +1,8 @@ +cnorm,current_iterate,function_value,gnorm,gradient_iterate,hessian_iterate,objective_iterate,xdiff +0,0,0.0027559190255043,0.0020980201625854,0,0,5,1 +0,1,0.0011436194086151,0.0012646984507389,0,0,25,341 +0,2,0.00066401427375986,0.002609319456353,0,0,5,1 +0,3,0.00030817664707322,0.0014179698042894,0,0,5,1 +0,4,2.0142789211171e-05,0.00036418217816026,0,0,5,1 +0,5,5.7454708065525e-06,7.8329160282155e-05,0,0,5,1 +0,6,4.8895578321809e-06,5.5213436174544e-06,0,0,5,1 diff --git a/modules/optimization/examples/bimaterial/gold/main_out_optInfo_0001_nm.csv b/modules/optimization/examples/bimaterial/gold/main_out_optInfo_0001_nm.csv new file mode 100644 index 000000000000..09b4375c8dd7 --- /dev/null +++ b/modules/optimization/examples/bimaterial/gold/main_out_optInfo_0001_nm.csv @@ -0,0 +1,16 @@ +cnorm,current_iterate,function_value,gnorm,gradient_iterate,hessian_iterate,objective_iterate,xdiff +0,1,0.0027559190255043,0.0020800280563548,0,0,3,1 +0,2,0.00089784067857055,0.0034672376867739,0,0,2,1 +0,3,0.00089784067857055,0.0018580783469338,0,0,1,1 +0,4,0.00089784067857055,0.00060014784528086,0,0,1,1 +0,5,0.00012032946950342,0.0012502070399019,0,0,2,1 +0,6,0.00012032946950342,0.00077751120906713,0,0,2,1 +0,7,0.00012032946950342,0.00018088323722629,0,0,2,1 +0,8,0.00010539609933997,6.1690449462141e-05,0,0,2,1 +0,9,2.8944854931607e-05,9.1384614571812e-05,0,0,2,1 +0,10,2.8944854931607e-05,7.6451244408368e-05,0,0,2,1 +0,11,5.7147045142231e-06,5.3978230616923e-05,0,0,2,1 +0,12,5.7147045142231e-06,2.3230150417384e-05,0,0,2,1 +0,13,5.7147045142231e-06,1.5221607202598e-05,0,0,1,1 +0,14,5.7147045142231e-06,1.3244963526617e-05,0,0,2,1 +0,15,5.7147045142231e-06,4.1115087138692e-06,0,0,2,1 diff --git a/modules/optimization/examples/bimaterial/gold/pest.rec b/modules/optimization/examples/bimaterial/gold/pest.rec new file mode 100644 index 000000000000..7f06a86521ba --- /dev/null +++ b/modules/optimization/examples/bimaterial/gold/pest.rec @@ -0,0 +1,380 @@ + PEST RUN RECORD: CASE pest + + + PEST Version: 17.3 + + + PEST run mode:- + + Parameter estimation mode + + + Case dimensions:- + + Number of parameters : 2 + Number of adjustable parameters : 2 + Number of parameter groups : 1 + Number of observations : 4 + Number of prior estimates : 0 + + + Model command line(s):- + + ./run_model.py + + + Jacobian command line:- + + na + + + Model interface files:- + + Templates: + pest_params.tpl + for model input files: + pest_params.csv + + (Parameter values written using single precision protocol.) + (Decimal point always included.) + + Instruction files: + pest_instruction.ins + for reading model output files: + pest_to_read.txt + + + PEST-to-model message file:- + + na + + + Derivatives calculation:- + + Param Increment Increment Increment Forward or Multiplier Method + group type low bound central (central) (central) + p relative 1.0000E-02 none switch 2.000 parabolic + + + Parameter definitions:- + + Name Trans- Change Initial Lower Upper + formation limit value bound bound + ppktop none factor 2.00000 0.100000 20.0000 + ppkbot none factor 2.00000 0.100000 20.0000 + + Name Group Scale Offset Model command number + ppktop p 1.00000 0.00000 1 + ppkbot p 1.00000 0.00000 1 + + + Prior information:- + + No prior information supplied + + + Observations:- + + Observation name Observation Weight Group + ttopmid 0.138000 1.000 obsgroup + tbotmid 4.000000E-02 1.000 obsgroup + tbotlef 2.200000E-02 1.000 obsgroup + tbotrig 2.200000E-02 1.000 obsgroup + + + Control settings:- + + Initial lambda : 10.000 + Lambda adjustment factor : 2.0000 + Sufficient new/old phi ratio per optimisation iteration : 0.30000 + Limiting relative phi reduction between lambdas : 3.00000E-02 + Maximum trial lambdas per iteration : 10 + Forgive model run failure during lamda testing : no + Forgive model run failure during Jacobian runs : no + + Perform Broyden's update of Jacobian matrix : no + Undertake observation re-referencing : no + + Maximum factor parameter change (factor-limited changes) : 5.0000 + Maximum relative parameter change (relative-limited changes) : na + Fraction of initial parameter values used in computing + change limit for near-zero parameters : 1.00000E-03 + Allow bending of parameter upgrade vector : no + Allow parameters to stick to their bounds : no + + Relative phi reduction below which to begin use of + central derivatives : 0.10000 + Iteration at which to first consider derivatives switch : 1 + + Relative phi reduction indicating convergence : 0.10000E-01 + Number of phi values required within this range : 4 + Maximum number of consecutive failures to lower phi : 3 + Minimal relative parameter change indicating convergence : 0.10000E-01 + Number of consecutive iterations with minimal param change : 3 + Maximum number of optimisation iterations : 100 + + Attempt automatic user intervention : no + + Attempt reuse of parameter sensitivities : no + + Scale parameters by their bounds : no + + + File saving options: - + + Save best JCO file : yes + Save multiple JCO files : no + Save multiple REI files : no + Save multiple PAR files : no + + + OPTIMISATION RECORD + + + INITIAL CONDITIONS: + Sum of squared weighted residuals (ie phi) = 5.62473E-03 + + Current parameter values + ppktop 2.00000 + ppkbot 2.00000 + + + OPTIMISATION ITERATION NO. : 1 + Model calls so far : 1 + Starting phi for this iteration: 5.62473E-03 + + Lambda = 10.000 -----> + Phi = 1.79112E-03 ( 0.318 of starting phi) + + Lambda = 5.0000 -----> + Phi = 1.68259E-03 ( 0.299 of starting phi) + + No more lambdas: phi is less than 0.3000 of starting phi + Lowest phi this iteration: 1.68259E-03 + + Current parameter values Previous parameter values + ppktop 3.49071 ppktop 2.00000 + ppkbot 0.400000 ppkbot 2.00000 + Maximum factor change: 5.000 ["ppkbot"] + Maximum relative change: 0.8000 ["ppkbot"] + + + OPTIMISATION ITERATION NO. : 2 + Model calls so far : 5 + Starting phi for this iteration: 1.68259E-03 + + Lambda = 2.5000 -----> + Phi = 3.27317E-04 ( 0.195 of starting phi) + + No more lambdas: phi is less than 0.3000 of starting phi + Lowest phi this iteration: 3.27317E-04 + + Current parameter values Previous parameter values + ppktop 3.88756 ppktop 3.49071 + ppkbot 0.774183 ppkbot 0.400000 + Maximum factor change: 1.935 ["ppkbot"] + Maximum relative change: 0.9355 ["ppkbot"] + + + OPTIMISATION ITERATION NO. : 3 + Model calls so far : 8 + Starting phi for this iteration: 3.27317E-04 + + Lambda = 1.2500 -----> + Phi = 2.67179E-05 ( 0.082 of starting phi) + + No more lambdas: phi is less than 0.3000 of starting phi + Lowest phi this iteration: 2.67179E-05 + + Current parameter values Previous parameter values + ppktop 4.69387 ppktop 3.88756 + ppkbot 0.842669 ppkbot 0.774183 + Maximum factor change: 1.207 ["ppktop"] + Maximum relative change: 0.2074 ["ppktop"] + + + OPTIMISATION ITERATION NO. : 4 + Model calls so far : 11 + Starting phi for this iteration: 2.67179E-05 + + Lambda = 0.62500 -----> + Phi = 1.00989E-05 ( 0.378 of starting phi) + + Lambda = 0.31250 -----> + Phi = 9.68202E-06 ( 0.362 of starting phi) + + Lambda = 0.15625 -----> + Phi = 9.51455E-06 ( 0.356 of starting phi) + + No more lambdas: relative phi reduction between lambdas less than 0.0300 + Lowest phi this iteration: 9.51455E-06 + + Current parameter values Previous parameter values + ppktop 4.81608 ppktop 4.69387 + ppkbot 0.785615 ppkbot 0.842669 + Maximum factor change: 1.073 ["ppkbot"] + Maximum relative change: 6.7706E-02 ["ppkbot"] + + + OPTIMISATION ITERATION NO. : 5 + Model calls so far : 16 + Starting phi for this iteration: 9.51455E-06 + + Lambda = 7.81250E-02 -----> + Phi = 9.42226E-06 ( 0.990 of starting phi) + + Lambda = 3.90625E-02 -----> + Phi = 9.42216E-06 ( 0.990 of starting phi) + + No more lambdas: relative phi reduction between lambdas less than 0.0300 + Lowest phi this iteration: 9.42216E-06 + Relative phi reduction between optimisation iterations less than 0.1000 + Switch to higher order derivatives calculation + + Current parameter values Previous parameter values + ppktop 4.82782 ppktop 4.81608 + ppkbot 0.788320 ppkbot 0.785615 + Maximum factor change: 1.003 ["ppkbot"] + Maximum relative change: 3.4430E-03 ["ppkbot"] + + + OPTIMISATION ITERATION NO. : 6 + Model calls so far : 20 + Starting phi for this iteration: 9.42216E-06 + + Lambda = 1.95312E-02 -----> + Phi = 9.42213E-06 ( 1.000 of starting phi) + + Lambda = 9.76562E-03 -----> + Phi = 9.42213E-06 ( 1.000 of starting phi) + + No more lambdas: relative phi reduction between lambdas less than 0.0300 + Lowest phi this iteration: 9.42213E-06 + + Current parameter values Previous parameter values + ppktop 4.82812 ppktop 4.82782 + ppkbot 0.788274 ppkbot 0.788320 + Maximum factor change: 1.000 ["ppktop"] + Maximum relative change: 6.1373E-05 ["ppktop"] + + + OPTIMISATION ITERATION NO. : 7 + Model calls so far : 26 + Starting phi for this iteration: 9.42213E-06 + + Lambda = 4.88281E-03 -----> + Phi = 9.42213E-06 ( 1.000 of starting phi) + + Lambda = 2.44141E-03 -----> + Phi = 9.42213E-06 ( 1.000 of starting phi) + + No more lambdas: relative phi reduction between lambdas less than 0.0300 + Lowest phi this iteration: 9.42213E-06 + + Current parameter values Previous parameter values + ppktop 4.82814 ppktop 4.82812 + ppkbot 0.788271 ppkbot 0.788274 + Maximum factor change: 1.000 ["ppktop"] + Maximum relative change: 3.8939E-06 ["ppktop"] + + Optimisation complete: relative parameter change less than 1.0000E-02 + over 3 successive iterations. + Total model calls: 32 + + The model has been run one final time using best parameters. + Thus all model input files contain best parameter values, and model + output files contain model results based on these parameters. + + + OPTIMISATION RESULTS + + + Parameters -----> + + Parameter Estimated 95% percent confidence limits + value lower limit upper limit + ppktop 4.82814 4.29574 5.36053 + ppkbot 0.788271 0.655277 0.921265 + + Note: confidence limits provide only an indication of parameter uncertainty. + They rely on a linearity assumption which may not extend as far in + parameter space as the confidence limits themselves - see PEST manual. + + See file pest.sen for parameter sensitivities. + + + Observations -----> + + Observation Measured Calculated Residual Weight Group + value value + ttopmid 0.138000 0.138151 -1.510200E-04 1.000 obsgroup + tbotmid 4.000000E-02 4.186851E-02 -1.868510E-03 1.000 obsgroup + tbotlef 2.200000E-02 2.028128E-02 1.718719E-03 1.000 obsgroup + tbotrig 2.200000E-02 2.028128E-02 1.718719E-03 1.000 obsgroup + + See file pest.res for more details of residuals in graph-ready format. + + See file pest.seo for composite observation sensitivities. + + + Objective function -----> + + Sum of squared weighted residuals (ie phi) = 9.4221E-06 + + + Correlation Coefficient -----> + + Correlation coefficient = 0.99959 + + + Analysis of residuals -----> + + All residuals:- + Number of residuals with non-zero weight = 4 + Mean value of non-zero weighted residuals = 3.5448E-04 + Maximum weighted residual [observation "tbotrig"] = 1.7187E-03 + Minimum weighted residual [observation "tbotmid"] = -1.8685E-03 + Standard variance of weighted residuals = 4.7111E-06 + Standard error of weighted residuals = 2.1705E-03 + + Note: the above variance was obtained by dividing the objective + function by the number of system degrees of freedom (ie. number of + observations with non-zero weight plus number of prior information + articles with non-zero weight minus the number of adjustable parameters.) + If the degrees of freedom is negative the divisor becomes + the number of observations with non-zero weight plus the number of + prior information items with non-zero weight. + + + K-L information statistics -----> + + + AIC = -45.83498 + AICC = indeterminate + BIC = -47.67609 + KIC = -42.97078 + + Parameter covariance matrix -----> + + ppktop ppkbot + ppktop 1.5308E-02 -5.4816E-04 + ppkbot -5.4816E-04 9.5526E-04 + + + Parameter correlation coefficient matrix -----> + + ppktop ppkbot + ppktop 1.000 -0.1433 + ppkbot -0.1433 1.000 + + + Normalized eigenvectors of parameter covariance matrix -----> + + Vector_1 Vector_2 + ppktop 3.8109E-02 -0.9993 + ppkbot 0.9993 3.8109E-02 + + + Eigenvalues -----> + + 9.3436E-04 1.5329E-02 diff --git a/modules/optimization/examples/bimaterial/main.i b/modules/optimization/examples/bimaterial/main.i index 9904b05d1fc9..a71e55e4716f 100644 --- a/modules/optimization/examples/bimaterial/main.i +++ b/modules/optimization/examples/bimaterial/main.i @@ -3,7 +3,7 @@ # PETSc-TAO optimisation is used to perform this inversion # # PETSc-TAO is set to uses "lmvm" and finite-difference approximations to the derivative of the objective function -# This means that for each PETSc-TAO iteration, 5 runs of model.i are performed: one for the current value of diffusivity_values, and four more for forward and backward finite differences around this point. (Hence, if there were n diffusivity values, PETSc-TAO would run model.i 1 + 2^n times per iteration). +# This means that for each PETSc-TAO iteration, 5 runs of model.i are performed: one for the current value of diffusivity_values, and four more for forward and backward finite differences around this point. (Hence, if there were n diffusivity values, PETSc-TAO would run model.i 1 + 2*n times per iteration). # # This file passes the diffusivity_values to model.i via the MultiAppReporterTransfer defined below # This file retrieves from model.i: the misfit of the temperature predicted by model.i and the experimental observations; the temperature at the observation points; the diffusivity_values (so they can be recorded in the output). This occurs through the MultiAppReporterTransfer defined below diff --git a/modules/optimization/examples/bimaterial/plot_results.py b/modules/optimization/examples/bimaterial/plot_results.py new file mode 100644 index 000000000000..1343c580360e --- /dev/null +++ b/modules/optimization/examples/bimaterial/plot_results.py @@ -0,0 +1,50 @@ +import os +import sys +import matplotlib.pyplot as plt + +plt.figure() + +model_calls = [] +phi = [] +with open("gold/pest.rec", "r") as f: + for line in f: + line = line.strip() + if line.startswith("Model calls so far"): + model_calls.append(int(line.split(":")[-1])) + if line.startswith("Starting phi for this iteration"): + phi.append(float(line.split(":")[-1])) +plt.semilogy(model_calls, phi, 'o-', label = "PEST") + +forwards = [0] +obj_fcn = [] +with open("gold/main_out_optInfo_0001_lmvm.csv", "r") as f: + data = f.readlines()[1:] +for line in data: + cnorm, current_iterate, function_value, gnorm, gradient_iterate, hessian_iterate, objective_iterate, xdiff = line.strip().split(",") + forwards.append(forwards[-1] + int(objective_iterate)) + obj_fcn.append(float(function_value)) +forwards.pop(0) +plt.semilogy(forwards, obj_fcn, 'o-', label = "TAO: LMVM") + +forwards = [0] +obj_fcn = [] +with open("gold/main_out_optInfo_0001_nm.csv", "r") as f: + data = f.readlines()[1:] +for line in data: + cnorm, current_iterate, function_value, gnorm, gradient_iterate, hessian_iterate, objective_iterate, xdiff = line.strip().split(",") + forwards.append(forwards[-1] + int(objective_iterate)) + obj_fcn.append(float(function_value)) +forwards.pop(0) +plt.semilogy(forwards, obj_fcn, 'o-', label = "TAO: NM") + + +plt.legend() +plt.grid() +plt.xlim([0, plt.xlim()[1]]) +plt.xlabel("Number of model solves") +plt.ylabel("Objective function") +plt.title("Comparison of inverse methods for the bimaterial example") +plt.show() +plt.savefig("comparison.png", bbox_inches = 'tight') + + diff --git a/modules/optimization/examples/bimaterial/tests b/modules/optimization/examples/bimaterial/tests index 0c6e03c21d9a..9cc053be8720 100644 --- a/modules/optimization/examples/bimaterial/tests +++ b/modules/optimization/examples/bimaterial/tests @@ -1,14 +1,14 @@ [Tests] -[nm] - type = CSVDiff - input = main.i - csvdiff = main_out_temperature_at_observation_points_0001.csv - heavy = true - cli_args = "Executioner/tao_solver=TAONM" - " Executioner/petsc_options_iname='-tao_gatol' " - " Executioner/petsc_options_value='0.0001'" - requirement = "XXXX shall be able to use Nelder Mead to invert for material properties" -[] + [nm] + type = CSVDiff + input = main.i + csvdiff = main_out_temperature_at_observation_points_0001.csv + heavy = true + cli_args = "Executioner/tao_solver=TAONM" + " Executioner/petsc_options_iname='-tao_gatol' " + " Executioner/petsc_options_value='0.0001'" + requirement = "XXXX shall be able to use Nelder Mead to invert for material properties" + [] [fd_grad] type = CSVDiff input = main.i