.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "advanced_examples/6061T6_anisotropic_calibration/plot_6061T6_h_laplace_study.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_advanced_examples_6061T6_anisotropic_calibration_plot_6061T6_h_laplace_study.py: 6061T6 aluminum calibration uncertainty quantification ------------------------------------------------------ In this example, we will use MatCal's :class:`~matcal.core.parameter_studies.LaplaceStudy` to estimate the parameter uncertainty for the calibration. This uncertainty quantification study differs from the one in :ref:`304L stainless steel viscoplastic calibration uncertainty quantification` because many of the parameters are directly correlated to other parameters in the model. Specifically, the temperature and anisotropy parameters, are multipliers of the yield stress and Voce hardening parameters. As a result, we will assume all parameter uncertainty can be attributed to the yield stress and Voce hardening parameters alone. This will significantly reduce the cost of the finite difference calculations needed for the :class:`~matcal.core.parameter_studies.LaplaceStudy` and ensure robustness of the method. We want the uncertainty in these three parameters to account for all uncertainty in the room temperature experiments, so we include these models in the uncertainty study. .. warning:: The LaplaceStudy is still in development and may not accurately attribute uncertainty to to the parameters. Always verify results before use. .. note:: Useful Documentation links: #. :ref:`Uniaxial Tension Models` #. :ref:`Top Hat Shear Model` #. :class:`~matcal.core.parameter_studies.LaplaceStudy` To begin, we import the tools we need for this study. .. GENERATED FROM PYTHON SOURCE LINES 35-48 .. code-block:: Python from matcal import * from site_matcal.sandia.computing_platforms import is_sandia_cluster, get_sandia_computing_platform from site_matcal.sandia.tests.utilities import MATCAL_WCID import numpy as np import matplotlib.pyplot as plt plt.rc('text', usetex=True) plt.rc('font', family='serif') plt.rc('font', size=12) figsize = (4,3) .. GENERATED FROM PYTHON SOURCE LINES 49-55 Next, we import the data and remove any uncalibrated data from the :class:`~matcal.core.data.DataCollection` objects. We do this in place of weighting, because zeros in the residuals can cause scaling and conditioning issues in the linear algebra required for the study. .. GENERATED FROM PYTHON SOURCE LINES 55-59 .. code-block:: Python tension_data_collection = BatchDataImporter("aluminum_6061_data/" "uniaxial_tension/processed_data/" "cleaned_[CANM]*.csv",).batch .. GENERATED FROM PYTHON SOURCE LINES 60-63 For the room temperature tension data, we remove data in the elastic region and in regions of unloading to match what was included in the objective for the calibration. .. GENERATED FROM PYTHON SOURCE LINES 63-79 .. code-block:: Python down_selected_tension_data = DataCollection("down selected data") for state in tension_data_collection.keys(): for index, data in enumerate(tension_data_collection[state]): stresses = data["engineering_stress"] strains = data["engineering_strain"] peak_index = np.argmax(stresses) peak_strain = strains[peak_index] peak_stress = stresses[peak_index] data_to_keep = (((strains>peak_strain) & (stresses > 0.89*peak_stress)) | (strains>0.005) & (strains < peak_strain)) down_selected_tension_data.add(data[data_to_keep]) down_selected_tension_data = scale_data_collection(down_selected_tension_data, "engineering_stress", 1000) down_selected_tension_data.remove_field("time") .. GENERATED FROM PYTHON SOURCE LINES 80-85 With the down-selected tension data created, we create the :class:`~matcal.sierra.models.RoundUniaxialTensionModel` as we did in :ref:`6061T6 aluminum calibration with anisotropic yield`, and add the :class:`~matcal.core.data.DataCollection` that we created as the model boundary condition data. .. GENERATED FROM PYTHON SOURCE LINES 85-119 .. code-block:: Python material_filename = "hill_plasticity_temperature_dependent.inc" material_model = "hill_plasticity" material_name = "ductile_failure_6061T6" sierra_material = Material(material_name, material_filename, material_model) gauge_radius = 0.125 element_size = gauge_radius/8 geo_params = {"extensometer_length": 1.0, "gauge_length": 1.25, "gauge_radius": gauge_radius, "grip_radius": 0.25, "total_length": 4, "fillet_radius": 0.188, "taper": 0.0015, "necking_region":0.375, "element_size": element_size, "mesh_method":3, "grip_contact_length":1} tension_model = RoundUniaxialTensionModel(sierra_material, **geo_params) tension_model.set_name("tension_model") tension_model.add_boundary_condition_data(down_selected_tension_data) tension_model.set_allowable_load_drop_factor(0.70) tension_model.set_boundary_condition_scale_factor(1.5) if is_sandia_cluster(): tension_model.run_in_queue(MATCAL_WCID, 1) tension_model.continue_when_simulation_fails() platform = get_sandia_computing_platform() num_cores = platform.get_processors_per_node() else: num_cores = 8 tension_model.set_number_of_cores(num_cores) .. GENERATED FROM PYTHON SOURCE LINES 120-122 Similarly, we import the top hat data and down select the data of interest for the residuals. .. GENERATED FROM PYTHON SOURCE LINES 122-137 .. code-block:: Python top_hat_data_collection = BatchDataImporter("aluminum_6061_data/" "top_hat_shear/processed_data/cleaned_*.csv").batch for state, state_data_list in top_hat_data_collection.items(): for index, data in enumerate(state_data_list): max_load_arg = np.argmax(data["load"]) # This slicing procedure removes the data after peak load # and before displacements of 0.005". data = data[data["time"] < data["time"][max_load_arg]] data = data[data["displacement"] > 0.005] # This one removes the data after a displacement of 0.02" # and reassigns the modified data to the # DataCollection top_hat_data_collection[state][index] = data[data["displacement"] < 0.02] top_hat_data_collection.remove_field("time") .. GENERATED FROM PYTHON SOURCE LINES 138-141 With the data prepared, we can build the model as we did in the previous example :ref:`6061T6 aluminum calibration with anisotropic yield`. .. GENERATED FROM PYTHON SOURCE LINES 141-159 .. code-block:: Python top_hat_geo_params = {"total_height":1.25, "base_height":0.75, "trapezoid_angle": 10.0, "top_width": 0.417*2, "base_width": 1.625, "base_bottom_height": (0.75-0.425), "thickness":0.375, "external_radius": 0.05, "internal_radius": 0.05, "hole_height": 0.3, "lower_radius_center_width":0.390*2, "localization_region_scale":0.0, "element_size":0.005, "numsplits":1} top_hat_model = TopHatShearModel(sierra_material, **top_hat_geo_params) top_hat_model.set_name('top_hat_shear') .. GENERATED FROM PYTHON SOURCE LINES 160-162 Next, we set its allowable load drop factor and provide boundary condition data. .. GENERATED FROM PYTHON SOURCE LINES 162-165 .. code-block:: Python top_hat_model.set_allowable_load_drop_factor(0.05) top_hat_model.add_boundary_condition_data(top_hat_data_collection) .. GENERATED FROM PYTHON SOURCE LINES 166-168 Lastly, we setup the platform information for running the model. .. GENERATED FROM PYTHON SOURCE LINES 168-173 .. code-block:: Python top_hat_model.set_number_of_cores(num_cores*2) if is_sandia_cluster(): top_hat_model.run_in_queue(MATCAL_WCID, 1) top_hat_model.continue_when_simulation_fails() .. GENERATED FROM PYTHON SOURCE LINES 174-180 We now create the objectives for the calibration. Both models are compared to the data using a :class:`~matcal.core.objective.CurveBasedInterpolatedObjective`. The tension specimen is calibrated to the engineering stress/strain data and the top hat specimen is calibrated to the load-displacement data. .. GENERATED FROM PYTHON SOURCE LINES 180-185 .. code-block:: Python tension_objective = CurveBasedInterpolatedObjective("engineering_strain", "engineering_stress") tension_objective.set_name("engineering_stress_strain_obj") top_hat_objective = CurveBasedInterpolatedObjective("displacement", "load") top_hat_objective.set_name("load_displacement_obj") .. GENERATED FROM PYTHON SOURCE LINES 186-191 We now create our parameters for the study. The study parameters are the ``yield_stress``, ``hardening`` and ``b`` parameters from :ref:`6061T6 aluminum calibration with anisotropic yield` with their current value set to their calibration values. .. GENERATED FROM PYTHON SOURCE LINES 191-198 .. code-block:: Python RT_calibrated_params = matcal_load("anisotropy_parameters.serialized") yield_stress = Parameter("yield_stress", 15, 50, RT_calibrated_params.pop("yield_stress")) hardening = Parameter("hardening", 0, 60, RT_calibrated_params.pop("hardening")) b = Parameter("b", 10, 40, RT_calibrated_params.pop("b")) .. GENERATED FROM PYTHON SOURCE LINES 199-201 To simplify setting up the laplace study, we put all the parameters in a :class:`~matcal.core.parameters.ParameterCollection`. .. GENERATED FROM PYTHON SOURCE LINES 201-204 .. code-block:: Python pc = ParameterCollection("uncertain_params", yield_stress, hardening, b) .. GENERATED FROM PYTHON SOURCE LINES 205-207 We also need the anisotropy so we store those parameters with the current value equal to the calibrated parameter values from the calibration step. .. GENERATED FROM PYTHON SOURCE LINES 207-218 .. code-block:: Python R22 = Parameter("R22", 0.8, 1.15, RT_calibrated_params["R22"]) R33 = Parameter("R33", 0.8, 1.15, RT_calibrated_params["R33"]) R12 = Parameter("R12", 0.8, 1.15, RT_calibrated_params["R12"]) R23 = Parameter("R23", 0.8, 1.15, RT_calibrated_params["R23"]) R31 = Parameter("R31", 0.8, 1.15, RT_calibrated_params["R31"]) .. GENERATED FROM PYTHON SOURCE LINES 219-224 The anisotropy parameters and temperature dependence parameters from :ref:`6061T6 aluminum temperature dependent calibration` will be added as model constants because they are being treated as deterministic and are still required for the models. They are added for the two models for this study. .. GENERATED FROM PYTHON SOURCE LINES 224-230 .. code-block:: Python high_temp_calibrated_params = matcal_load("temperature_dependent_parameters.serialized") tension_model.add_constants(**high_temp_calibrated_params, **RT_calibrated_params) top_hat_model.add_constants(**high_temp_calibrated_params, **RT_calibrated_params) .. GENERATED FROM PYTHON SOURCE LINES 231-233 Now we can create our laplace study and add our two evaluation sets. .. GENERATED FROM PYTHON SOURCE LINES 233-240 .. code-block:: Python laplace_study = LaplaceStudy(pc) laplace_study.set_parameter_center(**pc.get_current_value_dict()) laplace_study.set_working_directory("laplace_study", remove_existing=True) laplace_study.set_core_limit(250) laplace_study.add_evaluation_set(tension_model, tension_objective, down_selected_tension_data) laplace_study.add_evaluation_set(top_hat_model, top_hat_objective, top_hat_data_collection) .. GENERATED FROM PYTHON SOURCE LINES 241-250 Laplace study specific options include :meth:`~matcal.core.parameter_studies.LaplaceStudy.set_step_size` to set the finite difference step size and :meth:`~matcal.core.parameter_studies.LaplaceStudy.set_noise_estimate` for setting the estimated amount of noise in the data. We set the finite difference step size to one order of magnitude less than the default. Results are likely sensitive to this value for practical problems, and re-running the study with different values may be required. .. GENERATED FROM PYTHON SOURCE LINES 250-251 .. code-block:: Python laplace_study.set_step_size(1e-4) .. GENERATED FROM PYTHON SOURCE LINES 252-262 For this study type, providing an inaccurate noise estimate can result in unreasonable solutions. .. warning:: Appropriately handling the noise estimate is an active area of research. If attempted, some iteration may be required to find an valid estimate for noise. This can be done by running the study once to evaluate the model response and then re-running the study as a restart after changing the noise estimate or by calling :meth:`~matcal.core.parameter_studies.LaplaceStudy.update_laplace_estimate`. .. GENERATED FROM PYTHON SOURCE LINES 262-265 .. code-block:: Python laplace_study.set_noise_estimate(1e-2) results = laplace_study.launch() .. GENERATED FROM PYTHON SOURCE LINES 266-276 After the study completes, there are two results of concern: #. The estimated parameter covariance - calculated directly from the residual magnitude and the gradients of the residuals w.r.t. the parameters. #. The fitted parameter covariance - an optimized covariance that ensures the the covariance of the parameters is representative of the uncertainty due to model form error. This corrects the estimated parameter covariance using the objective described in :ref:`Laplace Approximation: Error Attributed to Model Error` We print both of these values below and save the results to be used in the next step of this example. .. GENERATED FROM PYTHON SOURCE LINES 276-280 .. code-block:: Python print("Initial covariance estimate:\n", results.estimated_parameter_covariance) print("Calibrated covariance estimate:\n", results.fitted_parameter_covariance) matcal_save("laplace_study_results.joblib", results) .. rst-class:: sphx-glr-script-out .. code-block:: none Initial covariance estimate: [[ 113.66447863 -20.90050619 -422.18930642] [ -20.90050619 5.62379061 70.41401601] [-422.18930642 70.41401601 1607.23724872]] Calibrated covariance estimate: [[ 125.77931866 -22.62114298 -485.48255954] [ -22.62114298 5.84238658 77.55667787] [-485.48255954 77.55667787 1915.39817462]] .. GENERATED FROM PYTHON SOURCE LINES 281-286 As noted above, the results can be sensitive to the estimated noise. To illustrate this point, we re-run the study results processing with updated noise estimates and print the results. Before updating the results, we save the previous results as copy of themselves because the update just updates the values on the results object. .. GENERATED FROM PYTHON SOURCE LINES 286-300 .. code-block:: Python import copy results = copy.deepcopy(results) results_high_noise = laplace_study.update_laplace_estimate(1e-1) results_high_noise = copy.deepcopy(results_high_noise) print("Initial covariance estimate noise set to 1e-2:\n", results.estimated_parameter_covariance) print("Calibrated covariance estimate noise set to 1e-2:\n", results.fitted_parameter_covariance) print("Initial covariance estimate noise set to 1e-1:\n", results_high_noise.estimated_parameter_covariance) print("Calibrated covariance estimate noise set to 1e-1:\n", results_high_noise.fitted_parameter_covariance) .. rst-class:: sphx-glr-script-out .. code-block:: none Initial covariance estimate noise set to 1e-2: [[ 113.66447863 -20.90050619 -422.18930642] [ -20.90050619 5.62379061 70.41401601] [-422.18930642 70.41401601 1607.23724872]] Calibrated covariance estimate noise set to 1e-2: [[ 125.77931866 -22.62114298 -485.48255954] [ -22.62114298 5.84238658 77.55667787] [-485.48255954 77.55667787 1915.39817462]] Initial covariance estimate noise set to 1e-1: [[ 113.66447862 -20.90050619 -422.18930638] [ -20.90050619 5.62379061 70.414016 ] [-422.18930638 70.414016 1607.23724857]] Calibrated covariance estimate noise set to 1e-1: [[ 116.61950646 -21.1831117 -430.32258055] [ -21.1831117 5.65767236 72.91811568] [-430.32258055 72.91811568 1626.45894062]] .. GENERATED FROM PYTHON SOURCE LINES 301-308 Note the difference in the result. This highlights the sensitivity of the method to the noise estimate. Some iteration may be required to obtain a useful result. Next, we sample the multivariate normal provided by the study covariance and calibrated values as the mean and visualize the results using seaborn's KDE pair plot .. GENERATED FROM PYTHON SOURCE LINES 308-321 .. code-block:: Python num_samples=5000 uncertain_param_sets = sample_multivariate_normal(num_samples, results.mean.to_list(), results.fitted_parameter_covariance, 12345, pc.get_item_names()) import seaborn as sns import pandas as pd sns.pairplot(data=pd.DataFrame(uncertain_param_sets), kind="kde" ) plt.show() # From this plot, we can see the uncertainty is considerably overestimated # and could result in unphysical values of the parameters. This method is still # work in progress for models with significant model form error. .. image-sg:: /advanced_examples/6061T6_anisotropic_calibration/images/sphx_glr_plot_6061T6_h_laplace_study_001.png :alt: plot 6061T6 h laplace study :srcset: /advanced_examples/6061T6_anisotropic_calibration/images/sphx_glr_plot_6061T6_h_laplace_study_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none /gpfs/knkarls/projects/matcal-stable/external_matcal/matcal/core/parameter_studies.py:919: RuntimeWarning: covariance is not symmetric positive-semidefinite. samples = rng.multivariate_normal(mean, sigma, nsamples).T .. rst-class:: sphx-glr-timing **Total running time of the script:** (25 minutes 24.911 seconds) .. _sphx_glr_download_advanced_examples_6061T6_anisotropic_calibration_plot_6061T6_h_laplace_study.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_6061T6_h_laplace_study.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_6061T6_h_laplace_study.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_6061T6_h_laplace_study.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_