{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Run the fusion EasyVVUQ campaign\n", "\n", "Run an EasyVVUQ campaign to analyze the sensitivity of the temperature\n", "profile predicted by a simplified model of heat conduction in a\n", "tokamak plasma.\n", "\n", "This is done with SC." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2021-06-04T20:25:31.713466Z", "start_time": "2021-06-04T20:24:40.688493Z" }, "code_folding": [ 0 ] }, "outputs": [], "source": [ "# import packages that we will use\n", "\n", "import os\n", "import easyvvuq as uq\n", "import chaospy as cp\n", "import pickle\n", "import time\n", "import numpy as np\n", "import pandas as pd\n", "import matplotlib\n", "if not os.getenv(\"DISPLAY\"): matplotlib.use('Agg')\n", "import matplotlib.pylab as plt\n", "from IPython.display import display\n", "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2021-06-04T20:25:34.272217Z", "start_time": "2021-06-04T20:25:31.716919Z" }, "code_folding": [ 0, 2, 4 ] }, "outputs": [], "source": [ "# we need fipy -- install if not already available\n", "\n", "try:\n", " import fipy\n", "except ModuleNotFoundError:\n", " ! pip install future\n", " ! pip install fipy\n", " import fipy" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2021-06-04T20:25:34.279881Z", "start_time": "2021-06-04T20:25:34.275176Z" }, "code_folding": [ 0, 2 ] }, "outputs": [], "source": [ "# routine to write out (if needed) the fusion .template file\n", "\n", "def write_template(params):\n", " str = \"\"\n", " first = True\n", " for k in params.keys():\n", " if first:\n", " str += '{\"%s\": \"$%s\"' % (k,k) ; first = False\n", " else:\n", " str += ', \"%s\": \"$%s\"' % (k,k)\n", " str += '}'\n", " print(str, file=open('fusion.template','w'))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2021-06-04T20:25:34.304016Z", "start_time": "2021-06-04T20:25:34.282193Z" }, "code_folding": [ 0 ] }, "outputs": [], "source": [ "# define parameters of the fusion model\n", "def define_params():\n", " return {\n", " \"Qe_tot\": {\"type\": \"float\", \"min\": 1.0e6, \"max\": 50.0e6, \"default\": 2e6},\n", " \"H0\": {\"type\": \"float\", \"min\": 0.00, \"max\": 1.0, \"default\": 0},\n", " \"Hw\": {\"type\": \"float\", \"min\": 0.01, \"max\": 100.0, \"default\": 0.1},\n", " \"Te_bc\": {\"type\": \"float\", \"min\": 10.0, \"max\": 1000.0, \"default\": 100},\n", " \"chi\": {\"type\": \"float\", \"min\": 0.01, \"max\": 100.0, \"default\": 1},\n", " \"a0\": {\"type\": \"float\", \"min\": 0.2, \"max\": 10.0, \"default\": 1},\n", " \"R0\": {\"type\": \"float\", \"min\": 0.5, \"max\": 20.0, \"default\": 3},\n", " \"E0\": {\"type\": \"float\", \"min\": 1.0, \"max\": 10.0, \"default\": 1.5},\n", " \"b_pos\": {\"type\": \"float\", \"min\": 0.95, \"max\": 0.99, \"default\": 0.98},\n", " \"b_height\": {\"type\": \"float\", \"min\": 3e19, \"max\": 10e19, \"default\": 6e19},\n", " \"b_sol\": {\"type\": \"float\", \"min\": 2e18, \"max\": 3e19, \"default\": 2e19},\n", " \"b_width\": {\"type\": \"float\", \"min\": 0.005, \"max\": 0.025, \"default\": 0.01},\n", " \"b_slope\": {\"type\": \"float\", \"min\": 0.0, \"max\": 0.05, \"default\": 0.01},\n", " \"nr\": {\"type\": \"integer\", \"min\": 10, \"max\": 1000, \"default\": 100},\n", " \"dt\": {\"type\": \"float\", \"min\": 1e-3, \"max\": 1e3, \"default\": 100},\n", " \"out_file\": {\"type\": \"string\", \"default\": \"output.csv\"}\n", " }" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2021-06-04T20:25:34.320552Z", "start_time": "2021-06-04T20:25:34.305858Z" }, "code_folding": [ 0, 2, 17, 21, 28 ] }, "outputs": [], "source": [ "# define varying quantities\n", "def define_vary():\n", " vary_all = {\n", " \"Qe_tot\": cp.Uniform(1.8e6, 2.2e6),\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", " \"Te_bc\": cp.Uniform(80.0, 120.0),\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", " }\n", " vary_2 = {\n", " \"Qe_tot\": cp.Uniform(1.8e6, 2.2e6),\n", " \"Te_bc\": cp.Uniform(80.0, 120.0)\n", " }\n", " vary_5 = {\n", " \"Qe_tot\": cp.Uniform(1.8e6, 2.2e6),\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", " \"Te_bc\": cp.Uniform(80.0, 120.0)\n", " }\n", " vary_10 = {\n", " \"Qe_tot\": cp.Uniform(1.8e6, 2.2e6),\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", " \"Te_bc\": cp.Uniform(80.0, 120.0),\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", " }\n", " return vary_10" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2021-06-04T20:25:34.335304Z", "start_time": "2021-06-04T20:25:34.322416Z" }, "code_folding": [ 0 ] }, "outputs": [], "source": [ "# define a model to run the fusion code directly from python, expecting a dictionary and returning a dictionary\n", "def run_fusion_model(input):\n", " import json\n", " import fusion\n", " qois = [\"te\", \"ne\", \"rho\", \"rho_norm\"]\n", " del input['out_file']\n", " return {q: v for q,v in zip(qois, [t.tolist() for t in fusion.solve_Te(**input, plots=False, output=False)])}" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2021-06-04T20:25:34.351660Z", "start_time": "2021-06-04T20:25:34.337059Z" }, "code_folding": [ 0 ] }, "outputs": [], "source": [ "# routine to run a SC campaign\n", "\n", "def run_sc_case(sc_order=2, local=True, dask=True, batch_size=os.cpu_count(), use_files=True):\n", " \"\"\"\n", " Inputs:\n", " sc_order: order of the sc expansion\n", " local: if using Dask, whether to use the local option (True)\n", " dask: whether to use dask (True)\n", " batch_size: for the non Dask option, number of cases to run in parallel (16)\n", " Outputs:\n", " results_df: Pandas dataFrame containing inputs to and output from the model\n", " results: Results of the sc analysis\n", " times: Information about the elapsed time for the various phases of the calculation\n", " sc_order: sc_order \n", " count: number of sc samples\n", " \"\"\"\n", " \n", " if dask:\n", " if local:\n", " print('Running locally')\n", " import multiprocessing.popen_spawn_posix\n", " from dask.distributed import Client, LocalCluster\n", " cluster = LocalCluster(threads_per_worker=1)\n", " client = Client(cluster) # processes=True, threads_per_worker=1)\n", " else:\n", " print('Running using SLURM')\n", " from dask.distributed import Client\n", " from dask_jobqueue import SLURMCluster\n", " cluster = SLURMCluster(\n", " job_extra=['--qos=p.tok.openmp.2h', '--mail-type=end', '--mail-user=dpc@rzg.mpg.de', '-t 2:00:00'], \n", " queue='p.tok.openmp', \n", " cores=8, \n", " memory='8 GB',\n", " processes=8)\n", " cluster.scale(32)\n", " print(cluster)\n", " print(cluster.job_script())\n", " client = Client(cluster)\n", " print(client)\n", "\n", " else:\n", " import concurrent.futures\n", "# client = concurrent.futures.ProcessPoolExecutor(max_workers=batch_size)\n", " client = concurrent.futures.ThreadPoolExecutor(max_workers=batch_size)\n", "# client = None\n", " \n", " times = np.zeros(7)\n", "\n", " time_start = time.time()\n", " time_start_whole = time_start\n", " # Set up a fresh campaign called \"fusion_sc.\"\n", " my_campaign = uq.Campaign(name='fusion_sc.') \n", "\n", " # Define parameter space\n", " params = define_params()\n", "\n", " # Create an encoder and decoder for sc test app\n", " if use_files:\n", " encoder = uq.encoders.GenericEncoder(template_fname='fusion.template',\n", " delimiter='$',\n", " target_filename='fusion_in.json')\n", "\n", "\n", " decoder = uq.decoders.SimpleCSV(target_filename=\"output.csv\",\n", " output_columns=[\"te\", \"ne\", \"rho\", \"rho_norm\"])\n", "\n", " execute = uq.actions.ExecuteLocal('python3 %s/fusion_model.py fusion_in.json' % (os.getcwd()))\n", "\n", " actions = uq.actions.Actions(uq.actions.CreateRunDirectory('/tmp'), \n", " uq.actions.Encode(encoder), execute, uq.actions.Decode(decoder))\n", " else:\n", " actions = uq.actions.Actions(uq.actions.ExecutePython(run_fusion_model))\n", "\n", "\n", " # Add the app (automatically set as current app)\n", " my_campaign.add_app(name=\"fusion\", params=params, actions=actions)\n", "\n", " time_end = time.time()\n", " times[1] = time_end-time_start\n", " print('Time for phase 1 = %.3f' % (times[1]))\n", "\n", " time_start = time.time()\n", " # Associate a sampler with the campaign\n", " my_campaign.set_sampler(uq.sampling.SCSampler(vary=define_vary(), polynomial_order=sc_order))\n", " my_campaign.draw_samples()\n", " print('Number of samples = %s' % my_campaign.get_active_sampler().count)\n", "\n", " time_end = time.time()\n", " times[2] = time_end-time_start\n", " print('Time for phase 2 = %.3f' % (times[2]))\n", "\n", " time_start = time.time()\n", " # Perform the actions\n", " my_campaign.execute(pool=client).collate(progress_bar=True)\n", "\n", " if dask:\n", " client.close()\n", " client.shutdown()\n", "\n", " time_end = time.time()\n", " times[3] = time_end-time_start\n", " print('Time for phase 3 = %.3f' % (times[3]))\n", "\n", " time_start = time.time()\n", " # Collate the results\n", " results_df = my_campaign.get_collation_result()\n", "\n", " time_end = time.time()\n", " times[4] = time_end-time_start\n", " print('Time for phase 4 = %.3f' % (times[4]))\n", "\n", " time_start = time.time()\n", " # Post-processing analysis\n", " results = my_campaign.analyse(qoi_cols=[\"te\", \"ne\", \"rho\", \"rho_norm\"])\n", "\n", " time_end = time.time()\n", " times[5] = time_end-time_start\n", " print('Time for phase 5 = %.3f' % (times[5]))\n", "\n", " time_start = time.time()\n", " # Save the results\n", " pickle.dump(results, open('fusion_results.pickle','bw'))\n", " time_end = time.time()\n", " times[6] = time_end-time_start\n", " print('Time for phase 6 = %.3f' % (times[6]))\n", "\n", " times[0] = time_end - time_start_whole\n", "\n", " return results_df, results, times, sc_order, my_campaign.get_active_sampler().count" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2021-06-04T20:25:34.378499Z", "start_time": "2021-06-04T20:25:34.354265Z" }, "code_folding": [ 0, 2, 24, 46, 58, 72, 84 ] }, "outputs": [], "source": [ "# routines for plotting the results\n", "\n", "def plot_Te(results, title=None):\n", " # plot the calculated Te: mean, with std deviation, 1, 10, 90 and 99%\n", " plt.figure()\n", " rho = results.describe('rho', 'mean')\n", " plt.plot(rho, results.describe('te', 'mean'), 'b-', label='Mean')\n", " plt.plot(rho, results.describe('te', 'mean')-results.describe('te', 'std'), 'b--', label='+1 std deviation')\n", " plt.plot(rho, results.describe('te', 'mean')+results.describe('te', 'std'), 'b--')\n", " 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)\n", " try:\n", " plt.plot(rho, results.describe('te', '10%'), 'b:', label='10 and 90 percentiles')\n", " plt.plot(rho, results.describe('te', '90%'), 'b:')\n", " plt.fill_between(rho, results.describe('te', '10%'), results.describe('te', '90%'), color='b', alpha=0.1)\n", " plt.fill_between(rho, results.describe('te', '1%'), results.describe('te', '99%'), color='b', alpha=0.05)\n", " except:\n", " print('Problem with some of the percentiles')\n", " plt.legend(loc=0)\n", " plt.xlabel('rho [$m$]')\n", " plt.ylabel('Te [$eV$]')\n", " if not title is None: plt.title(title)\n", " plt.savefig('Te.png')\n", " plt.savefig('Te.pdf')\n", "\n", "def plot_ne(results, title=None):\n", " # plot the calculated ne: mean, with std deviation, 1, 10, 90 and 99%\n", " plt.figure()\n", " rho = results.describe('rho', 'mean')\n", " plt.plot(rho, results.describe('ne', 'mean'), 'b-', label='Mean')\n", " plt.plot(rho, results.describe('ne', 'mean')-results.describe('ne', 'std'), 'b--', label='+1 std deviation')\n", " plt.plot(rho, results.describe('ne', 'mean')+results.describe('ne', 'std'), 'b--')\n", " plt.fill_between(rho, results.describe('ne', 'mean')-results.describe('ne', 'std'), results.describe('ne', 'mean')+results.describe('ne', 'std'), color='b', alpha=0.2)\n", " try:\n", " plt.plot(rho, results.describe('ne', '10%'), 'b:', label='10 and 90 percentiles')\n", " plt.plot(rho, results.describe('ne', '90%'), 'b:')\n", " plt.fill_between(rho, results.describe('ne', '10%'), results.describe('ne', '90%'), color='b', alpha=0.1)\n", " plt.fill_between(rho, results.describe('ne', '1%'), results.describe('ne', '99%'), color='b', alpha=0.05)\n", " except:\n", " print('Problem with some of the percentiles')\n", " plt.legend(loc=0)\n", " plt.xlabel('rho [$m$]')\n", " plt.ylabel('ne [$m^{-3}$]')\n", " if not title is None: plt.title(title)\n", " plt.savefig('ne.png')\n", " plt.savefig('ne.pdf')\n", "\n", "def plot_sobols_first(results, title=None, field='te'):\n", " # plot the first Sobol results\n", " plt.figure()\n", " rho = results.describe('rho', 'mean')\n", " for k in results.sobols_first()[field].keys(): plt.plot(rho, results.sobols_first()[field][k], label=k)\n", " plt.legend(loc=0)\n", " plt.xlabel('rho [$m$]')\n", " plt.ylabel('sobols_first')\n", " if not title is None: plt.title(field + ': ' + title)\n", " plt.savefig('sobols_first_%s.png' % (field))\n", " plt.savefig('sobols_first_%s.pdf' % (field))\n", "\n", "def plot_sobols_second(results, title=None, field='te'):\n", " # plot the second Sobol results\n", " plt.figure()\n", " rho = results.describe('rho', 'mean')\n", " for k1 in results.sobols_second()[field].keys():\n", " for k2 in results.sobols_second()[field][k1].keys():\n", " plt.plot(rho, results.sobols_second()[field][k1][k2], label=k1+'/'+k2)\n", " plt.legend(loc=0, ncol=2)\n", " plt.xlabel('rho [$m$]')\n", " plt.ylabel('sobols_second')\n", " if not title is None: plt.title(field + ': ' + title)\n", " plt.savefig('sobols_second_%s.png' % (field))\n", " plt.savefig('sobols_second_%s.pdf' % (field))\n", "\n", "def plot_sobols_total(results, title=None, field='te'):\n", " # plot the total Sobol results\n", " plt.figure()\n", " rho = results.describe('rho', 'mean')\n", " for k in results.sobols_total()[field].keys(): plt.plot(rho, results.sobols_total()[field][k], label=k)\n", " plt.legend(loc=0)\n", " plt.xlabel('rho [$m$]')\n", " plt.ylabel('sobols_total')\n", " if not title is None: plt.title(field + ': ' + title)\n", " plt.savefig('sobols_total_%s.png' % (field))\n", " plt.savefig('sobols_total_%s.pdf' % (field))\n", "\n", "def plot_distribution(results, results_df, title=None):\n", " te_dist = results.raw_data['output_distributions']['te']\n", " rho_norm = results.describe('rho_norm', 'mean')\n", " for i in [np.maximum(0, int(i-1)) \n", " for i in np.linspace(0,1,5) * rho_norm.shape]:\n", " plt.figure()\n", " pdf_raw_samples = cp.GaussianKDE(results_df.te[i])\n", " pdf_kde_samples = cp.GaussianKDE(te_dist.samples[i])\n", " plt.hist(results_df.te[i], density=True, bins=50, label='histogram of raw samples', alpha=0.25)\n", " if hasattr(te_dist, 'samples'):\n", " plt.hist(te_dist.samples[i], density=True, bins=50, label='histogram of kde samples', alpha=0.25)\n", "\n", " plt.plot(np.linspace(pdf_raw_samples.lower, pdf_raw_samples.upper), pdf_raw_samples.pdf(np.linspace(pdf_raw_samples.lower, pdf_raw_samples.upper)), label='PDF (raw samples)')\n", " plt.plot(np.linspace(pdf_kde_samples.lower, pdf_kde_samples.upper), pdf_kde_samples.pdf(np.linspace(pdf_kde_samples.lower, pdf_kde_samples.upper)), label='PDF (kde samples)')\n", "\n", " plt.legend(loc=0)\n", " plt.xlabel('Te [$eV$]')\n", " if title is None:\n", " plt.title('Distributions for rho_norm = %0.4f' % (rho_norm[i]))\n", " else:\n", " plt.title('%s\\nDistributions for rho_norm = %0.4f' % (title, rho_norm[i]))\n", " plt.savefig('distribution_function_rho_norm=%0.4f.png' % (rho_norm[i]))\n", " plt.savefig('distribution_function_rho_norm=%0.4f.pdf' % (rho_norm[i]))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2021-07-16T11:18:30.251249Z", "start_time": "2021-06-04T20:25:34.380415Z" }, "code_folding": [ 0 ], "scrolled": true }, "outputs": [], "source": [ "# Calculate the stochastic collocation expansion for a range of orders\n", "\n", "if __name__ == '__main__':\n", " local = False # if True, use local cores; if False, use SLURM\n", " dask = False # if True, use DASK; if False, use a fall-back non-DASK option\n", "\n", " R = {}\n", " for sc_order in range(1, 4):\n", " R[sc_order] = {}\n", " (R[sc_order]['results_df'], \n", " R[sc_order]['results'], \n", " R[sc_order]['times'], \n", " R[sc_order]['order'], \n", " R[sc_order]['number_of_samples']) = run_sc_case(sc_order=sc_order, \n", " local=local, dask=dask, \n", " batch_size=7, use_files=False)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2021-07-16T11:20:48.322921Z", "start_time": "2021-07-16T11:18:30.253551Z" }, "code_folding": [] }, "outputs": [], "source": [ "# save the results\n", "\n", "if __name__ == '__main__':\n", "\n", " pickle.dump(R, open('collected_results.pickle','bw'))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2021-07-16T11:20:52.573368Z", "start_time": "2021-07-16T11:20:48.326442Z" }, "code_folding": [ 0 ], "scrolled": false }, "outputs": [], "source": [ "# produce a table of the time taken for various phases\n", "# the phases are:\n", "# 1: creation of campaign\n", "# 2: creation of samples\n", "# 3: running the cases\n", "# 4: calculation of statistics including Sobols\n", "# 5: returning of analysed results\n", "# 6: saving campaign and pickled results\n", "\n", "if __name__ == '__main__':\n", "\n", " Timings = pd.DataFrame(np.array([R[r]['times'] for r in list(R.keys())]), \n", " columns=['Total', 'Phase 1', 'Phase 2', 'Phase 3', 'Phase 4', 'Phase 5', 'Phase 6'], \n", " index=[R[r]['order'] for r in list(R.keys())])\n", " Timings.to_csv(open('Timings.csv', 'w'))\n", " display(Timings)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2021-07-16T11:21:05.954235Z", "start_time": "2021-07-16T11:20:52.578061Z" }, "code_folding": [ 0 ] }, "outputs": [], "source": [ "# plot the convergence of the mean and standard deviation to that of the highest order\n", "\n", "if __name__ == '__main__':\n", " last = -1\n", " O = [R[r]['order'] for r in list(R.keys())]\n", " if len(O[0:last]) > 0:\n", " plt.figure()\n", " plt.semilogy([o for o in O[0:last]],\n", " [np.sqrt(np.mean((R[o]['results'].describe('te', 'mean') - \n", " R[O[last]]['results'].describe('te', 'mean'))**2)) for o in O[0:last]],\n", " 'o-', label='mean')\n", " plt.semilogy([o for o in O[0:last]],\n", " [np.sqrt(np.mean((R[o]['results'].describe('te', 'std') - \n", " R[O[last]]['results'].describe('te', 'std'))**2)) for o in O[0:last]],\n", " 'o-', label='std')\n", " plt.xlabel('SC order')\n", " plt.ylabel('RMSerror compared to order=%s' % (O[last]))\n", " plt.legend(loc=0)\n", " plt.savefig('Convergence_mean_std.png')\n", " plt.savefig('Convergence_mean_std.pdf')" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2021-07-16T11:21:07.272635Z", "start_time": "2021-07-16T11:21:05.956222Z" }, "code_folding": [ 0 ] }, "outputs": [], "source": [ "# plot the convergence of the first sobol to that of the highest order\n", "\n", "if __name__ == '__main__':\n", " last = -1\n", " O = [R[r]['order'] for r in list(R.keys())]\n", " if len(O[0:last]) > 0:\n", " plt.figure()\n", " O = [R[r]['order'] for r in list(R.keys())]\n", " for v in list(R[O[last]]['results'].sobols_first('te').keys()):\n", " plt.semilogy([o for o in O[0:last]],\n", " [np.sqrt(np.mean((R[o]['results'].sobols_first('te')[v] - \n", " R[O[last]]['results'].sobols_first('te')[v])**2)) for o in O[0:last]],\n", " 'o-',\n", " label=v)\n", " plt.xlabel('SC order')\n", " plt.ylabel('RMSerror for 1st sobol compared to order=%s' % (O[last]))\n", " plt.legend(loc=0)\n", " plt.savefig('Convergence_sobol_first.png')\n", " plt.savefig('Convergence_sobol_first.pdf')" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2021-07-16T11:21:12.363744Z", "start_time": "2021-07-16T11:21:07.278208Z" }, "code_folding": [ 0 ], "scrolled": false }, "outputs": [], "source": [ "# plot a standard set of graphs for the highest order case\n", "\n", "if __name__ == '__main__':\n", " last = -1\n", " title = 'fusion test case with SC order = %i' % list(R.values())[last]['order']\n", " plot_Te(list(R.values())[last]['results'], title=title,)\n", " plot_ne(list(R.values())[last]['results'], title=title)\n", " plot_sobols_first(list(R.values())[last]['results'], title=title)\n", " try:\n", " plot_sobols_second(list(R.values())[last]['results'], title=title)\n", " except:\n", " print('Problem with sobols_second')\n", " try:\n", " plot_sobols_total(list(R.values())[last]['results'], title=title)\n", " except:\n", " print('Problem with sobols_total')\n", " try:\n", " plot_distribution(list(R.values())[last]['results'], list(R.values())[last]['results_df'], title=title)\n", " except:\n", " print('Problem with distribution')\n", " plot_sobols_first(list(R.values())[last]['results'], title=title, field='ne')\n", " try:\n", " plot_sobols_second(list(R.values())[last]['results'], title=title, field='ne')\n", " except:\n", " print('Problem with sobols_second')\n", " try:\n", " plot_sobols_total(list(R.values())[last]['results'], title=title, field='ne')\n", " except:\n", " print('Problem with sobols_total')" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2021-07-16T11:23:21.604268Z", "start_time": "2021-07-16T11:21:12.365988Z" }, "code_folding": [ 0 ] }, "outputs": [], "source": [ "# prepare the test data\n", "\n", "if __name__ == '__main__':\n", "\n", " test_campaign = uq.Campaign(name='fusion_pce.') \n", "\n", " # Add the app (automatically set as current app)\n", " test_campaign.add_app(name=\"fusion\", params=define_params(), \n", " actions=uq.actions.Actions(uq.actions.ExecutePython(run_fusion_model)))\n", "\n", " # Associate a sampler with the campaign\n", " test_campaign.set_sampler(uq.sampling.quasirandom.LHCSampler(vary=define_vary(), count=100))\n", "\n", " # Perform the actions\n", " test_campaign.execute(nsamples=1000).collate(progress_bar=True)\n", "\n", " # Collate the results\n", " test_df = test_campaign.get_collation_result()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2021-07-18T12:46:20.388796Z", "start_time": "2021-07-16T11:23:21.606413Z" }, "code_folding": [] }, "outputs": [], "source": [ "# calculate the SC surrogates\n", "if __name__ == '__main__':\n", "\n", " test_points = test_df[test_campaign.get_active_sampler().vary.get_keys()]\n", " test_results = test_df['te'].values\n", " test_predictions = {}\n", " for i in list(R.keys()):\n", " test_predictions[i] = np.squeeze(np.array(R[i]['results'].surrogate()(test_points)['te']))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2021-07-18T12:46:21.334748Z", "start_time": "2021-07-18T12:46:20.391462Z" }, "code_folding": [ 0 ] }, "outputs": [], "source": [ "# plot the performance of the SC surrogates\n", "if __name__ == '__main__':\n", "\n", " for i in list(R.keys()):\n", " plt.semilogy(R[i]['results'].describe('rho', 'mean'), \n", " np.sqrt(((test_predictions[i] - test_results)**2).mean(axis=0)) / test_results.mean(axis=0), \n", " label='SC order %s with %s samples' % (R[i]['order'], R[i]['number_of_samples']))\n", " plt.xlabel('rho [m]') ; plt.ylabel('fractional RMS for predicted Te') ; plt.legend(loc=0)\n", " plt.title('Performance of the SC surrogate')\n", " plt.savefig('SC_surrogate.png')\n", " plt.savefig('SC_surrogate.pdf')" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2021-07-18T12:46:21.873529Z", "start_time": "2021-07-18T12:46:21.336591Z" }, "code_folding": [] }, "outputs": [], "source": [ "# plot the convergence of the surrogate based on 1000 random points\n", "if __name__ == '__main__':\n", " _o = []\n", " _RMS = []\n", " for r in R.values():\n", " _RMS.append((np.sqrt((((test_predictions[r['order']] - test_results) / test_results)**2).mean())))\n", " _o.append(r['order'])\n", "\n", " plt.figure()\n", " plt.semilogy(_o, _RMS, 'o-')\n", " plt.xlabel('SC order')\n", " plt.ylabel('fractional RMS error for the SC surrogate')\n", " plt.legend(loc=0)\n", " plt.savefig('Convergence_SC_surrogate.png')\n", " plt.savefig('Convergence_SC_surrogate.pdf')" ] } ], "metadata": { "jupytext": { "cell_metadata_filter": "-all", "executable": " /usr/bin/env python", "main_language": "python", "notebook_metadata_filter": "-all" }, "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.7" }, "latex_envs": { "LaTeX_envs_menu_present": true, "autoclose": false, "autocomplete": true, "bibliofile": "biblio.bib", "cite_by": "apalike", "current_citInitial": 1, "eqLabelWithNumbers": true, "eqNumInitial": 1, "hotkeys": { "equation": "Ctrl-E", "itemize": "Ctrl-I" }, "labels_anchors": false, "latex_user_defs": false, "report_style_numbering": false, "user_envs_cfg": false } }, "nbformat": 4, "nbformat_minor": 4 }