#!/usr/bin/env python # coding: utf-8 # # EasyVVUQ fusion tutorial # # Run an EasyVVUQ campaign to analyze the sensitivity of the temperature # profile predicted by a simplified model of heat conduction in a # tokamak plasma. # # This is done with PCE. # In[1]: #!pip install EasyVVUQ #!pip install future #!pip install fipy # In[2]: # import packages that we will use import os import easyvvuq as uq import chaospy as cp import pickle import time import numpy as np import matplotlib.pylab as plt # In[3]: # Define parameter space params = { "Qe_tot": {"type": "float", "min": 1.0e6, "max": 50.0e6, "default": 2e6}, "H0": {"type": "float", "min": 0.00, "max": 1.0, "default": 0}, "Hw": {"type": "float", "min": 0.01, "max": 100.0, "default": 0.1}, "Te_bc": {"type": "float", "min": 10.0, "max": 1000.0, "default": 100}, "chi": {"type": "float", "min": 0.01, "max": 100.0, "default": 1}, "a0": {"type": "float", "min": 0.2, "max": 10.0, "default": 1}, "R0": {"type": "float", "min": 0.5, "max": 20.0, "default": 3}, "E0": {"type": "float", "min": 1.0, "max": 10.0, "default": 1.5}, "b_pos": {"type": "float", "min": 0.95, "max": 0.99, "default": 0.98}, "b_height": {"type": "float", "min": 3e19, "max": 10e19, "default": 6e19}, "b_sol": {"type": "float", "min": 2e18, "max": 3e19, "default": 2e19}, "b_width": {"type": "float", "min": 0.005, "max": 0.025, "default": 0.01}, "b_slope": {"type": "float", "min": 0.0, "max": 0.05, "default": 0.01}, "nr": {"type": "integer", "min": 10, "max": 1000, "default": 100}, "dt": {"type": "float", "min": 1e-3, "max": 1e3, "default": 100}, "out_file": {"type": "string", "default": "output.csv"} } # In[4]: # Create an encoder, decoder and collater for PCE test app encoder = uq.encoders.GenericEncoder(template_fname='fusion.template', delimiter='$', target_filename='fusion_in.json') decoder = uq.decoders.SimpleCSV(target_filename="output.csv", output_columns=["te", "ne", "rho", "rho_norm"]) execute = uq.actions.ExecuteLocal('python3 %s/fusion_model.py fusion_in.json' % (os.getcwd())) actions = uq.actions.Actions(uq.actions.CreateRunDirectory('/tmp'), uq.actions.Encode(encoder), execute, uq.actions.Decode(decoder)) #collater = uq.collate.AggregateSamples(average=False) # In[5]: # Create the sampler (here simplified to two uncertain quantities) vary = { "Qe_tot": cp.Uniform(1.8e6, 2.2e6), "Te_bc": cp.Uniform(80.0, 120.0) } """ other possible quantities to vary "H0": cp.Uniform(0.0, 0.2), "Hw": cp.Uniform(0.1, 0.5), "chi": cp.Uniform(0.8, 1.2), "a0": cp.Uniform(0.9, 1.1), "R0": cp.Uniform(2.7, 3.3), "E0": cp.Uniform(1.4, 1.6), "b_pos": cp.Uniform(0.95, 0.99), "b_height": cp.Uniform(5e19, 7e19), "b_sol": cp.Uniform(1e19, 3e19), "b_width": cp.Uniform(0.015, 0.025), "b_slope": cp.Uniform(0.005, 0.020) """ # In[6]: # Set up a fresh campaign called "fusion_pce." my_campaign = uq.Campaign(name='fusion_pce.') # In[7]: # Add the app (automatically set as current app) my_campaign.add_app(name="fusion", params=params, actions=actions) # In[8]: # Associate a sampler with the campaign my_campaign.set_sampler(uq.sampling.PCESampler(vary=vary, polynomial_order=3)) # In[9]: my_campaign.execute().collate(progress_bar=True) # In[10]: results_df = my_campaign.get_collation_result() # In[11]: results = my_campaign.analyse(qoi_cols=["te", "ne", "rho", "rho_norm"]) # In[12]: # Get Descriptive Statistics rho = results.describe('rho', 'mean') rho_norm = results.describe('rho_norm', 'mean') # In[13]: # plot the calculated Te: mean, with std deviation, 10 and 90% and range plt.figure() plt.plot(rho, results.describe('te', 'mean'), 'b-', label='Mean') plt.plot(rho, results.describe('te', 'mean')-results.describe('te', 'std'), 'b--', label='+1 std deviation') plt.plot(rho, results.describe('te', 'mean')+results.describe('te', 'std'), 'b--') plt.fill_between(rho, results.describe('te', 'mean')-results.describe('te', 'std'), results.describe('te', 'mean')+results.describe('te', 'std'), color='b', alpha=0.2) plt.plot(rho, results.describe('te', '10%'), 'b:', label='10 and 90 percentiles') plt.plot(rho, results.describe('te', '90%'), 'b:') plt.fill_between(rho, results.describe('te', '10%'), results.describe('te', '90%'), color='b', alpha=0.1) plt.fill_between(rho, results.describe('te', 'min'), results.describe('te', 'max'), color='b', alpha=0.05) plt.legend(loc=0) plt.xlabel('rho [m]') plt.ylabel('Te [eV]') plt.title(my_campaign.campaign_dir); # In[14]: # plot the first Sobol results plt.figure() for k in results.sobols_first()['te'].keys(): plt.plot(rho, results.sobols_first()['te'][k], label=k) plt.legend(loc=0) plt.xlabel('rho [m]') plt.ylabel('sobols_first') plt.title(my_campaign.campaign_dir); # In[15]: # plot the second Sobol results plt.figure() for k1 in results.sobols_second()['te'].keys(): for k2 in results.sobols_second()['te'][k1].keys(): plt.plot(rho, results.sobols_second()['te'][k1][k2], label=k1+'/'+k2) plt.legend(loc=0) plt.xlabel('rho [m]') plt.ylabel('sobols_second') plt.title(my_campaign.campaign_dir+'\n'); # In[16]: # plot the total Sobol results plt.figure() for k in results.sobols_total()['te'].keys(): plt.plot(rho, results.sobols_total()['te'][k], label=k) plt.legend(loc=0) plt.xlabel('rho [m]') plt.ylabel('sobols_total') plt.title(my_campaign.campaign_dir); # In[17]: # plot the distributions plt.figure(figsize=(12,6)) for i, D in enumerate(results.raw_data['output_distributions']['te'].samples): pdf_kde_samples = cp.GaussianKDE(D) _Te = np.linspace(pdf_kde_samples.lower, pdf_kde_samples.upper[0], 101) plt.loglog(_Te, pdf_kde_samples.pdf(_Te), 'b-', alpha=0.25) plt.loglog(results.describe('te', 'mean')[i], pdf_kde_samples.pdf(results.describe('te', 'mean')[i]), 'bo') plt.loglog(results.describe('te', 'mean')[i]-results.describe('te', 'std')[i], pdf_kde_samples.pdf(results.describe('te', 'mean')[i]-results.describe('te', 'std')[i]), 'b*') plt.loglog(results.describe('te', 'mean')[i]+results.describe('te', 'std')[i], pdf_kde_samples.pdf(results.describe('te', 'mean')[i]+results.describe('te', 'std')[i]), 'b*') plt.loglog(results.describe('te', '10%')[i], pdf_kde_samples.pdf(results.describe('te', '10%')[i]), 'b+') plt.loglog(results.describe('te', '90%')[i], pdf_kde_samples.pdf(results.describe('te', '90%')[i]), 'b+') plt.loglog(results.describe('te', '1%')[i], pdf_kde_samples.pdf(results.describe('te', '1%')[i]), 'bs') plt.loglog(results.describe('te', '99%')[i], pdf_kde_samples.pdf(results.describe('te', '99%')[i]), 'bs') plt.xlabel('Te') plt.ylabel('distribution function') plt.title(my_campaign.campaign_dir) plt.ylim(1e-4,1e-1) plt.savefig('distribution_functions.png');