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.
#!pip install EasyVVUQ
#!pip install future
#!pip install fipy
# 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
# 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"}
}
# 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)
# 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)
"""
' other possible quantities to vary\n "H0": cp.Uniform(0.0, 0.2),\n "Hw": cp.Uniform(0.1, 0.5),\n "chi": cp.Uniform(0.8, 1.2), \n\n "a0": cp.Uniform(0.9, 1.1), \n "R0": cp.Uniform(2.7, 3.3), \n "E0": cp.Uniform(1.4, 1.6), \n "b_pos": cp.Uniform(0.95, 0.99), \n "b_height": cp.Uniform(5e19, 7e19), \n "b_sol": cp.Uniform(1e19, 3e19), \n "b_width": cp.Uniform(0.015, 0.025), \n "b_slope": cp.Uniform(0.005, 0.020)\n'
# Set up a fresh campaign called "fusion_pce."
my_campaign = uq.Campaign(name='fusion_pce.')
# Add the app (automatically set as current app)
my_campaign.add_app(name="fusion", params=params, actions=actions)
# Associate a sampler with the campaign
my_campaign.set_sampler(uq.sampling.PCESampler(vary=vary, polynomial_order=3))
my_campaign.execute().collate(progress_bar=True)
100%|██████████| 16/16 [00:01<00:00, 10.25it/s]
results_df = my_campaign.get_collation_result()
results = my_campaign.analyse(qoi_cols=["te", "ne", "rho", "rho_norm"])
# Get Descriptive Statistics
rho = results.describe('rho', 'mean')
rho_norm = results.describe('rho_norm', 'mean')
# 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);
# 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);
# 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');
# 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);
# 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');