diff --git a/CHANGELOG.md b/CHANGELOG.md index b3a3a23a06..9bcc47938b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -44,6 +44,7 @@ ## Breaking changes +- Moved `results` to separate repositories ([#761](https://github.com/pybamm-team/PyBaMM/pull/761)) - The parameters "Bruggeman coefficient" must now be specified separately as "Bruggeman coefficient (electrolyte)" and "Bruggeman coefficient (electrode)" - The current classes (`GetConstantCurrent`, `GetUserCurrent` and `GetUserData`) have now been removed. Please refer to the [`change-input-current` notebook](https://github.com/pybamm-team/PyBaMM/blob/master/examples/notebooks/change-input-current.ipynb) for information on how to specify an input current - Parameter functions must now use pybamm functions instead of numpy functions (e.g. `pybamm.exp` instead of `numpy.exp`), as these are then used to construct the expression tree directly. Generally, pybamm syntax follows numpy syntax; please get in touch if a function you need is missing. diff --git a/README.md b/README.md index e302eeb08a..789f4eae29 100644 --- a/README.md +++ b/README.md @@ -29,6 +29,8 @@ hosted on [Read The Docs](readthedocs.io). A set of slides giving an overview of can be found [here](https://github.com/pybamm-team/pybamm_summary_slides/blob/master/pybamm.pdf). +For further examples, see the list of repositories that use PyBaMM [here](https://github.com/pybamm-team/pybamm-example-results) + ## How can I obtain & install PyBaMM? ### Linux diff --git a/examples/README.md b/examples/README.md new file mode 100644 index 0000000000..c717c5d98a --- /dev/null +++ b/examples/README.md @@ -0,0 +1,4 @@ +# Examples + +A collection of Python scripts and Jupyter notebooks that demonstrate how to use PyBaMM. +For further examples, see the list of repositories that use PyBaMM [here](https://github.com/pybamm-team/pybamm-example-results) diff --git a/results/2019_08_sulzer_thesis/README.md b/results/2019_08_sulzer_thesis/README.md deleted file mode 100644 index b0cb9bac69..0000000000 --- a/results/2019_08_sulzer_thesis/README.md +++ /dev/null @@ -1,48 +0,0 @@ -# Results: Valentin Sulzer's Thesis - -This folder contains the scripts used to generate results for Valentin Sulzer's thesis -The plots were formatted using a formatting file `matplotlibrc` identical to [this one](_matplotlibrc) (but not included in the GitHub repo to avoid clashes with different formatting files). - -## Chapter 2 - Model - -- (Dimensionless and dimensional) [parameters](print_lead_acid_parameters.py) - - Function to print out all standard parameters to a text file, including dependence on C-rate for dimensionless parameters - -## Chapter 3 - Simplified models for slow discharge - -- [Effect of capacitance](effect_of_capacitance.py): comparison of the one-dimensional model with and without capacitance terms included - - Voltages - - Time taken for full solution -- [Discharge asymptotics results](lead_acid_discharge.py): - - Comparison of voltage curves for the full porous electrode model and a hierarchy of simplified models (Leading-Order Quasi-Static, First-Order Quasi-Static and Composite) - - Comparison of variable profiles at various times for the full and reduced-order models - - Electrolyte concentration - - Electrolyte potential - - Interfacial current density - - Errors compared to full model and time taken to solve each model - - Decomposition of voltage into constituent overpotentials -- [Effect of convection](effect_of_convection.py): - - Voltage at various C-rates with and without convection - - Velocity profiles - - Increasing the volume changes to see more of an effect -- [Effect of side reactions](effect_of_side_reactions.py): - - Voltage at various C-rates with and without side reactions - - Interfacial current densities - - Electrolyte concentrations -- [Charge asymptotics results](lead_acid_charge.py): - - Comparison of voltage curves for the full porous electrode model and a hierarchy of simplified models (Leading-Order Quasi-Static, First-Order Quasi-Static and Composite) - - Comparison of average interfacial current densities for each of the side reactions for the full porous electrode model and a hierarchy of simplified models (Leading-Order Quasi-Static, First-Order Quasi-Static and Composite) - - Comparison of variable profiles at various times for the full and reduced-order models - - Electrolyte concentration - - Oxygen concentration - - Decomposition of voltage into constituent overpotentials -- [Self-discharge](self_discharge.py): - - Self-discharge voltages - -## Chapter 4 - Small aspect ratio cells - -- 2+1D model - - Model and capacitance formulation - - Concentrations and potentials as functions of x, y, z - - Times taken -- Further asymptotics diff --git a/results/2019_08_sulzer_thesis/_matplotlibrc b/results/2019_08_sulzer_thesis/_matplotlibrc deleted file mode 100644 index 7e8738810d..0000000000 --- a/results/2019_08_sulzer_thesis/_matplotlibrc +++ /dev/null @@ -1,27 +0,0 @@ -# Set a 2/1 aspect ratio with automatic layout. -figure.figsize: 6.4, 4. - -# Use LaTeX fonts -font.family: serif -font.serif: ComputerModern -text.usetex: true - -# Backend options (42: use TrueType fonts) -pdf.fonttype: 42 -ps.fonttype: 42 - -# Set large font sizes for paper figures -axes.labelsize: 11 -axes.titlesize: 11 -xtick.labelsize: 11 -ytick.labelsize: 11 -legend.fontsize: 11 -legend.numpoints: 1 -legend.scatterpoints: 1 - -# Set lines and ticks according to the `seaborn-paper` style sheet -grid.linewidth: 0.8 -lines.linewidth: 1.4 -patch.linewidth: 0.24 -lines.markersize: 5.6 -lines.markeredgewidth: 0 diff --git a/results/2019_08_sulzer_thesis/effect_of_capacitance.py b/results/2019_08_sulzer_thesis/effect_of_capacitance.py deleted file mode 100644 index c182b928a9..0000000000 --- a/results/2019_08_sulzer_thesis/effect_of_capacitance.py +++ /dev/null @@ -1,159 +0,0 @@ -# -# Simulations: discharge of a lead-acid battery -# -import argparse -import matplotlib.pyplot as plt -import numpy as np -import pickle -import pybamm -import shared_plotting -from config import OUTPUT_DIR -from mpl_toolkits.axes_grid1.inset_locator import inset_axes -from shared_solutions import model_comparison, convergence_study - -save_folder = "results/2019_08_sulzer_thesis/data/capacitance_results/" - - -def plot_voltages(all_variables, t_eval): - linestyles = ["k-", "b-.", "r--"] - _, axes = shared_plotting.plot_voltages( - all_variables, t_eval, linestyles=linestyles, figsize=(6.4, 4) - ) - - # Add inset plot - for k, (Crate, models_variables) in enumerate(all_variables.items()): - ax = axes.flat[k] - y_min = ax.get_ylim()[0] - ax.set_ylim([y_min, 13.6]) - inset = inset_axes(ax, width="40%", height="30%", loc=1, borderpad=0) - for j, variables in enumerate(models_variables.values()): - time = variables["Time [s]"](t_eval) - capacitance_indices = np.where(time < 50) - time = time[capacitance_indices] - voltage = variables["Battery voltage [V]"](t_eval)[capacitance_indices] - inset.plot(time, voltage, linestyles[j]) - inset.set_xlabel("Time [s]", fontsize=9) - inset.set_xlim([0, 3]) - inset.tick_params(axis="both", which="major", labelsize=9) - - file_name = "capacitance_voltage_comparison.eps" - if OUTPUT_DIR is not None: - plt.savefig(OUTPUT_DIR + file_name, format="eps", dpi=1000) - - -def plot_errors(all_variables, t_eval, Crates): - # Linestyles - linestyles = ["k-", "b-.", "r--"] - # Only use some Crates - all_variables = {k: v for k, v in all_variables.items() if k in Crates} - # Plot - fig, ax = plt.subplots(1, 1, figsize=(6.4, 4)) - for k, (Crate, models_variables) in enumerate(all_variables.items()): - ax.set_xlabel("Time [h]") - ax.set_ylabel("Error [V]") - - for j, (model, variables) in enumerate(models_variables.items()): - if model == "direct form": - base_model_results = models_variables[model] - continue - error = np.abs( - variables["Battery voltage [V]"](t_eval) - - base_model_results["Battery voltage [V]"](t_eval) - ) - ax.loglog(variables["Time [h]"](t_eval), error, linestyles[j], label=model) - ax.legend(loc="upper right") - fig.tight_layout() - file_name = "capacitance_errors_voltages.eps".format(Crate) - if OUTPUT_DIR is not None: - plt.savefig(OUTPUT_DIR + file_name, format="eps", dpi=1000) - - -def discharge_states(compute): - savefile = "effect_of_capacitance_data.pickle" - if compute: - models = [ - pybamm.lead_acid.Full(name="direct form"), - pybamm.lead_acid.Full( - {"surface form": "differential"}, - name="capacitance form\n(differential)", - ), - pybamm.lead_acid.Full( - {"surface form": "algebraic"}, name="capacitance form\n(algebraic)" - ), - ] - Crates = [0.1, 1, 5] - t_eval = np.concatenate( - [np.logspace(-6, -3, 50), np.linspace(0.001, 1, 100)[1:]] - ) - all_variables, t_eval = model_comparison(models, Crates, t_eval) - with open(savefile, "wb") as f: - data = (all_variables, t_eval) - pickle.dump(data, f, pickle.HIGHEST_PROTOCOL) - else: - try: - with open(savefile, "rb") as f: - (all_variables, t_eval) = pickle.load(f) - except FileNotFoundError: - raise FileNotFoundError( - "Run script with '--compute' first to generate results" - ) - plot_voltages(all_variables, t_eval) - plot_errors(all_variables, t_eval, [5]) - - -def plot_times(models_times_and_voltages): - shared_plotting.plot_times( - models_times_and_voltages, Crate=1, linestyles=["k-", "b-.", "r--"] - ) - file_name = "capacitance_solver_times.eps" - if OUTPUT_DIR is not None: - plt.savefig(OUTPUT_DIR + file_name, format="eps", dpi=1000) - - -def discharge_times_and_errors(compute): - savefile = "capacitance_times_and_errors.pickle" - if compute: - try: - with open(savefile, "rb") as f: - models_times_and_voltages = pickle.load(f) - except FileNotFoundError: - models_times_and_voltages = pybamm.get_infinite_nested_dict() - models = [ - pybamm.lead_acid.Full(name="direct form"), - pybamm.lead_acid.Full( - {"surface form": "differential"}, - name="capacitance form\n(differential)", - ), - pybamm.lead_acid.Full( - {"surface form": "algebraic"}, name="capacitance form\n(algebraic)" - ), - ] - Crates = [1] - all_npts = np.linspace(10, 100, 2) - t_eval = np.linspace(0, 1, 100) - new_models_times_and_voltages = convergence_study( - models, Crates, all_npts, t_eval - ) - models_times_and_voltages.update(new_models_times_and_voltages) - with open(savefile, "wb") as f: - pickle.dump(models_times_and_voltages, f, pickle.HIGHEST_PROTOCOL) - else: - try: - with open(savefile, "rb") as f: - models_times_and_voltages = pickle.load(f) - except FileNotFoundError: - raise FileNotFoundError( - "Run script with '--compute' first to generate results" - ) - plot_errors(models_times_and_voltages) - # plot_times(models_times_and_voltages) - - -if __name__ == "__main__": - pybamm.set_logging_level("INFO") - parser = argparse.ArgumentParser() - parser.add_argument("--compute", action="store_true", help="(Re)-compute results.") - args = parser.parse_args() - discharge_states(args.compute) - # discharge_times_and_errors(args.compute) - plt.show() diff --git a/results/2019_08_sulzer_thesis/effect_of_convection.py b/results/2019_08_sulzer_thesis/effect_of_convection.py deleted file mode 100644 index ab1a3000e9..0000000000 --- a/results/2019_08_sulzer_thesis/effect_of_convection.py +++ /dev/null @@ -1,125 +0,0 @@ -# -# Simulations: effect of side reactions for charge of a lead-acid battery -# -import argparse -import matplotlib.pyplot as plt -import numpy as np -import pickle -import pybamm -import shared_plotting -from shared_solutions import model_comparison - -try: - from config import OUTPUT_DIR -except ImportError: - OUTPUT_DIR = None - - -def plot_voltages(all_variables, t_eval, bigger_beta=False): - linestyles = ["k-", "b--"] - shared_plotting.plot_voltages(all_variables, t_eval, linestyles, figsize=(6.4, 2.5)) - if bigger_beta: - file_name = "convection_voltage_comparison_bigger_beta.eps" - else: - file_name = "convection_voltage_comparison.eps" - plt.subplots_adjust(bottom=0.4) - if OUTPUT_DIR is not None: - plt.savefig(OUTPUT_DIR + file_name, format="eps", dpi=1000) - - -def plot_variables(all_variables, t_eval, bigger_beta=False): - # Set up - times = np.array([0.195]) - linestyles = ["k-", "b--"] - if bigger_beta: - var_file_names = { - "Volume-averaged velocity [m.s-1]" - + "": "convection_velocity_comparison_bigger_beta.eps", - "Electrolyte concentration [Molar]" - + "": "convection_electrolyte_concentration_comparison_bigger_beta.eps", - } - else: - var_file_names = { - "Volume-averaged velocity [m.s-1]": "convection_velocity_comparison.eps", - "Electrolyte concentration [Molar]" - + "": "convection_electrolyte_concentration_comparison.eps", - } - for var, file_name in var_file_names.items(): - fig, axes = shared_plotting.plot_variable( - all_variables, times, var, linestyles=linestyles, figsize=(6.4, 3) - ) - for ax in axes.flat: - title = ax.get_title() - ax.set_title(title, y=1.08) - plt.subplots_adjust( - bottom=0.3, top=0.85, left=0.1, right=0.9, hspace=0.08, wspace=0.05 - ) - if OUTPUT_DIR is not None: - plt.savefig(OUTPUT_DIR + file_name, format="eps", dpi=1000) - - -def charge_states(compute): - savefile = "effect_of_convection_data.pickle" - if compute: - models = [ - pybamm.lead_acid.Full( - {"convection": True}, name="With convection" - ), - pybamm.lead_acid.Full(name="Without convection"), - ] - Crates = [0.5, 1, 5] - t_eval = np.linspace(0, 1, 100) - all_variables, t_eval = model_comparison(models, Crates, t_eval) - with open(savefile, "wb") as f: - data = (all_variables, t_eval) - pickle.dump(data, f, pickle.HIGHEST_PROTOCOL) - else: - try: - with open(savefile, "rb") as f: - (all_variables, t_eval) = pickle.load(f) - except FileNotFoundError: - raise FileNotFoundError( - "Run script with '--compute' first to generate results" - ) - plot_voltages(all_variables, t_eval) - plot_variables(all_variables, t_eval) - - -def charge_states_bigger_volume_change(compute): - savefile = "effect_of_convection_bigger_beta_data.pickle" - if compute: - models = [ - pybamm.lead_acid.Full( - {"convection": True}, name="With convection" - ), - pybamm.lead_acid.Full(name="Without convection"), - ] - Crates = [0.5, 1, 5] - t_eval = np.linspace(0, 1, 100) - extra_parameter_values = {"Volume change factor": 10} - all_variables, t_eval = model_comparison( - models, Crates, t_eval, extra_parameter_values=extra_parameter_values - ) - with open(savefile, "wb") as f: - data = (all_variables, t_eval) - pickle.dump(data, f, pickle.HIGHEST_PROTOCOL) - else: - try: - with open(savefile, "rb") as f: - (all_variables, t_eval) = pickle.load(f) - except FileNotFoundError: - raise FileNotFoundError( - "Run script with '--compute' first to generate results" - ) - plot_voltages(all_variables, t_eval, bigger_beta=True) - plot_variables(all_variables, t_eval, bigger_beta=True) - - -if __name__ == "__main__": - pybamm.set_logging_level("DEBUG") - parser = argparse.ArgumentParser() - parser.add_argument("--compute", action="store_true", help="(Re)-compute results.") - args = parser.parse_args() - charge_states(args.compute) - charge_states_bigger_volume_change(args.compute) - plt.show() diff --git a/results/2019_08_sulzer_thesis/effect_of_side_reactions.py b/results/2019_08_sulzer_thesis/effect_of_side_reactions.py deleted file mode 100644 index 7c7cdbbc1a..0000000000 --- a/results/2019_08_sulzer_thesis/effect_of_side_reactions.py +++ /dev/null @@ -1,111 +0,0 @@ -# -# Simulations: effect of side reactions for charge of a lead-acid battery -# -import argparse -import matplotlib.pyplot as plt -import numpy as np -import pickle -import pybamm -import shared_plotting -from shared_solutions import model_comparison - -try: - from config import OUTPUT_DIR -except ImportError: - OUTPUT_DIR = None - - -def plot_voltages(all_variables, t_eval): - linestyles = ["k-", "b--"] - shared_plotting.plot_voltages(all_variables, t_eval, linestyles, figsize=(6.4, 2.5)) - file_name = "side_reactions_voltage_comparison.eps" - plt.subplots_adjust(bottom=0.4) - if OUTPUT_DIR is not None: - plt.savefig(OUTPUT_DIR + file_name, format="eps", dpi=1000) - - -def plot_interfacial_currents(all_variables, t_eval): - file_name = "side_reactions_interfacial_current_density_comparison.eps" - output_vars = [ - "Average positive electrode interfacial current density", - "Average positive electrode oxygen interfacial current density", - "Average negative electrode oxygen interfacial current density", - "Average negative electrode interfacial current density", - ] - labels = [ - "Pos electrode\n(main)", - "Pos electrode\n(oxygen)", - "Neg electrode\n(oxygen)", - "Neg electrode\n(main)", - ] - shared_plotting.plot_time_dependent_variables( - all_variables, t_eval, output_vars, labels - ) - plt.subplots_adjust(bottom=0.4, right=0.95, wspace=0.3) - if OUTPUT_DIR is not None: - plt.savefig(OUTPUT_DIR + file_name, format="eps", dpi=1000) - - -def charge_states(compute): - savefile1 = "effect_of_side_reactions_data.pickle" - savefile2 = "effect_of_side_reactions_loqs_data.pickle" - if compute: - models1 = [ - pybamm.lead_acid.Full( - {"surface form": "algebraic", "side reactions": ["oxygen"]}, - name="With oxygen", - ), - pybamm.lead_acid.Full( - {"surface form": "algebraic"}, name="Without oxygen" - ), - ] - Crates = [-0.1, -1, -5] - t_eval = np.linspace(0, 5, 100) - extra_parameter_values = { - "Positive electrode" - + "reference exchange-current density (oxygen) [A.m-2]": 1e-24, - "Initial State of Charge": 0.5, - } - all_variables1, t_eval1 = model_comparison( - models1, Crates, t_eval, extra_parameter_values=extra_parameter_values - ) - # Use LOQS without voltage cut-off for interfacial current densities, so that - # the current goes all the way - models2 = [ - pybamm.lead_acid.Full( - {"surface form": "algebraic", "side reactions": ["oxygen"]}, - name="With oxygen", - ), - pybamm.lead_acid.LOQS({"surface form": "algebraic"}, name="Without oxygen"), - ] - extra_parameter_values["Upper voltage cut-off [V]"] = 100 - all_variables2, t_eval2 = model_comparison( - models2, Crates, t_eval, extra_parameter_values=extra_parameter_values - ) - with open(savefile1, "wb") as f: - data1 = (all_variables1, t_eval1) - pickle.dump(data1, f, pickle.HIGHEST_PROTOCOL) - with open(savefile2, "wb") as f: - data2 = (all_variables2, t_eval2) - pickle.dump(data2, f, pickle.HIGHEST_PROTOCOL) - else: - try: - with open(savefile1, "rb") as f: - (all_variables1, t_eval1) = pickle.load(f) - with open(savefile2, "rb") as f: - (all_variables2, t_eval2) = pickle.load(f) - except FileNotFoundError: - raise FileNotFoundError( - "Run script with '--compute' first to generate results" - ) - plot_voltages(all_variables1, t_eval1) - plot_interfacial_currents(all_variables2, t_eval2) - - -if __name__ == "__main__": - pybamm.set_logging_level("DEBUG") - parser = argparse.ArgumentParser() - parser.add_argument("--compute", action="store_true", help="(Re)-compute results.") - args = parser.parse_args() - charge_states(args.compute) - plt.show() diff --git a/results/2019_08_sulzer_thesis/lead_acid_charge.py b/results/2019_08_sulzer_thesis/lead_acid_charge.py deleted file mode 100644 index c4f59a9010..0000000000 --- a/results/2019_08_sulzer_thesis/lead_acid_charge.py +++ /dev/null @@ -1,145 +0,0 @@ -# -# Simulations: charge of a lead-acid battery -# -import argparse -import matplotlib.pyplot as plt -import numpy as np -import pickle -import pybamm -import shared_plotting -from shared_solutions import model_comparison - -try: - from config import OUTPUT_DIR -except ImportError: - OUTPUT_DIR = None - - -def plot_voltages(all_variables, t_eval): - Crates = [-0.1, -0.2, -0.5, -1, -2, -5] - all_variables = {k: v for k, v in all_variables.items() if k in Crates} - shared_plotting.plot_voltages(all_variables, t_eval) - file_name = "charge_voltage_comparison.eps" - if OUTPUT_DIR is not None: - plt.savefig(OUTPUT_DIR + file_name, format="eps", dpi=1000) - - -def plot_interfacial_currents(all_variables, t_eval): - Crates = [-0.1, -2, -5] - all_variables = {Crate: v for Crate, v in all_variables.items() if Crate in Crates} - file_name = "charge_interfacial_current_density_comparison.eps" - output_vars = [ - "Average positive electrode interfacial current density", - "Average positive electrode oxygen interfacial current density", - "Average negative electrode oxygen interfacial current density", - "Average negative electrode interfacial current density", - ] - labels = [ - "Pos electrode\n(main)", - "Pos electrode\n(oxygen)", - "Neg electrode\n(oxygen)", - "Neg electrode\n(main)", - ] - shared_plotting.plot_time_dependent_variables( - all_variables, - t_eval, - output_vars, - labels, - colors=["k", "g", "r", "b"], - figsize=(6.4, 6.4), - ) - plt.subplots_adjust( - bottom=0.15, left=0.15, right=0.95, wspace=0.3, hspace=0.4, top=0.95 - ) - if OUTPUT_DIR is not None: - plt.savefig(OUTPUT_DIR + file_name, format="eps", dpi=1000) - - -def plot_variables(all_variables, t_eval): - # Set up - Crates = [-0.1, -2, -5] - times = np.linspace(0, 2, 4) - var_file_names = { - "Electrolyte concentration [Molar]" - + "": "charge_electrolyte_concentration_comparison.eps", - "Oxygen concentration [Molar]": "charge_oxygen_concentration_comparison.eps", - } - limits_exceptions = {"Electrolyte concentration [Molar]": {"min": 0}} - all_variables = {k: v for k, v in all_variables.items() if k in Crates} - for var, file_name in var_file_names.items(): - if var in limits_exceptions: - exceptions = limits_exceptions[var] - else: - exceptions = {} - shared_plotting.plot_variable( - all_variables, times, var, exceptions, yaxis="FCI" - ) - if OUTPUT_DIR is not None: - plt.savefig(OUTPUT_DIR + file_name, format="eps", dpi=1000) - - -def plot_voltage_components(all_variables, t_eval): - Crates = [-0.1, -2, -5] - model = "Composite" - shared_plotting.plot_voltage_components(all_variables, t_eval, model, Crates) - file_name = "charge_voltage_components.eps" - if OUTPUT_DIR is not None: - plt.savefig(OUTPUT_DIR + file_name, format="eps", dpi=1000) - - -def charge_states(compute): - if compute: - models = [ - pybamm.lead_acid.Full( - {"side reactions": ["oxygen"]}, name="Full" - ), - pybamm.lead_acid.LOQS( - {"surface form": "algebraic", "side reactions": ["oxygen"]}, name="LOQS" - ), - pybamm.lead_acid.FOQS( - {"surface form": "algebraic", "side reactions": ["oxygen"]}, name="FOQS" - ), - # pybamm.lead_acid.Composite( - # {"surface form": "algebraic", "side reactions": ["oxygen"]}, - # name="Composite", - # ), - pybamm.lead_acid.CompositeExtended( - {"surface form": "algebraic", "side reactions": ["oxygen"]}, - name="Composite", - ), - ] - Crates = [-0.1, -0.2, -0.5, -1, -2, -5] - t_eval = np.linspace(0, 3, 100) - extra_parameter_values = { - "Positive electrode" - + "reference exchange-current density (oxygen) [A.m-2]": 1e-24, - "Initial State of Charge": 0.5, - } - all_variables, t_eval = model_comparison( - models, Crates, t_eval, extra_parameter_values=extra_parameter_values - ) - with open("charge_asymptotics_data.pickle", "wb") as f: - data = (all_variables, t_eval) - pickle.dump(data, f, pickle.HIGHEST_PROTOCOL) - else: - try: - with open("charge_asymptotics_data.pickle", "rb") as f: - (all_variables, t_eval) = pickle.load(f) - except FileNotFoundError: - raise FileNotFoundError( - "Run script with '--compute' first to generate results" - ) - plot_voltages(all_variables, t_eval) - plot_interfacial_currents(all_variables, t_eval) - plot_variables(all_variables, t_eval) - plot_voltage_components(all_variables, t_eval) - - -if __name__ == "__main__": - pybamm.set_logging_level("INFO") - parser = argparse.ArgumentParser() - parser.add_argument("--compute", action="store_true", help="(Re)-compute results.") - args = parser.parse_args() - charge_states(args.compute) - # charge_times_and_errors(args.compute) - plt.show() diff --git a/results/2019_08_sulzer_thesis/lead_acid_discharge.py b/results/2019_08_sulzer_thesis/lead_acid_discharge.py deleted file mode 100644 index 887521e0a4..0000000000 --- a/results/2019_08_sulzer_thesis/lead_acid_discharge.py +++ /dev/null @@ -1,169 +0,0 @@ -# -# Simulations: discharge of a lead-acid battery -# -import argparse -import matplotlib.pyplot as plt -import numpy as np -import pickle -import pybamm -import shared_plotting -from collections import defaultdict -from shared_solutions import model_comparison, convergence_study - -try: - from config import OUTPUT_DIR -except ImportError: - OUTPUT_DIR = None - - -def plot_voltages(all_variables, t_eval): - Crates = [0.1, 0.2, 0.5, 1, 2, 5] - all_variables = {k: v for k, v in all_variables.items() if k in Crates} - shared_plotting.plot_voltages(all_variables, t_eval) - file_name = "discharge_voltage_comparison.eps" - if OUTPUT_DIR is not None: - plt.savefig(OUTPUT_DIR + file_name, format="eps", dpi=1000) - - -def plot_variables(all_variables, t_eval): - # Set up - Crates = [0.1, 1, 4] - times = np.array([0, 0.195, 0.375, 0.545]) - var_file_names = { - "Electrolyte concentration [Molar]" - + "": "discharge_electrolyte_concentration_comparison.eps", - "Electrolyte potential [V]": "discharge_electrolyte_potential_comparison.eps", - "Interfacial current density" - + "": "discharge_interfacial_current_density_comparison.eps", - } - limits_exceptions = {"Electrolyte concentration [Molar]": {"min": 0}} - all_variables = {k: v for k, v in all_variables.items() if k in Crates} - for var, file_name in var_file_names.items(): - if var in limits_exceptions: - exceptions = limits_exceptions[var] - else: - exceptions = {} - shared_plotting.plot_variable(all_variables, times, var, exceptions) - if OUTPUT_DIR is not None: - plt.savefig(OUTPUT_DIR + file_name, format="eps", dpi=1000) - - -def plot_voltage_components(all_variables, t_eval): - Crates = [0.1, 2, 5] - model = "Composite" - shared_plotting.plot_voltage_components(all_variables, t_eval, model, Crates) - file_name = "discharge_voltage_components.eps" - if OUTPUT_DIR is not None: - plt.savefig(OUTPUT_DIR + file_name, format="eps", dpi=1000) - - -def discharge_states(compute): - savefile = "discharge_asymptotics_data.pickle" - if compute: - models = [ - pybamm.lead_acid.Full(name="Full"), - pybamm.lead_acid.LOQS(name="LOQS"), - pybamm.lead_acid.FOQS(name="FOQS"), - pybamm.lead_acid.Composite(name="Composite"), - ] - Crates = [0.1, 0.2, 0.5, 1, 2, 4, 5, 10, 20] - t_eval = np.linspace(0, 1, 100) - extra_parameter_values = {"Bruggeman coefficient": 0.001} - all_variables, t_eval = model_comparison( - models, Crates, t_eval, extra_parameter_values=extra_parameter_values - ) - with open(savefile, "wb") as f: - data = (all_variables, t_eval) - pickle.dump(data, f, pickle.HIGHEST_PROTOCOL) - else: - try: - with open(savefile, "rb") as f: - (all_variables, t_eval) = pickle.load(f) - except FileNotFoundError: - raise FileNotFoundError( - "Run script with '--compute' first to generate results" - ) - plot_voltages(all_variables, t_eval) - plot_variables(all_variables, t_eval) - plot_voltage_components(all_variables, t_eval) - - -def plot_errors(models_times_and_voltages): - npts = 20 - linestyles = ["k-", "g--", "r:", "b-."] - Crates = defaultdict(list) - voltage_errors = defaultdict(list) - fig, ax = plt.subplots(1, 1) - for i, (model, times_and_voltages) in enumerate(models_times_and_voltages.items()): - if model != "Full": - for Crate, variables in times_and_voltages[npts].items(): - Crates[model].append(Crate) - full_voltage = models_times_and_voltages["Full"][npts][Crate][ - "Battery voltage [V]" - ] - reduced_voltage = variables["Battery voltage [V]"] - voltage_errors[model].append(pybamm.rmse(full_voltage, reduced_voltage)) - ax.semilogx( - Crates[model], voltage_errors[model], linestyles[i], label=model - ) - ax.set_xlabel("C-rate") - ax.set_ylabel("RMSE [V]") - ax.legend(loc="best") - fig.tight_layout() - file_name = "discharge_asymptotics_rmse.eps" - if OUTPUT_DIR is not None: - plt.savefig(OUTPUT_DIR + file_name, format="eps", dpi=1000) - - -def plot_times(models_times_and_voltages): - shared_plotting.plot_times(models_times_and_voltages, Crate=1) - file_name = "discharge_asymptotics_solver_times.eps" - if OUTPUT_DIR is not None: - plt.savefig(OUTPUT_DIR + file_name, format="eps", dpi=1000) - - -def discharge_times_and_errors(compute): - savefile = "discharge_asymptotics_times_and_errors.pickle" - if compute: - try: - with open(savefile, "rb") as f: - models_times_and_voltages = pickle.load(f) - except FileNotFoundError: - models_times_and_voltages = pybamm.get_infinite_nested_dict() - models = [ - pybamm.lead_acid.Full( - {"surface form": "algebraic"}, name="Full" - ), - pybamm.lead_acid.LOQS(name="LOQS"), - # pybamm.lead_acid.FOQS(name="FOQS"), - # pybamm.lead_acid.Composite(name="Composite"), - ] - Crates = np.linspace(0.01, 5, 2) - all_npts = [20] - t_eval = np.linspace(0, 1, 100) - new_models_times_and_voltages = convergence_study( - models, Crates, all_npts, t_eval - ) - models_times_and_voltages.update(new_models_times_and_voltages) - with open(savefile, "wb") as f: - pickle.dump(models_times_and_voltages, f, pickle.HIGHEST_PROTOCOL) - else: - try: - with open(savefile, "rb") as f: - models_times_and_voltages = pickle.load(f) - except FileNotFoundError: - raise FileNotFoundError( - "Run script with '--compute' first to generate results" - ) - plot_errors(models_times_and_voltages) - plot_times(models_times_and_voltages) - - -if __name__ == "__main__": - pybamm.set_logging_level("INFO") - parser = argparse.ArgumentParser() - parser.add_argument("--compute", action="store_true", help="(Re)-compute results.") - args = parser.parse_args() - discharge_states(args.compute) - # discharge_times_and_errors(args.compute) - plt.show() diff --git a/results/2019_08_sulzer_thesis/print_lead_acid_parameters.py b/results/2019_08_sulzer_thesis/print_lead_acid_parameters.py deleted file mode 100644 index be13842e3a..0000000000 --- a/results/2019_08_sulzer_thesis/print_lead_acid_parameters.py +++ /dev/null @@ -1,10 +0,0 @@ -# -# Print parameters for lead-acid models -# -import pybamm - -parameters = pybamm.standard_parameters_lead_acid -parameter_values = pybamm.lead_acid.BaseModel().default_parameter_values -output_file = "results/2019_08_sulzer_thesis/parameters.txt" - -pybamm.print_parameters(parameters, parameter_values, output_file) diff --git a/results/2019_08_sulzer_thesis/self_discharge.py b/results/2019_08_sulzer_thesis/self_discharge.py deleted file mode 100644 index bee360be25..0000000000 --- a/results/2019_08_sulzer_thesis/self_discharge.py +++ /dev/null @@ -1,61 +0,0 @@ -# -# Simulations: self-discharge -# -import argparse -import matplotlib.pyplot as plt -import numpy as np -import pickle -import pybamm -import shared_plotting -from config import OUTPUT_DIR -from shared_solutions import model_comparison - - -def plot_voltages(all_variables, t_eval): - shared_plotting.plot_voltages(all_variables, t_eval) - file_name = "sefl_discharge_voltage_comparison.eps" - if OUTPUT_DIR is not None: - plt.savefig(OUTPUT_DIR + file_name, format="eps", dpi=1000) - - -def self_discharge_states(compute): - save_file = "self_discharge_data.pickle" - if compute: - models = [ - pybamm.lead_acid.Full(name="Full, without oxygen"), - pybamm.lead_acid.Full( - {"side reactions": ["oxygen"]}, name="Full, with oxygen" - ), - pybamm.lead_acid.LOQS( - {"surface form": "algebraic", "side reactions": ["oxygen"]}, - name="LOQS, with oxygen", - ), - ] - extra_parameter_values = { - "Current function": "[zero]" - } - t_eval = np.linspace(0, 1000, 100) - all_variables, t_eval = model_comparison( - models, [1], t_eval, extra_parameter_values=extra_parameter_values - ) - with open(save_file, "wb") as f: - data = (all_variables, t_eval) - pickle.dump(data, f, pickle.HIGHEST_PROTOCOL) - else: - try: - with open(save_file, "rb") as f: - (all_variables, t_eval) = pickle.load(f) - except FileNotFoundError: - raise FileNotFoundError( - "Run script with '--compute' first to generate results" - ) - plot_voltages(all_variables, t_eval) - - -if __name__ == "__main__": - pybamm.set_logging_level("INFO") - parser = argparse.ArgumentParser() - parser.add_argument("--compute", action="store_true", help="(Re)-compute results.") - args = parser.parse_args() - self_discharge_states(args.compute) - plt.show() diff --git a/results/2019_08_sulzer_thesis/shared_plotting.py b/results/2019_08_sulzer_thesis/shared_plotting.py deleted file mode 100644 index 313aa4a185..0000000000 --- a/results/2019_08_sulzer_thesis/shared_plotting.py +++ /dev/null @@ -1,335 +0,0 @@ -# -# Shared plotting -# -import matplotlib.pyplot as plt -import numpy as np -import pybamm -from collections import defaultdict - - -def plot_voltages(all_variables, t_eval, linestyles=None, figsize=(6.4, 4.5)): - # Plot - linestyles = linestyles or ["k-", "g--", "r:", "b-."] - n = int(len(all_variables) // np.sqrt(len(all_variables))) - m = int(np.ceil(len(all_variables) / n)) - fig, axes = plt.subplots(n, m, figsize=figsize) - labels = [model for model in [x for x in all_variables.values()][0].keys()] - y_min = 0.98 * min( - np.nanmin(variables["Battery voltage [V]"](t_eval)) - for models_variables in all_variables.values() - for variables in models_variables.values() - ) - y_max = 1.02 * max( - np.nanmax(variables["Battery voltage [V]"](t_eval)) - for models_variables in all_variables.values() - for variables in models_variables.values() - ) - # Strict voltage cut-offs - y_min = max(y_min, 10.5) - y_max = min(y_max, 14.6) - for k, (Crate, models_variables) in enumerate(all_variables.items()): - if len(all_variables) == 1: - ax = axes - else: - ax = axes.flat[k] - t_max = max( - np.nanmax(var["Time [h]"](t_eval)) for var in models_variables.values() - ) - ax.set_xlim([0, t_max]) - ax.set_ylim([y_min, y_max]) - ax.set_xlabel("Time [h]") - if len(all_variables) > 1: - ax.set_title( - "\\textbf{{({})}} {}C ($\\mathcal{{C}}_e={}$)".format( - chr(97 + k), abs(Crate), abs(Crate) * 0.6 - ) - ) - # # Hide the right and top spines - # ax.spines["right"].set_visible(False) - # ax.spines["top"].set_visible(False) - # - # # Only show ticks on the left and bottom spines - # ax.yaxis.set_ticks_position("left") - # ax.xaxis.set_ticks_position("bottom") - ax.xaxis.set_major_locator(plt.MaxNLocator(3)) - if k % m == 0: - ax.set_ylabel("Voltage [V]") - # else: - # ax.set_yticklabels([]) - - for j, variables in enumerate(models_variables.values()): - ax.plot( - variables["Time [h]"](t_eval), - variables["Battery voltage [V]"](t_eval), - linestyles[j], - ) - if len(all_variables) == 1: - leg = ax.legend(labels, loc="best") - fig.tight_layout() - else: - leg = fig.legend(labels, loc="lower center", ncol=len(labels)) - plt.subplots_adjust(bottom=0.25, right=0.95, hspace=1.1, wspace=0.3) - leg.get_frame().set_edgecolor("k") - return fig, axes - - -def plot_variable( - all_variables, - times, - variable, - limits_exceptions=None, - yaxis="SOC", - linestyles=None, - figsize=(6.4, 5), -): - limits_exceptions = limits_exceptions or {} - linestyles = linestyles or ["k-", "g--", "r:", "b-."] - n = len(times) - m = len(all_variables) - Crates = list(all_variables.keys()) - labels = [model for model in [x for x in all_variables.values()][0].keys()] - x = all_variables[Crates[0]][labels[0]]["x"](0, np.linspace(0, 1))[:, 0] - x_dim = all_variables[Crates[0]][labels[0]]["x [m]"](0, np.linspace(0, 1))[:, 0] - - fig, axes = plt.subplots(n, m, figsize=figsize) - - # Default limits - y_min = pybamm.ax_min( - [ - np.nanmin(variables[variable](times, x)) - for Crate, models_variables in all_variables.items() - for variables in models_variables.values() - ] - ) - y_max = pybamm.ax_max( - [ - np.nanmax(variables[variable](times, x)) - for Crate, models_variables in all_variables.items() - for variables in models_variables.values() - ] - ) - # Exceptions - if "min" in limits_exceptions: - y_min = limits_exceptions["min"] - if "max" in limits_exceptions: - y_max = limits_exceptions["max"] - - # Plot - for i, (Crate, models_variables) in enumerate(all_variables.items()): - for j, time in enumerate(times): - if len(times) == 1: - ax = axes[i] - else: - ax = axes[j, i] - ax.set_xlim([x_dim[0], x_dim[-1]]) - ax.set_ylim([y_min, y_max]) - ax.yaxis.set_major_locator(plt.MaxNLocator(3)) - - # Title - if j == 0: - ax.set_title( - "\\textbf{{({})}} {}C ($\\mathcal{{C}}_e={}$)".format( - chr(97 + i), abs(Crate), abs(Crate) * 0.6 - ) - ) - # x-axis - if j == len(times) - 1: - ax.set_xlabel("x [m]") - else: - ax.set_xticklabels([]) - - # y-axis - if i == 0: - # If we only want to plot one time the y label is the variable - if len(times) == 1: - ax.set_ylabel(variable) - # Otherwise the y label is the time - else: - for variables in models_variables.values(): - try: - if yaxis == "SOC": - soc = variables["State of Charge"](time) - ax.set_ylabel( - "{}\% SoC".format(int(soc)), rotation=0, labelpad=30 - ) - elif yaxis == "FCI": - fci = variables["Fractional Charge Input"](time) - ax.set_ylabel( - "{}\% FCI".format(int(fci)), rotation=0, labelpad=30 - ) - ax.yaxis.get_label().set_verticalalignment("center") - except ValueError: - pass - else: - ax.set_yticklabels([]) - - # Plot - for j, variables in enumerate(models_variables.values()): - ax.plot(x_dim, variables[variable](time, x), linestyles[j]) - leg = fig.legend(labels, loc="lower center", ncol=len(labels), frameon=True) - leg.get_frame().set_edgecolor("k") - plt.subplots_adjust( - bottom=0.17, top=0.95, left=0.18, right=0.97, hspace=0.08, wspace=0.05 - ) - return fig, axes - - -def plot_time_dependent_variables( - all_variables, t_eval, output_vars, labels, colors=None, figsize=(6.4, 3.5) -): - models = list(list(all_variables.values())[0].keys()) - full_model = models[0] - fig, axes = plt.subplots(len(models) - 1, len(all_variables), figsize=figsize) - y_min = pybamm.ax_min( - [ - np.nanmin(variables[var](t_eval)) - for models_variables in all_variables.values() - for variables in models_variables.values() - for var in output_vars - ] - ) - y_max = pybamm.ax_max( - [ - np.nanmax(variables[var](t_eval)) - for models_variables in all_variables.values() - for variables in models_variables.values() - for var in output_vars - ] - ) - linestyles = ["--", ":", "-.", "-"] - colors = colors or ["k", "b"] - for i, (Crate, models_variables) in enumerate(all_variables.items()): - for j, model in enumerate(models[1:]): - full_variables = models_variables[full_model] - variables = models_variables[model] - if len(models) == 2: - ax = axes[i] - else: - ax = axes[j, i] - t_max = max( - np.nanmax(var["Time [h]"](t_eval)) for var in models_variables.values() - ) - ax.set_xlim([0, t_max]) - ax.set_ylim([y_min, y_max]) - if i == 0: - if len(models) == 2: - ax.set_ylabel("Interfacial current densities") - else: - ax.set_ylabel("{}\nvs Full".format(model), rotation=0, labelpad=30) - ax.yaxis.get_label().set_verticalalignment("center") - if j == 0 and len(all_variables) > 1: - ax.set_title( - "\\textbf{{({})}} {}C ($\\mathcal{{C}}_e={}$)".format( - chr(97 + i), abs(Crate), abs(Crate) * 0.6 - ) - ) - if j == len(models) - 2: - ax.set_xlabel("Time [h]") - plots = {} - for k, var in enumerate(output_vars): - plots[(full_model, k)], = ax.plot( - full_variables["Time [h]"](t_eval), - full_variables[var](t_eval), - linestyle=linestyles[k], - color=colors[0], - ) - for k, var in enumerate(output_vars): - plots[(model, k)], = ax.plot( - variables["Time [h]"](t_eval), - variables[var](t_eval), - linestyle=linestyles[k], - color=colors[j + 1], - ) - if len(models) == 2: - leg1 = fig.legend( - [plots[(model, len(linestyles) - 1)] for model in models], - models, - loc="lower center", - ncol=len(models), - bbox_to_anchor=(0.5, 0), - ) - fig.legend( - labels, loc="lower center", ncol=len(labels), bbox_to_anchor=(0.5, 0.1) - ) - fig.add_artist(leg1) - else: - fig.legend(labels, loc="lower center", ncol=len(labels)) - - -def plot_voltage_components(all_variables, t_eval, model, Crates): - n = int(len(Crates) // np.sqrt(len(Crates))) - m = int(np.ceil(len(Crates) / n)) - fig, axes = plt.subplots(n, m, figsize=(6.4, 2.3)) - labels = ["V", "$V_U$", "$V_k$", "$V_c$", "$V_o$"] - overpotentials = [ - "Average battery reaction overpotential [V]", - "Average battery concentration overpotential [V]", - "Average battery electrolyte ohmic losses [V]", - ] - y_min = 0.95 * min( - np.nanmin(models_variables[model]["Battery voltage [V]"](t_eval)) - for models_variables in all_variables.values() - ) - y_max = 1.05 * max( - np.nanmax(models_variables[model]["Battery voltage [V]"](t_eval)) - for models_variables in all_variables.values() - ) - for k, Crate in enumerate(Crates): - variables = all_variables[Crate][model] - ax = axes.flat[k] - - # Set up - t_max = np.nanmax(variables["Time [h]"](t_eval)) - ax.set_xlim([0, t_max]) - ax.set_ylim([y_min, y_max]) - ax.set_xlabel("Time [h]") - ax.set_title( - "\\textbf{{({})}} {}C ($\\mathcal{{C}}_e={}$)".format( - chr(97 + k), abs(Crate), abs(Crate) * 0.6 - ) - ) - ax.xaxis.set_major_locator(plt.MaxNLocator(3)) - if k % m == 0: - ax.set_ylabel("Voltage [V]") - - # Plot - # Initialise - # for lead-acid we multiply everything by 6 to - time = variables["Time [h]"](t_eval) - initial_ocv = variables["Average battery open circuit voltage [V]"](0) - ocv = variables["Average battery open circuit voltage [V]"](t_eval) - ax.fill_between(time, ocv, initial_ocv) - top = ocv - # Plot - for overpotential in overpotentials: - bottom = top + variables[overpotential](t_eval) - ax.fill_between(time, bottom, top) - top = bottom - ax.plot(time, variables["Battery voltage [V]"](t_eval), "k--") - leg = axes.flat[-1].legend( - labels, bbox_to_anchor=(1.05, 0.5), loc="center left", frameon=True - ) - leg.get_frame().set_edgecolor("k") - fig.tight_layout() - - -def plot_times(models_times_and_voltages, Crate=1, linestyles=None): - linestyles = linestyles or ["k-", "g--", "r:", "b-."] - all_npts = defaultdict(list) - solver_times = defaultdict(list) - fig, ax = plt.subplots(1, 1) - for i, (model, times_and_voltages) in enumerate(models_times_and_voltages.items()): - for npts in times_and_voltages.keys(): - try: - solver_time = times_and_voltages[npts][Crate][ - "solution object" - ].solve_time - except KeyError: - continue - all_npts[model].append(npts * 3) - solver_times[model].append(solver_time) - ax.loglog(all_npts[model], solver_times[model], linestyles[i], label=model) - ax.set_xlabel("Number of grid points") - ax.set_ylabel("Solver time [s]") - ax.legend(loc="best") - fig.tight_layout() diff --git a/results/2019_08_sulzer_thesis/shared_solutions.py b/results/2019_08_sulzer_thesis/shared_solutions.py deleted file mode 100644 index 26369949d4..0000000000 --- a/results/2019_08_sulzer_thesis/shared_solutions.py +++ /dev/null @@ -1,147 +0,0 @@ -# -# Simulations -# -import pybamm - - -def model_comparison(models, Crates, t_eval, extra_parameter_values=None): - " Solve models at a range of Crates " - # load parameter values and geometry - geometry = models[0].default_geometry - extra_parameter_values = extra_parameter_values or {} - param = models[0].default_parameter_values - param.update(extra_parameter_values) - - # Process parameters (same parameters for all models) - for model in models: - param.process_model(model) - param.process_geometry(geometry) - - # set mesh - var = pybamm.standard_spatial_vars - var_pts = {var.x_n: 20, var.x_s: 20, var.x_p: 20} - mesh = pybamm.Mesh(geometry, models[-1].default_submesh_types, var_pts) - - # discretise models - discs = {} - for model in models: - disc = pybamm.Discretisation(mesh, model.default_spatial_methods) - disc.process_model(model) - # Store discretisation - discs[model] = disc - - # solve model for range of Crates - all_variables = {} - for Crate in Crates: - all_variables[Crate] = {} - current = Crate * 17 - pybamm.logger.info("Setting typical current to {} A".format(current)) - param.update({"Typical current [A]": current}) - for model in models: - param.update_model(model, discs[model]) - solution = model.default_solver.solve(model, t_eval) - variables = pybamm.post_process_variables( - model.variables, solution.t, solution.y, mesh - ) - variables["solution"] = solution - all_variables[Crate][model.name] = variables - - return all_variables, t_eval - - -def convergence_study(models, Crates, all_npts, t_eval, extra_parameter_values=None): - " Solve models at a range of number of grid points " - # load parameter values and geometry - geometry = models[0].default_geometry - param = models[0].default_parameter_values - # Update parameters - extra_parameter_values = extra_parameter_values or {} - param.update(extra_parameter_values) - - # Process parameters (same parameters for all models) - for model in models: - param.process_model(model) - param.process_geometry(geometry) - - # set mesh - var = pybamm.standard_spatial_vars - - # solve model for range of Crates and npts - models_times_and_voltages = {model.name: {} for model in models} - for npts in all_npts: - pybamm.logger.info("Setting number of grid points to {}".format(npts)) - var_pts = {var.x_n: npts, var.x_s: npts, var.x_p: npts} - mesh = pybamm.Mesh(geometry, models[-1].default_submesh_types, var_pts) - - # discretise models, store discretised model and discretisation - models_disc = {} - discs = {} - for model in models: - disc = pybamm.Discretisation(mesh, model.default_spatial_methods) - models_times_and_voltages[model.name][npts] = {} - models_disc[model.name] = disc.process_model(model, inplace=False) - discs[model.name] = disc - - # Solve for a range of C-rates - for Crate in Crates: - current = Crate * 17 - pybamm.logger.info("Setting typical current to {} A".format(current)) - param.update({"Typical current [A]": current}) - for model in models: - model_disc = models_disc[model.name] - disc = discs[model.name] - param.update_model(model_disc, disc) - try: - solution = model.default_solver.solve(model_disc, t_eval) - except pybamm.SolverError: - pybamm.logger.error( - "Could not solve {!s} at {} A with {} points".format( - model.name, current, npts - ) - ) - continue - voltage = pybamm.ProcessedVariable( - model_disc.variables["Battery voltage [V]"], solution.t, solution.y - )(t_eval) - variables = { - "Battery voltage [V]": voltage, - "solution object": solution, - } - models_times_and_voltages[model.name][npts][Crate] = variables - - return models_times_and_voltages - - -def simulation(models, t_eval, extra_parameter_values=None, disc_only=False): - - # create geometry - geometry = models[-1].default_geometry - - # load parameter values and process models and geometry - param = models[0].default_parameter_values - extra_parameter_values = extra_parameter_values or {} - param.update(extra_parameter_values) - for model in models: - param.process_model(model) - param.process_geometry(geometry) - - # set mesh - var = pybamm.standard_spatial_vars - var_pts = {var.x_n: 25, var.x_s: 41, var.x_p: 34} - mesh = pybamm.Mesh(geometry, models[-1].default_submesh_types, var_pts) - - # discretise models - for model in models: - disc = pybamm.Discretisation(mesh, model.default_spatial_methods) - disc.process_model(model) - - if disc_only: - return model, mesh - - # solve model - solutions = [None] * len(models) - for i, model in enumerate(models): - solution = model.default_solver.solve(model, t_eval) - solutions[i] = solution - - return models, mesh, solutions diff --git a/results/2plus1D/README.md b/results/2plus1D/README.md deleted file mode 100644 index 1bdf963769..0000000000 --- a/results/2plus1D/README.md +++ /dev/null @@ -1,9 +0,0 @@ -# Results: "2plus1D" models - -This folder contains the scripts used to generate results for the "2plus1D" papers: - -1. Asymptotic Reduction of a Li-ion Pouch Cell Model: Part 1 - Scott Marquis, Robert Timms, Valentin Sulzer, Colin Please, S Jon Chapman - -2. Asymptotic Reduction of a Li-ion Pouch Cell Model: Part 2 - Robert Timms, Scott Marquis, Valentin Sulzer, Colin Please, S Jon Chapman diff --git a/results/2plus1D/_matplotlibrc b/results/2plus1D/_matplotlibrc deleted file mode 100644 index 7e8738810d..0000000000 --- a/results/2plus1D/_matplotlibrc +++ /dev/null @@ -1,27 +0,0 @@ -# Set a 2/1 aspect ratio with automatic layout. -figure.figsize: 6.4, 4. - -# Use LaTeX fonts -font.family: serif -font.serif: ComputerModern -text.usetex: true - -# Backend options (42: use TrueType fonts) -pdf.fonttype: 42 -ps.fonttype: 42 - -# Set large font sizes for paper figures -axes.labelsize: 11 -axes.titlesize: 11 -xtick.labelsize: 11 -ytick.labelsize: 11 -legend.fontsize: 11 -legend.numpoints: 1 -legend.scatterpoints: 1 - -# Set lines and ticks according to the `seaborn-paper` style sheet -grid.linewidth: 0.8 -lines.linewidth: 1.4 -patch.linewidth: 0.24 -lines.markersize: 5.6 -lines.markeredgewidth: 0 diff --git a/results/2plus1D/compare_lead_acid_1plus1D.py b/results/2plus1D/compare_lead_acid_1plus1D.py deleted file mode 100644 index a6b015437e..0000000000 --- a/results/2plus1D/compare_lead_acid_1plus1D.py +++ /dev/null @@ -1,73 +0,0 @@ -import pybamm -import numpy as np -import matplotlib.pyplot as plt -import sys - -# set logging level and increase recursion limit -pybamm.set_logging_level("INFO") -sys.setrecursionlimit(10000) - -# load models -models = [ - pybamm.lead_acid.Full(name="1D Full"), - pybamm.lead_acid.Composite(name="1D composite"), - pybamm.lead_acid.LOQS(name="1D LOQS"), - pybamm.lead_acid.Full( - {"current collector": "potential pair", "dimensionality": 1}, name="1+1D Full" - ), - pybamm.lead_acid.Composite( - {"current collector": "potential pair", "dimensionality": 1}, - name="1+1D composite", - ), - pybamm.lead_acid.LOQS( - {"current collector": "potential pair", "dimensionality": 1}, name="1+1D LOQS" - ), -] - -# load parameter values and process models -param = models[0].default_parameter_values -for model in models: - param.process_model(model) - -# process geometry and discretise models -meshes = [None] * len(models) -for i, model in enumerate(models): - geometry = model.default_geometry - param.process_geometry(geometry) - var = pybamm.standard_spatial_vars - var_pts = { - var.x_n: 5, - var.x_s: 5, - var.x_p: 5, - var.r_n: 5, - var.r_p: 5, - var.y: 5, - var.z: 5, - } - meshes[i] = pybamm.Mesh(geometry, model.default_submesh_types, var_pts) - disc = pybamm.Discretisation(meshes[i], model.default_spatial_methods) - disc.process_model(model) - -# solve models and process time and voltage for plotting on different meshes -solutions = [None] * len(models) -times = [None] * len(models) -voltages = [None] * len(models) -t_eval = np.linspace(0, 1, 1000) -for i, model in enumerate(models): - solution = model.default_solver.solve(model, t_eval) - solutions[i] = solution - times[i] = pybamm.ProcessedVariable( - model.variables["Time [h]"], solution.t, solution.y - ) - voltages[i] = pybamm.ProcessedVariable( - model.variables["Terminal voltage [V]"], solution.t, solution.y, mesh=meshes[i] - ) - -# plot terminal voltage -t = np.linspace(0, solution.t[-1], 100) -for i, model in enumerate(models): - plt.plot(times[i](t), voltages[i](t), lw=2, label=model.name) -plt.xlabel("Time [h]", fontsize=15) -plt.ylabel("Terminal voltage [V]", fontsize=15) -plt.legend(fontsize=15) -plt.show() diff --git a/results/2plus1D/compare_lead_acid_2plus1D.py b/results/2plus1D/compare_lead_acid_2plus1D.py deleted file mode 100644 index afd7fdc6de..0000000000 --- a/results/2plus1D/compare_lead_acid_2plus1D.py +++ /dev/null @@ -1,75 +0,0 @@ -import pybamm -import numpy as np -import matplotlib.pyplot as plt -import sys - -# set logging level and increase recursion limit -pybamm.set_logging_level("INFO") -sys.setrecursionlimit(10000) - -# load models -models = [ - pybamm.lead_acid.Full(name="1D Full"), - pybamm.lead_acid.Composite(name="1D composite"), - pybamm.lead_acid.LOQS(name="1D LOQS"), - pybamm.lead_acid.Full( - {"current collector": "potential pair", "dimensionality": 2}, name="2+1D Full" - ), - pybamm.lead_acid.Composite( - {"current collector": "potential pair", "dimensionality": 2}, - name="2+1D composite", - ), - pybamm.lead_acid.LOQS( - {"current collector": "potential pair", "dimensionality": 2}, name="2+1D LOQS" - ), -] - -# load parameter values and process models -param = models[0].default_parameter_values -for model in models: - param.process_model(model) - -# process geometry and discretise models -meshes = [None] * len(models) -for i, model in enumerate(models): - geometry = model.default_geometry - param.process_geometry(geometry) - var = pybamm.standard_spatial_vars - var_pts = { - var.x_n: 5, - var.x_s: 5, - var.x_p: 5, - var.r_n: 5, - var.r_p: 5, - var.y: 5, - var.z: 5, - } - meshes[i] = pybamm.Mesh(geometry, model.default_submesh_types, var_pts) - disc = pybamm.Discretisation(meshes[i], model.default_spatial_methods) - disc.process_model(model) - -# solve models and process time and voltage for plotting on different meshes -solutions = [None] * len(models) -times = [None] * len(models) -voltages = [None] * len(models) -t_eval = np.linspace(0, 1, 1000) -for i, model in enumerate(models): - if "2+1D" in model.name: - model.use_simplify = False # simplifying jacobian slow for large systems - solution = model.default_solver.solve(model, t_eval) - solutions[i] = solution - times[i] = pybamm.ProcessedVariable( - model.variables["Time [h]"], solution.t, solution.y - ) - voltages[i] = pybamm.ProcessedVariable( - model.variables["Terminal voltage [V]"], solution.t, solution.y, mesh=meshes[i] - ) - -# plot terminal voltage -t = np.linspace(0, solution.t[-1], 100) -for i, model in enumerate(models): - plt.plot(times[i](t), voltages[i](t), lw=2, label=model.name) -plt.xlabel("Time [h]", fontsize=15) -plt.ylabel("Terminal voltage [V]", fontsize=15) -plt.legend(fontsize=15) -plt.show() diff --git a/results/2plus1D/compare_lithium_ion_2plus1D.py b/results/2plus1D/compare_lithium_ion_2plus1D.py deleted file mode 100644 index dfe33bba10..0000000000 --- a/results/2plus1D/compare_lithium_ion_2plus1D.py +++ /dev/null @@ -1,96 +0,0 @@ -import pybamm -import numpy as np -import matplotlib.pyplot as plt -import sys - -# set logging level and increase recursion limit -pybamm.set_logging_level("INFO") -sys.setrecursionlimit(10000) - -# load models -models = [ - pybamm.lithium_ion.SPM(name="1D SPM"), - pybamm.lithium_ion.SPMe(name="1D SPMe"), - pybamm.lithium_ion.DFN(name="1D DFN"), - pybamm.lithium_ion.SPM( - {"current collector": "potential pair", "dimensionality": 2}, name="2+1D SPM" - ), - pybamm.lithium_ion.SPMe( - {"current collector": "potential pair", "dimensionality": 2}, name="2+1D SPMe" - ), - pybamm.lithium_ion.DFN( - {"current collector": "potential pair", "dimensionality": 2}, name="2+1D DFN" - ), -] - -# load parameter values -param = models[0].default_parameter_values -C_rate = 1 -param.update({"C-rate": C_rate}) -# make current collectors not so conductive, just for illustrative purposes -param.update( - { - "Negative current collector conductivity [S.m-1]": 5.96e6, - "Positive current collector conductivity [S.m-1]": 3.55e6, - } -) - -# process models -for model in models: - param.process_model(model) - -# process geometry and discretise models -meshes = [None] * len(models) -for i, model in enumerate(models): - geometry = model.default_geometry - param.process_geometry(geometry) - var = pybamm.standard_spatial_vars - var_pts = { - var.x_n: 5, - var.x_s: 5, - var.x_p: 5, - var.r_n: 5, - var.r_p: 5, - var.y: 5, - var.z: 5, - } - meshes[i] = pybamm.Mesh(geometry, model.default_submesh_types, var_pts) - disc = pybamm.Discretisation(meshes[i], model.default_spatial_methods) - disc.process_model(model) - -# solve models and process time and voltage for plotting on different meshes -solutions = [None] * len(models) -times = [None] * len(models) -voltages = [None] * len(models) -t_eval = np.linspace(0, 1, 1000) -for i, model in enumerate(models): - if "2+1D" in model.name: - model.use_simplify = False # simplifying jacobian slow for large systems - solution = model.default_solver.solve(model, t_eval) - solutions[i] = solution - times[i] = pybamm.ProcessedVariable( - model.variables["Time [h]"], solution.t, solution.y - ) - voltages[i] = pybamm.ProcessedVariable( - model.variables["Terminal voltage [V]"], solution.t, solution.y, mesh=meshes[i] - ) - -# plot terminal voltage -t = np.linspace(0, solution.t[-1], 100) -for i, model in enumerate(models): - plt.plot(times[i](t), voltages[i](t), label=model.name) -plt.xlabel("Time [h]") -plt.ylabel("Terminal voltage [V]") -plt.legend() -# add C-rate, delta, and alpha to title -delta = param.evaluate(pybamm.standard_parameters_lithium_ion.delta) -alpha = param.evaluate(pybamm.standard_parameters_lithium_ion.alpha) -plt.title( - r"C-rate = {:3d}, $\alpha$ = {:.6f} , $\delta$ = {:.6f}".format( - C_rate, alpha, delta - ) -) -# save and show -file_name = "discharge_curve_2plus1D_comparison.eps" -plt.savefig(file_name, format="eps", dpi=1000) -plt.show() diff --git a/results/2plus1D/compare_spmecc.py b/results/2plus1D/compare_spmecc.py deleted file mode 100644 index 2f4cd2f136..0000000000 --- a/results/2plus1D/compare_spmecc.py +++ /dev/null @@ -1,196 +0,0 @@ -import pybamm -import numpy as np -import matplotlib.pyplot as plt -import sys - -# set logging level and increase recursion limit -pybamm.set_logging_level("INFO") -sys.setrecursionlimit(10000) - -# load current collector and SPMe models -cc_model = pybamm.current_collector.EffectiveResistance2D() -spme_av = pybamm.lithium_ion.SPMe(name="Average SPMe") -spme = pybamm.lithium_ion.SPMe( - {"current collector": "potential pair", "dimensionality": 2}, name="2+1D SPMe" -) -models = {"Current collector": cc_model, "Average SPMe": spme_av, "2+1D SPMe": spme} - -# set parameters based on the spme -param = spme.default_parameter_values - -# set mesh -var = pybamm.standard_spatial_vars -var_pts = { - var.x_n: 5, - var.x_s: 5, - var.x_p: 5, - var.r_n: 5, - var.r_p: 5, - var.y: 5, - var.z: 5, -} - -# process model and geometry, and discretise -meshes = {} -for name, model in models.items(): - param.process_model(model) - geometry = model.default_geometry - param.process_geometry(geometry) - meshes[name] = pybamm.Mesh(geometry, model.default_submesh_types, var_pts) - disc = pybamm.Discretisation(meshes[name], model.default_spatial_methods) - disc.process_model(model) - -# solve models -- simulate one hour discharge -tau = param.process_symbol(pybamm.standard_parameters_lithium_ion.tau_discharge) -t_end = 3600 / tau.evaluate(0) -t_eval = np.linspace(0, t_end, 120) -solutions = {} -for name, model in models.items(): - if name == "Current collector": - solutions[name] = model.default_solver.solve(model) - else: - solutions[name] = model.default_solver.solve(model, t_eval) - -# plot terminal voltage -for name in ["Average SPMe", "2+1D SPMe"]: - t, y = solutions[name].t, solutions[name].y - model = models[name] - time = pybamm.ProcessedVariable(model.variables["Time [h]"], t, y)(t) - voltage = pybamm.ProcessedVariable( - model.variables["Terminal voltage [V]"], t, y, mesh=meshes[name] - )(t) - - # add current collector Ohmic losses to average SPMEe to get SPMeCC voltage - if model.name == "Average SPMe": - current = pybamm.ProcessedVariable(model.variables["Current [A]"], t, y)(t) - delta = param.evaluate(pybamm.standard_parameters_lithium_ion.delta) - R_cc = param.process_symbol( - cc_model.variables["Effective current collector resistance [Ohm]"] - ).evaluate( - t=solutions["Current collector"].t, y=solutions["Current collector"].y - )[ - 0 - ][ - 0 - ] - cc_ohmic_losses = -delta * current * R_cc - voltage = voltage + cc_ohmic_losses - - # plot - plt.plot(time, voltage, label=model.name) -plt.xlabel("Time [h]") -plt.ylabel("Terminal voltage [V]") -plt.legend() - - -# plot potentials in current collector - -# get processed potentials from SPMeCC -V_av = pybamm.ProcessedVariable( - spme_av.variables["Terminal voltage"], - solutions["Average SPMe"].t, - solutions["Average SPMe"].y, - mesh=meshes["Average SPMe"], -) -I_av = pybamm.ProcessedVariable( - spme_av.variables["Total current density"], - solutions["Average SPMe"].t, - solutions["Average SPMe"].y, - mesh=meshes["Average SPMe"], -) -potentials = cc_model.get_processed_potentials( - solutions["Current collector"], meshes["Current collector"], param, V_av, I_av -) -phi_s_cn_spmecc = potentials["Negative current collector potential [V]"] -phi_s_cp_spmecc = potentials["Positive current collector potential [V]"] - -# get processed potentials from 2+1D SPMe -phi_s_cn = pybamm.ProcessedVariable( - model.variables["Negative current collector potential [V]"], - solutions["2+1D SPMe"].t, - solutions["2+1D SPMe"].y, - mesh=meshes["2+1D SPMe"], -) -phi_s_cp = pybamm.ProcessedVariable( - model.variables["Positive current collector potential [V]"], - solutions["2+1D SPMe"].t, - solutions["2+1D SPMe"].y, - mesh=meshes["2+1D SPMe"], -) - -# make plot -l_y = phi_s_cp.y_sol[-1] -l_z = phi_s_cp.z_sol[-1] -y_plot = np.linspace(0, l_y, 21) -z_plot = np.linspace(0, l_z, 21) - - -def plot(t): - plt.subplots(figsize=(15, 8)) - plt.tight_layout() - plt.subplots_adjust(left=-0.1) - - # negative current collector potential - plt.subplot(221) - phi_s_cn_plot = plt.pcolormesh( - y_plot, - z_plot, - np.transpose(phi_s_cn(y=y_plot, z=z_plot, t=t)), - shading="gouraud", - ) - plt.axis([0, l_y, 0, l_z]) - plt.xlabel(r"$y$") - plt.ylabel(r"$z$") - plt.title(r"$\phi_{s,cn}$") - plt.set_cmap("cividis") - plt.colorbar(phi_s_cn_plot) - plt.subplot(222) - phi_s_cn_spmecc_plot = plt.pcolormesh( - y_plot, - z_plot, - np.transpose(phi_s_cn_spmecc(y=y_plot, z=z_plot, t=t)), - shading="gouraud", - ) - plt.axis([0, l_y, 0, l_z]) - plt.xlabel(r"$y$") - plt.ylabel(r"$z$") - plt.title(r"$\phi_{s,cn}$ SPMeCC") - plt.set_cmap("cividis") - plt.colorbar(phi_s_cn_spmecc_plot) - - # positive current collector potential - plt.subplot(223) - phi_s_cp_plot = plt.pcolormesh( - y_plot, - z_plot, - np.transpose(phi_s_cp(y=y_plot, z=z_plot, t=t)), - shading="gouraud", - ) - - plt.axis([0, l_y, 0, l_z]) - plt.xlabel(r"$y$") - plt.ylabel(r"$z$") - plt.title(r"$\phi_{s,cp}$") - plt.set_cmap("viridis") - plt.colorbar(phi_s_cp_plot) - plt.subplot(224) - phi_s_cp_spmecc_plot = plt.pcolormesh( - y_plot, - z_plot, - np.transpose(phi_s_cp_spmecc(y=y_plot, z=z_plot, t=t)), - shading="gouraud", - ) - plt.axis([0, l_y, 0, l_z]) - plt.xlabel(r"$y$") - plt.ylabel(r"$z$") - plt.title(r"$\phi_{s,cp}$ SPMeCC") - plt.set_cmap("viridis") - plt.colorbar(phi_s_cp_spmecc_plot) - - plt.subplots_adjust( - top=0.92, bottom=0.15, left=0.10, right=0.9, hspace=0.5, wspace=0.5 - ) - - -plot(solutions["2+1D SPMe"].t[-1] / 2) -plt.show() diff --git a/results/2plus1D/compare_thermal_lithium_ion_2plus1D.py b/results/2plus1D/compare_thermal_lithium_ion_2plus1D.py deleted file mode 100644 index c1ac708e98..0000000000 --- a/results/2plus1D/compare_thermal_lithium_ion_2plus1D.py +++ /dev/null @@ -1,133 +0,0 @@ -import pybamm -import numpy as np -import matplotlib.pyplot as plt -import sys - -# set logging level and increase recursion limit -pybamm.set_logging_level("INFO") -sys.setrecursionlimit(10000) - -# load models -models = [ - pybamm.lithium_ion.SPM({"thermal": "x-lumped"}, name="1D SPM (lumped)"), - pybamm.lithium_ion.SPMe({"thermal": "x-lumped"}, name="1D SPMe (lumped)"), - pybamm.lithium_ion.DFN({"thermal": "x-lumped"}, name="1D DFN (lumped)"), - pybamm.lithium_ion.SPM({"thermal": "x-full"}, name="1D SPM (full)"), - pybamm.lithium_ion.SPMe({"thermal": "x-full"}, name="1D SPMe (full)"), - pybamm.lithium_ion.DFN({"thermal": "x-full"}, name="1D DFN (full)"), - pybamm.lithium_ion.SPM( - { - "current collector": "potential pair", - "dimensionality": 2, - "thermal": "xyz-lumped", - }, - name="2+1D SPM (lumped)", - ), - pybamm.lithium_ion.SPMe( - { - "current collector": "potential pair", - "dimensionality": 2, - "thermal": "xyz-lumped", - }, - name="2+1D SPMe (lumped)", - ), - pybamm.lithium_ion.DFN( - { - "current collector": "potential pair", - "dimensionality": 2, - "thermal": "xyz-lumped", - }, - name="2+1D DFN (lumped)", - ), - pybamm.lithium_ion.SPM( - { - "current collector": "potential pair", - "dimensionality": 2, - "thermal": "x-lumped", - }, - name="2+1D SPM (full)", - ), - pybamm.lithium_ion.SPMe( - { - "current collector": "potential pair", - "dimensionality": 2, - "thermal": "x-lumped", - }, - name="2+1D SPMe (full)", - ), - pybamm.lithium_ion.DFN( - { - "current collector": "potential pair", - "dimensionality": 2, - "thermal": "x-lumped", - }, - name="2+1D DFN (full)", - ), -] - -# load parameter values -param = models[0].default_parameter_values - -# process models -for model in models: - param.process_model(model) - -# process geometry and discretise models -meshes = [None] * len(models) -for i, model in enumerate(models): - geometry = model.default_geometry - param.process_geometry(geometry) - var = pybamm.standard_spatial_vars - var_pts = { - var.x_n: 5, - var.x_s: 5, - var.x_p: 5, - var.r_n: 5, - var.r_p: 5, - var.y: 5, - var.z: 5, - } - meshes[i] = pybamm.Mesh(geometry, model.default_submesh_types, var_pts) - disc = pybamm.Discretisation(meshes[i], model.default_spatial_methods) - disc.process_model(model) - -# solve models and process time and voltage for plotting on different meshes -solutions = [None] * len(models) -times = [None] * len(models) -voltages = [None] * len(models) -temperatures = [None] * len(models) - -t_eval = np.linspace(0, 1, 1000) -for i, model in enumerate(models): - if "2+1D" in model.name: - model.use_simplify = False # simplifying jacobian slow for large systems - solution = model.default_solver.solve(model, t_eval) - solutions[i] = solution - times[i] = pybamm.ProcessedVariable( - model.variables["Time [h]"], solution.t, solution.y - ) - voltages[i] = pybamm.ProcessedVariable( - model.variables["Terminal voltage [V]"], solution.t, solution.y, mesh=meshes[i] - ) - temperatures[i] = pybamm.ProcessedVariable( - model.variables["Volume-averaged cell temperature [K]"], - solution.t, - solution.y, - mesh=meshes[i], - ) - -# plot terminal voltage and temperature -t = np.linspace(0, solution.t[-1], 100) -plt.subplot(121) -for i, model in enumerate(models): - plt.plot(times[i](t), voltages[i](t), label=model.name) -plt.xlabel("Time [h]") -plt.ylabel("Terminal voltage [V]") -plt.legend() -plt.subplot(122) -for i, model in enumerate(models): - plt.plot(times[i](t), temperatures[i](t), label=model.name) -plt.xlabel("Time [h]") -plt.ylabel("Temperature [K]") -plt.tight_layout() -plt.show() diff --git a/results/2plus1D/dfn_2plus1D.py b/results/2plus1D/dfn_2plus1D.py deleted file mode 100644 index 32daf9ed58..0000000000 --- a/results/2plus1D/dfn_2plus1D.py +++ /dev/null @@ -1,137 +0,0 @@ -import pybamm -import numpy as np -import matplotlib.pyplot as plt -import sys - -# set logging level -pybamm.set_logging_level("INFO") - -# load (2+1D) DFN model -options = { - "current collector": "potential pair", - "dimensionality": 2, - "thermal": "x-lumped", -} -model = pybamm.lithium_ion.DFN(options) -model.use_simplify = False # simplifying jacobian slow for large systems - -# create geometry -geometry = model.default_geometry - -# load parameter values and process model and geometry -param = model.default_parameter_values -param.process_model(model) -param.process_geometry(geometry) - -# set mesh -var = pybamm.standard_spatial_vars -var_pts = { - var.x_n: 5, - var.x_s: 5, - var.x_p: 5, - var.r_n: 5, - var.r_p: 5, - var.y: 5, - var.z: 5, -} -# depending on number of points in y-z plane may need to increase recursion depth... -sys.setrecursionlimit(10000) -mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts) - -# discretise model -disc = pybamm.Discretisation(mesh, model.default_spatial_methods) -disc.process_model(model) - -# solve model -- simulate one hour discharge -tau = param.process_symbol(pybamm.standard_parameters_lithium_ion.tau_discharge) -t_end = 3600 / tau.evaluate(0) -t_eval = np.linspace(0, t_end, 120) -solution = pybamm.IDAKLUSolver().solve(model, t_eval) - -# TO DO: 2+1D automated plotting -phi_s_cn = pybamm.ProcessedVariable( - model.variables["Negative current collector potential [V]"], - solution.t, - solution.y, - mesh=mesh, -) -phi_s_cp = pybamm.ProcessedVariable( - model.variables["Positive current collector potential [V]"], - solution.t, - solution.y, - mesh=mesh, -) -T = pybamm.ProcessedVariable( - model.variables["X-averaged cell temperature [K]"], - solution.t, - solution.y, - mesh=mesh, -) -l_y = phi_s_cp.y_sol[-1] -l_z = phi_s_cp.z_sol[-1] -y_plot = np.linspace(0, l_y, 21) -z_plot = np.linspace(0, l_z, 21) - - -def plot(t): - fig, ax = plt.subplots(figsize=(15, 8)) - plt.tight_layout() - plt.subplots_adjust(left=-0.1) - - # find t index - ind = (np.abs(solution.t - t)).argmin() - - # negative current collector potential - plt.subplot(131) - phi_s_cn_plot = plt.pcolormesh( - y_plot, - z_plot, - np.transpose(phi_s_cn(y=y_plot, z=z_plot, t=solution.t[ind])), - shading="gouraud", - ) - plt.axis([0, l_y, 0, l_z]) - plt.xlabel(r"$y$") - plt.ylabel(r"$z$") - plt.title(r"$\phi_{s,cn}$ [V]") - plt.set_cmap("cividis") - plt.colorbar(phi_s_cn_plot) - - # positive current collector potential - plt.subplot(132) - phi_s_cp_plot = plt.pcolormesh( - y_plot, - z_plot, - np.transpose(phi_s_cp(y=y_plot, z=z_plot, t=solution.t[ind])), - shading="gouraud", - ) - - plt.axis([0, l_y, 0, l_z]) - plt.xlabel(r"$y$") - plt.ylabel(r"$z$") - plt.title(r"$\phi_{s,cp}$ [V]") - plt.set_cmap("viridis") - plt.colorbar(phi_s_cp_plot) - - # temperature - plt.subplot(133) - T_plot = plt.pcolormesh( - y_plot, - z_plot, - np.transpose(T(y=y_plot, z=z_plot, t=solution.t[ind])), - shading="gouraud", - ) - - plt.axis([0, l_y, 0, l_z]) - plt.xlabel(r"$y$") - plt.ylabel(r"$z$") - plt.title(r"$T$ [K]") - plt.set_cmap("inferno") - plt.colorbar(T_plot) - - plt.subplots_adjust( - top=0.92, bottom=0.15, left=0.10, right=0.9, hspace=0.5, wspace=0.5 - ) - plt.show() - - -plot(solution.t[-1] / 2) diff --git a/results/2plus1D/set_temperature_spm_1plus1D.py b/results/2plus1D/set_temperature_spm_1plus1D.py deleted file mode 100644 index c469d76079..0000000000 --- a/results/2plus1D/set_temperature_spm_1plus1D.py +++ /dev/null @@ -1,51 +0,0 @@ -# -# Example of 1+1D SPM where the temperature can be set by the user -# - -import pybamm -import numpy as np - -# set logging level -pybamm.set_logging_level("INFO") - -model_options = { - "current collector": "potential pair", - "dimensionality": 1, - "thermal": "x-lumped", - "external submodels": ["thermal"], -} -model = pybamm.lithium_ion.SPMe(model_options) - -var = pybamm.standard_spatial_vars -z_pts = 20 -var_pts = {var.x_n: 5, var.x_s: 5, var.x_p: 5, var.r_n: 5, var.r_p: 5, var.z: z_pts} - -sim = pybamm.Simulation(model, var_pts=var_pts, C_rate=2) - -# Set the temperature (in dimensionless form) -# T_av = np.linspace(0, 1, z_pts)[:, np.newaxis] - -z = np.linspace(0, 1, z_pts) -t_eval = np.linspace(0, 0.13, 50) -# step through the solver, setting the temperature at each timestep -for i in np.arange(1, len(t_eval) - 1): - dt = t_eval[i + 1] - t_eval[i] - T_av = (np.sin(2 * np.pi * z) * np.sin(2 * np.pi * 100 * t_eval[i]))[ - :, np.newaxis - ] - external_variables = {"X-averaged cell temperature": T_av} - sim.step(dt, external_variables=external_variables) - -sim.plot( - [ - "Terminal voltage [V]", - "X-averaged total heating [W.m-3]", - "X-averaged cell temperature [K]", - "X-averaged negative particle surface concentration [mol.m-3]", - "X-averaged positive particle surface concentration [mol.m-3]", - "Negative current collector potential [V]", - "Positive current collector potential [V]", - "Local voltage [V]", - ] -) - diff --git a/results/2plus1D/spm_1plus1D.py b/results/2plus1D/spm_1plus1D.py deleted file mode 100644 index 993ab68ae9..0000000000 --- a/results/2plus1D/spm_1plus1D.py +++ /dev/null @@ -1,63 +0,0 @@ -import pybamm -import numpy as np -import sys - -# set logging level -pybamm.set_logging_level("INFO") - -# load (1+1D) SPMe model -options = { - "current collector": "potential pair", - "dimensionality": 1, - "thermal": "lumped", -} -model = pybamm.lithium_ion.SPM(options) - -# create geometry -geometry = model.default_geometry - -# load parameter values and process model and geometry -param = model.default_parameter_values -C_rate = 1 -current_1C = 24 * param.process_symbol(pybamm.geometric_parameters.A_cc).evaluate() -param.update( - { - "Typical current [A]": C_rate * current_1C, - "Initial temperature [K]": 298.15, - "Negative current collector conductivity [S.m-1]": 1e7, - "Positive current collector conductivity [S.m-1]": 1e7, - "Heat transfer coefficient [W.m-2.K-1]": 1, - } -) -param.process_model(model) -param.process_geometry(geometry) - -# set mesh -var = pybamm.standard_spatial_vars -var_pts = {var.x_n: 5, var.x_s: 5, var.x_p: 5, var.r_n: 10, var.r_p: 10, var.z: 15} -# depending on number of points in y-z plane may need to increase recursion depth... -sys.setrecursionlimit(10000) -mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts) - -# discretise model -disc = pybamm.Discretisation(mesh, model.default_spatial_methods) -disc.process_model(model) - -# solve model -- simulate one hour discharge -tau = param.process_symbol(pybamm.standard_parameters_lithium_ion.tau_discharge) -t_end = 3600 / tau.evaluate(0) -t_eval = np.linspace(0, t_end, 120) -solution = model.default_solver.solve(model, t_eval) - -# plot -output_variables = [ - "X-averaged negative particle surface concentration [mol.m-3]", - "X-averaged positive particle surface concentration [mol.m-3]", - # "X-averaged cell temperature [K]", - "Local potenital difference [V]", - "Current collector current density [A.m-2]", - "Terminal voltage [V]", - "Volume-averaged cell temperature [K]", -] -plot = pybamm.QuickPlot(model, mesh, solution, output_variables) -plot.dynamic_plot() diff --git a/results/2plus1D/spm_2plus1D.py b/results/2plus1D/spm_2plus1D.py deleted file mode 100644 index 5f95cd31d0..0000000000 --- a/results/2plus1D/spm_2plus1D.py +++ /dev/null @@ -1,137 +0,0 @@ -import pybamm -import numpy as np -import matplotlib.pyplot as plt -import sys - -# set logging level -pybamm.set_logging_level("DEBUG") - -# load (2+1D) SPM model -options = { - "current collector": "potential pair", - "dimensionality": 2, - "thermal": "x-lumped", -} -model = pybamm.lithium_ion.SPM(options) -model.use_simplify = False # simplifying jacobian slow for large systems - -# create geometry -geometry = model.default_geometry - -# load parameter values and process model and geometry -param = model.default_parameter_values -param.process_model(model) -param.process_geometry(geometry) - -# set mesh -var = pybamm.standard_spatial_vars -var_pts = { - var.x_n: 5, - var.x_s: 5, - var.x_p: 5, - var.r_n: 5, - var.r_p: 5, - var.y: 5, - var.z: 5, -} -# depending on number of points in y-z plane may need to increase recursion depth... -sys.setrecursionlimit(10000) -mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts) - -# discretise model -disc = pybamm.Discretisation(mesh, model.default_spatial_methods) -disc.process_model(model) - -# solve model -- simulate one hour discharge -tau = param.process_symbol(pybamm.standard_parameters_lithium_ion.tau_discharge) -t_end = 3600 / tau.evaluate(0) -t_eval = np.linspace(0, t_end, 120) -solution = pybamm.IDAKLUSolver().solve(model, t_eval) - -# TO DO: 2+1D automated plotting -phi_s_cn = pybamm.ProcessedVariable( - model.variables["Negative current collector potential [V]"], - solution.t, - solution.y, - mesh=mesh, -) -phi_s_cp = pybamm.ProcessedVariable( - model.variables["Positive current collector potential [V]"], - solution.t, - solution.y, - mesh=mesh, -) -T = pybamm.ProcessedVariable( - model.variables["X-averaged cell temperature [K]"], - solution.t, - solution.y, - mesh=mesh, -) -l_y = phi_s_cp.y_sol[-1] -l_z = phi_s_cp.z_sol[-1] -y_plot = np.linspace(0, l_y, 21) -z_plot = np.linspace(0, l_z, 21) - - -def plot(t): - fig, ax = plt.subplots(figsize=(15, 8)) - plt.tight_layout() - plt.subplots_adjust(left=-0.1) - - # find t index - ind = (np.abs(solution.t - t)).argmin() - - # negative current collector potential - plt.subplot(131) - phi_s_cn_plot = plt.pcolormesh( - y_plot, - z_plot, - np.transpose(phi_s_cn(y=y_plot, z=z_plot, t=solution.t[ind])), - shading="gouraud", - ) - plt.axis([0, l_y, 0, l_z]) - plt.xlabel(r"$y$") - plt.ylabel(r"$z$") - plt.title(r"$\phi_{s,cn}$ [V]") - plt.set_cmap("cividis") - plt.colorbar(phi_s_cn_plot) - - # positive current collector potential - plt.subplot(132) - phi_s_cp_plot = plt.pcolormesh( - y_plot, - z_plot, - np.transpose(phi_s_cp(y=y_plot, z=z_plot, t=solution.t[ind])), - shading="gouraud", - ) - - plt.axis([0, l_y, 0, l_z]) - plt.xlabel(r"$y$") - plt.ylabel(r"$z$") - plt.title(r"$\phi_{s,cp}$ [V]") - plt.set_cmap("viridis") - plt.colorbar(phi_s_cp_plot) - - # temperature - plt.subplot(133) - T_plot = plt.pcolormesh( - y_plot, - z_plot, - np.transpose(T(y=y_plot, z=z_plot, t=solution.t[ind])), - shading="gouraud", - ) - - plt.axis([0, l_y, 0, l_z]) - plt.xlabel(r"$y$") - plt.ylabel(r"$z$") - plt.title(r"$T$ [K]") - plt.set_cmap("inferno") - plt.colorbar(T_plot) - - plt.subplots_adjust( - top=0.92, bottom=0.15, left=0.10, right=0.9, hspace=0.5, wspace=0.5 - ) - plt.show() - - -plot(solution.t[-1] / 2) diff --git a/results/2plus1D/spm_2plus1D_tab_grid.py b/results/2plus1D/spm_2plus1D_tab_grid.py deleted file mode 100644 index 19835003c5..0000000000 --- a/results/2plus1D/spm_2plus1D_tab_grid.py +++ /dev/null @@ -1,138 +0,0 @@ -import pybamm -import numpy as np -import matplotlib.pyplot as plt -import sys - -# set logging level -pybamm.set_logging_level("INFO") - -# load (2+1D) SPM model -options = { - "current collector": "potential pair", - "dimensionality": 2, - "thermal": "x-lumped", -} -model = pybamm.lithium_ion.SPM(options) - -# create geometry -geometry = model.default_geometry - -# load parameter values and process model and geometry -param = model.default_parameter_values -param.process_model(model) -param.process_geometry(geometry) - -# set mesh -var = pybamm.standard_spatial_vars -var_pts = { - var.x_n: 5, - var.x_s: 5, - var.x_p: 5, - var.r_n: 5, - var.r_p: 5, - var.y: 5, - var.z: 5, -} -submesh_types = model.default_submesh_types -submesh_types["current collector"] = pybamm.ScikitExponential2DSubMesh -# depnding on number of points in y-z plane may need to increase recursion depth... -sys.setrecursionlimit(10000) -mesh = pybamm.Mesh(geometry, submesh_types, var_pts) - -# discretise model -disc = pybamm.Discretisation(mesh, model.default_spatial_methods) -disc.process_model(model) - -# solve model -- simulate one hour discharge -tau = param.process_symbol(pybamm.standard_parameters_lithium_ion.tau_discharge) -t_end = 3600 / tau.evaluate(0) -t_eval = np.linspace(0, t_end, 120) -solution = model.default_solver.solve(model, t_eval) - -# TO DO: 2+1D automated plotting -phi_s_cn = pybamm.ProcessedVariable( - model.variables["Negative current collector potential [V]"], - solution.t, - solution.y, - mesh=mesh, -) -phi_s_cp = pybamm.ProcessedVariable( - model.variables["Positive current collector potential [V]"], - solution.t, - solution.y, - mesh=mesh, -) -T = pybamm.ProcessedVariable( - model.variables["X-averaged cell temperature [K]"], - solution.t, - solution.y, - mesh=mesh, -) -l_y = phi_s_cp.y_sol[-1] -l_z = phi_s_cp.z_sol[-1] -y_plot = np.linspace(0, l_y, 21) -z_plot = np.linspace(0, l_z, 21) - - -def plot(t): - fig, ax = plt.subplots(figsize=(15, 8)) - plt.tight_layout() - plt.subplots_adjust(left=-0.1) - - # find t index - ind = (np.abs(solution.t - t)).argmin() - - # negative current collector potential - plt.subplot(131) - phi_s_cn_plot = plt.pcolormesh( - y_plot, - z_plot, - np.transpose(phi_s_cn(y=y_plot, z=z_plot, t=solution.t[ind])), - shading="gouraud", - ) - plt.axis([0, l_y, 0, l_z]) - plt.xlabel(r"$y$") - plt.ylabel(r"$z$") - plt.title(r"$\phi_{s,cn}$ [V]") - plt.set_cmap("cividis") - plt.colorbar(phi_s_cn_plot) - - # positive current collector potential - plt.subplot(132) - phi_s_cp_plot = plt.pcolormesh( - y_plot, - z_plot, - np.transpose(phi_s_cp(y=y_plot, z=z_plot, t=solution.t[ind])), - shading="gouraud", - ) - - plt.axis([0, l_y, 0, l_z]) - plt.xlabel(r"$y$") - plt.ylabel(r"$z$") - plt.title(r"$\phi_{s,cp}$ [V]") - plt.set_cmap("viridis") - plt.colorbar(phi_s_cp_plot) - - # temperature - plt.subplot(133) - T_plot = plt.pcolormesh( - y_plot, - z_plot, - np.transpose(T(y=y_plot, z=z_plot, t=solution.t[ind])), - shading="gouraud", - ) - - plt.axis([0, l_y, 0, l_z]) - plt.xlabel(r"$y$") - plt.ylabel(r"$z$") - plt.title(r"$T$ [K]") - plt.set_cmap("inferno") - plt.colorbar(T_plot) - - plt.subplots_adjust( - top=0.92, bottom=0.15, left=0.10, right=0.9, hspace=0.5, wspace=0.5 - ) - plt.show() - - -plot(solution.t[-1] / 2) diff --git a/results/2plus1D/spme_2plus1D.py b/results/2plus1D/spme_2plus1D.py deleted file mode 100644 index c7256b993a..0000000000 --- a/results/2plus1D/spme_2plus1D.py +++ /dev/null @@ -1,137 +0,0 @@ -import pybamm -import numpy as np -import matplotlib.pyplot as plt -import sys - -# set logging level -pybamm.set_logging_level("INFO") - -# load (2+1D) SPMe model -options = { - "current collector": "potential pair", - "dimensionality": 2, - "thermal": "x-lumped", -} -model = pybamm.lithium_ion.SPMe(options) -model.use_simplify = False # simplifying jacobian slow for large systems - -# create geometry -geometry = model.default_geometry - -# load parameter values and process model and geometry -param = model.default_parameter_values -param.process_model(model) -param.process_geometry(geometry) - -# set mesh -var = pybamm.standard_spatial_vars -var_pts = { - var.x_n: 5, - var.x_s: 5, - var.x_p: 5, - var.r_n: 5, - var.r_p: 5, - var.y: 5, - var.z: 5, -} -# depending on number of points in y-z plane may need to increase recursion depth... -sys.setrecursionlimit(10000) -mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts) - -# discretise model -disc = pybamm.Discretisation(mesh, model.default_spatial_methods) -disc.process_model(model) - -# solve model -- simulate one hour discharge -tau = param.process_symbol(pybamm.standard_parameters_lithium_ion.tau_discharge) -t_end = 3600 / tau.evaluate(0) -t_eval = np.linspace(0, t_end, 120) -solution = pybamm.IDAKLUSolver().solve(model, t_eval) - -# TO DO: 2+1D automated plotting -phi_s_cn = pybamm.ProcessedVariable( - model.variables["Negative current collector potential [V]"], - solution.t, - solution.y, - mesh=mesh, -) -phi_s_cp = pybamm.ProcessedVariable( - model.variables["Positive current collector potential [V]"], - solution.t, - solution.y, - mesh=mesh, -) -T = pybamm.ProcessedVariable( - model.variables["X-averaged cell temperature [K]"], - solution.t, - solution.y, - mesh=mesh, -) -l_y = phi_s_cp.y_sol[-1] -l_z = phi_s_cp.z_sol[-1] -y_plot = np.linspace(0, l_y, 21) -z_plot = np.linspace(0, l_z, 21) - - -def plot(t): - fig, ax = plt.subplots(figsize=(15, 8)) - plt.tight_layout() - plt.subplots_adjust(left=-0.1) - - # find t index - ind = (np.abs(solution.t - t)).argmin() - - # negative current collector potential - plt.subplot(131) - phi_s_cn_plot = plt.pcolormesh( - y_plot, - z_plot, - np.transpose(phi_s_cn(y=y_plot, z=z_plot, t=solution.t[ind])), - shading="gouraud", - ) - plt.axis([0, l_y, 0, l_z]) - plt.xlabel(r"$y$") - plt.ylabel(r"$z$") - plt.title(r"$\phi_{s,cn}$ [V]") - plt.set_cmap("cividis") - plt.colorbar(phi_s_cn_plot) - - # positive current collector potential - plt.subplot(132) - phi_s_cp_plot = plt.pcolormesh( - y_plot, - z_plot, - np.transpose(phi_s_cp(y=y_plot, z=z_plot, t=solution.t[ind])), - shading="gouraud", - ) - - plt.axis([0, l_y, 0, l_z]) - plt.xlabel(r"$y$") - plt.ylabel(r"$z$") - plt.title(r"$\phi_{s,cp}$ [V]") - plt.set_cmap("viridis") - plt.colorbar(phi_s_cp_plot) - - # temperature - plt.subplot(133) - T_plot = plt.pcolormesh( - y_plot, - z_plot, - np.transpose(T(y=y_plot, z=z_plot, t=solution.t[ind])), - shading="gouraud", - ) - - plt.axis([0, l_y, 0, l_z]) - plt.xlabel(r"$y$") - plt.ylabel(r"$z$") - plt.title(r"$T$ [K]") - plt.set_cmap("inferno") - plt.colorbar(T_plot) - - plt.subplots_adjust( - top=0.92, bottom=0.15, left=0.10, right=0.9, hspace=0.5, wspace=0.5 - ) - plt.show() - - -plot(solution.t[-1] / 2) diff --git a/results/2plus1D/spmecc.py b/results/2plus1D/spmecc.py deleted file mode 100644 index 0939667226..0000000000 --- a/results/2plus1D/spmecc.py +++ /dev/null @@ -1,59 +0,0 @@ -import pybamm -import numpy as np -import matplotlib.pyplot as plt - -# set logging level -pybamm.set_logging_level("INFO") - -# load current collector and SPMe models -cell_model = pybamm.lithium_ion.SPMe() -cc_model = pybamm.current_collector.EffectiveResistance2D() -models = [cell_model, cc_model] - -# set parameters based on the cell model -param = cell_model.default_parameter_values - -# make current collectors not so conductive, just for illustrative purposes -param.update( - { - "Negative current collector conductivity [S.m-1]": 5.96e6, - "Positive current collector conductivity [S.m-1]": 3.55e6, - } -) - -# process model and geometry, and discretise -for model in models: - param.process_model(model) - geometry = model.default_geometry - param.process_geometry(geometry) - mesh = pybamm.Mesh(geometry, model.default_submesh_types, model.default_var_pts) - disc = pybamm.Discretisation(mesh, model.default_spatial_methods) - disc.process_model(model) - - -# solve current collector model -cc_solution = cc_model.default_solver.solve(cc_model) - -# solve SPMe -- simulate one hour discharge -tau = param.process_symbol(pybamm.standard_parameters_lithium_ion.tau_discharge) -t_end = 3600 / tau.evaluate(0) -t_eval = np.linspace(0, t_end, 120) -solution = cell_model.default_solver.solve(cell_model, t_eval) - -# plot terminal voltage -t, y = solution.t, solution.y -time = pybamm.ProcessedVariable(cell_model.variables["Time [h]"], t, y)(t) -voltage = pybamm.ProcessedVariable(cell_model.variables["Terminal voltage [V]"], t, y) -current = pybamm.ProcessedVariable(cell_model.variables["Current [A]"], t, y)(t) -delta = param.evaluate(pybamm.standard_parameters_lithium_ion.delta) -R_cc = param.process_symbol( - cc_model.variables["Effective current collector resistance [Ohm]"] -).evaluate(t=cc_solution.t, y=cc_solution.y)[0][0] -cc_ohmic_losses = -delta * current * R_cc - -plt.plot(time, voltage(t), label="SPMe") -plt.plot(time, voltage(t) + cc_ohmic_losses, label="SPMeCC") -plt.xlabel("Time [h]") -plt.ylabel("Terminal voltage [V]") -plt.legend() -plt.show() diff --git a/results/2plus1D/user_mesh_spm_1plus1D.py b/results/2plus1D/user_mesh_spm_1plus1D.py deleted file mode 100644 index 9d479d6a66..0000000000 --- a/results/2plus1D/user_mesh_spm_1plus1D.py +++ /dev/null @@ -1,73 +0,0 @@ -import pybamm -import numpy as np -import sys - -# set logging level -pybamm.set_logging_level("INFO") - -# load (1+1D) SPM model -options = { - "current collector": "potential pair", - "dimensionality": 1, - "thermal": "lumped", -} -model = pybamm.lithium_ion.SPM(options) - -# create geometry -geometry = model.default_geometry - -# load parameter values and process model and geometry -param = model.default_parameter_values -C_rate = 1 -current_1C = 24 * param.evaluate(pybamm.geometric_parameters.A_cc) -param.update( - { - "Typical current [A]": C_rate * current_1C, - "Initial temperature [K]": 298.15, - "Negative current collector conductivity [S.m-1]": 1e5, - "Positive current collector conductivity [S.m-1]": 1e5, - "Heat transfer coefficient [W.m-2.K-1]": 1, - "Negative tab centre z-coordinate [m]": 0, # negative tab at bottom - "Positive tab centre z-coordinate [m]": 0.137, # positive tab at top - } -) -param.process_model(model) -param.process_geometry(geometry) - -# set mesh using user-supplied edges in z -z_edges = np.array([0, 0.025, 0.05, 0.1, 0.3, 0.5, 0.7, 0.9, 0.95, 0.975, 1]) -submesh_types = model.default_submesh_types -submesh_types["current collector"] = pybamm.MeshGenerator( - pybamm.UserSupplied1DSubMesh, submesh_params={"edges": z_edges} -) -# Need to make sure var_pts for z is one less than number of edges (variables are -# evaluated at cell centres) -npts_z = len(z_edges) - 1 -var = pybamm.standard_spatial_vars -var_pts = {var.x_n: 5, var.x_s: 5, var.x_p: 5, var.r_n: 10, var.r_p: 10, var.z: npts_z} -# depending on number of points in y-z plane may need to increase recursion depth... -sys.setrecursionlimit(10000) -mesh = pybamm.Mesh(geometry, submesh_types, var_pts) - -# discretise model -disc = pybamm.Discretisation(mesh, model.default_spatial_methods) -disc.process_model(model) - -# solve model -- simulate one hour discharge -tau = param.process_symbol(pybamm.standard_parameters_lithium_ion.tau_discharge) -t_end = 3600 / tau.evaluate(0) -t_eval = np.linspace(0, t_end, 120) -solution = model.default_solver.solve(model, t_eval) - -# plot -output_variables = [ - "X-averaged negative particle surface concentration [mol.m-3]", - "X-averaged positive particle surface concentration [mol.m-3]", - # "X-averaged cell temperature [K]", - "Local current collector potential difference [V]", - "Current collector current density [A.m-2]", - "Terminal voltage [V]", - "Volume-averaged cell temperature [K]", -] -plot = pybamm.QuickPlot(model, mesh, solution, output_variables) -plot.dynamic_plot() diff --git a/results/change_settings/change_solver_tolerances.py b/results/change_settings/change_solver_tolerances.py deleted file mode 100644 index 81261e796d..0000000000 --- a/results/change_settings/change_solver_tolerances.py +++ /dev/null @@ -1,52 +0,0 @@ -# -# Compare solution of li-ion battery models when changing solver tolerances -# -import numpy as np -import pybamm - -pybamm.set_logging_level("INFO") - -# load model -model = pybamm.lithium_ion.DFN() - - -# process and discretise -param = model.default_parameter_values -param.process_model(model) -geometry = model.default_geometry -param.process_geometry(geometry) -mesh = pybamm.Mesh(geometry, model.default_submesh_types, model.default_var_pts) -disc = pybamm.Discretisation(mesh, model.default_spatial_methods) -disc.process_model(model) - -# tolerances (rtol, atol) -tols = [[1e-8, 1e-8], [1e-6, 1e-6], [1e-3, 1e-6], [1e-3, 1e-3]] - -# solve model -solutions = [None] * len(tols) -voltages = [None] * len(tols) -voltage_rmse = [None] * len(tols) -labels = [None] * len(tols) -t_eval = np.linspace(0, 0.17, 100) -for i, tol in enumerate(tols): - solver = pybamm.ScikitsDaeSolver(rtol=tol[0], atol=tol[1]) - solutions[i] = solver.solve(model, t_eval) - voltages[i] = pybamm.ProcessedVariable( - model.variables["Terminal voltage [V]"], - solutions[i].t, - solutions[i].y, - mesh=mesh, - )(solutions[i].t) - voltage_rmse[i] = pybamm.rmse(voltages[0], voltages[i]) - labels[i] = "rtol = {}, atol = {}".format(tol[0], tol[1]) - -# print RMSE voltage errors vs tighest tolerance -for i, tol in enumerate(tols): - print( - "rtol = {}, atol = {}, solve time = {} s, Voltage RMSE = {}".format( - tol[0], tol[1], solutions[i].solve_time, voltage_rmse[i] - ) - ) -# plot -plot = pybamm.QuickPlot([model] * len(solutions), mesh, solutions, labels=labels) -plot.dynamic_plot() diff --git a/results/change_settings/compare_var_pts.py b/results/change_settings/compare_var_pts.py deleted file mode 100644 index 2153a96f62..0000000000 --- a/results/change_settings/compare_var_pts.py +++ /dev/null @@ -1,68 +0,0 @@ -# -# Compare solution of li-ion battery models when varying the number of grid points -# -import numpy as np -import pybamm -import matplotlib.pyplot as plt - -pybamm.set_logging_level("INFO") - -# choose number of points per domain (all domains will have same npts) -Npts = [30, 20, 10, 5] - -# create models -models = [None] * len(Npts) -for i, npts in enumerate(Npts): - models[i] = pybamm.lithium_ion.DFN(name="Npts = {}".format(npts)) - -# load parameter values and process models and geometry -param = models[0].default_parameter_values -for model in models: - param.process_model(model) - -# set mesh -meshes = [None] * len(models) - -# create geometry and discretise models -var = pybamm.standard_spatial_vars -for i, model in enumerate(models): - geometry = model.default_geometry - param.process_geometry(geometry) - var_pts = { - var.x_n: Npts[i], - var.x_s: Npts[i], - var.x_p: Npts[i], - var.r_n: Npts[i], - var.r_p: Npts[i], - } - meshes[i] = pybamm.Mesh(geometry, models[-1].default_submesh_types, var_pts) - disc = pybamm.Discretisation(meshes[i], model.default_spatial_methods) - disc.process_model(model) - -# solve model and plot voltage -solutions = [None] * len(models) -voltages = [None] * len(models) -voltage_rmse = [None] * len(models) -t_eval = np.linspace(0, 0.17, 100) -for i, model in enumerate(models): - solutions[i] = model.default_solver.solve(model, t_eval) - voltages[i] = pybamm.ProcessedVariable( - model.variables["Terminal voltage [V]"], - solutions[i].t, - solutions[i].y, - mesh=meshes[i], - )(solutions[i].t) - voltage_rmse[i] = pybamm.rmse(voltages[0], voltages[i]) - plt.plot(solutions[i].t, voltages[i], label=model.name) - -for i, npts in enumerate(Npts): - print( - "npts = {}, solve time = {} s, Voltage RMSE = {}".format( - npts, solutions[i].solve_time, voltage_rmse[i] - ) - ) - -plt.xlabel(r"$t$") -plt.ylabel("Voltage [V]") -plt.legend() -plt.show() diff --git a/results/drive_cycles/US06_simulation.py b/results/drive_cycles/US06_simulation.py deleted file mode 100644 index 64f47fd5e5..0000000000 --- a/results/drive_cycles/US06_simulation.py +++ /dev/null @@ -1,42 +0,0 @@ -# -# Simulate drive cycle loaded from csv file -# -import pybamm -import numpy as np - -# load model -pybamm.set_logging_level("INFO") -model = pybamm.lithium_ion.SPMe() - -# create geometry -geometry = model.default_geometry - -# load parameter values and process model and geometry -param = model.default_parameter_values -param["Current function"] = "[current data]US06" -param.process_model(model) -param.process_geometry(geometry) - -# set mesh -var = pybamm.standard_spatial_vars -var_pts = {var.x_n: 10, var.x_s: 10, var.x_p: 10, var.r_n: 5, var.r_p: 5} -mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts) - -# discretise model -disc = pybamm.Discretisation(mesh, model.default_spatial_methods) -disc.process_model(model) - -# simulate US06 drive cycle -tau = param.evaluate(pybamm.standard_parameters_lithium_ion.tau_discharge) -t_eval = np.linspace(0, 600 / tau, 600) - -# need to increase max solver steps if solving DAEs along with an erratic drive cycle -solver = pybamm.CasadiSolver() -if isinstance(solver, pybamm.DaeSolver): - solver.max_steps = 10000 - -solution = solver.solve(model, t_eval) - -# plot -plot = pybamm.QuickPlot(model, mesh, solution) -plot.dynamic_plot() diff --git a/results/drive_cycles/car_current_simulation.py b/results/drive_cycles/car_current_simulation.py deleted file mode 100644 index 84d37e4f17..0000000000 --- a/results/drive_cycles/car_current_simulation.py +++ /dev/null @@ -1,63 +0,0 @@ -# -# Simulate user-defined current profile -# -import pybamm -import numpy as np - - -def car_current(t): - """ - Piecewise constant current as a function of time in seconds. This is adapted - from the file getCarCurrent.m, which is part of the LIONSIMBA toolbox [1]_. - - References - ---------- - .. [1] M Torchio, L Magni, R Bushan Gopaluni, RD Braatz, and D. Raimondoa. - LIONSIMBA: A Matlab framework based on a finite volume model suitable - for Li-ion battery design, simulation, and control. Journal of The - Electrochemical Society, 163(7):1192-1205, 2016. - """ - - current = ( - 1 * (t >= 0) * (t <= 50) - - 0.5 * (t > 50) * (t <= 60) - + 0.5 * (t > 60) * (t <= 210) - + 1 * (t > 210) * (t <= 410) - + 2 * (t > 410) * (t <= 415) - + 1.25 * (t > 415) * (t <= 615) - - 0.5 * (t > 615) - ) - - return current - - -# load model -pybamm.set_logging_level("INFO") -model = pybamm.lithium_ion.SPMe() - -# create geometry -geometry = model.default_geometry - -# load parameter values and process model and geometry -param = model.default_parameter_values -param["Current function"] = car_current -param.process_model(model) -param.process_geometry(geometry) - -# set mesh -mesh = pybamm.Mesh(geometry, model.default_submesh_types, model.default_var_pts) - -# discretise model -disc = pybamm.Discretisation(mesh, model.default_spatial_methods) -disc.process_model(model) - -# simulate car current for 30 minutes -tau = param.process_symbol( - pybamm.standard_parameters_lithium_ion.tau_discharge -).evaluate(0) -t_eval = np.linspace(0, 1800 / tau, 600) -solution = model.default_solver.solve(model, t_eval) - -# plot -plot = pybamm.QuickPlot(model, mesh, solution) -plot.dynamic_plot() diff --git a/results/drive_cycles/discharge_rest.py b/results/drive_cycles/discharge_rest.py deleted file mode 100644 index 9089648767..0000000000 --- a/results/drive_cycles/discharge_rest.py +++ /dev/null @@ -1,76 +0,0 @@ -# -# Simulate discharge followed by rest -# -import pybamm -import numpy as np -import matplotlib.pyplot as plt - -# load model -pybamm.set_logging_level("INFO") -model = pybamm.lithium_ion.SPMe() - -# create geometry -geometry = model.default_geometry - -# load parameter values and process model and geometry -param = model.default_parameter_values -param.process_model(model) -param.process_geometry(geometry) - -# set mesh -mesh = pybamm.Mesh(geometry, model.default_submesh_types, model.default_var_pts) - -# discretise model -disc = pybamm.Discretisation(mesh, model.default_spatial_methods) -disc.process_model(model) - -# solve model during discharge stage (1 hour) -tau = param.process_symbol(pybamm.standard_parameters_lithium_ion.tau_discharge) -t_end = 3600 / tau.evaluate(0) -t_eval1 = np.linspace(0, t_end, 120) -solution1 = model.default_solver.solve(model, t_eval1) - -# process variables for later plotting -time1 = pybamm.ProcessedVariable(model.variables["Time [h]"], solution1.t, solution1.y) -voltage1 = pybamm.ProcessedVariable( - model.variables["Terminal voltage [V]"], solution1.t, solution1.y, mesh=mesh -) -current1 = pybamm.ProcessedVariable( - model.variables["Current [A]"], solution1.t, solution1.y, mesh=mesh -) - -# solve again with zero current, using last step of solution1 as initial conditions -# update the current to be zero -param["Current function"] = "[zero]" -param.update_model(model, disc) -# Note: need to update model.concatenated_initial_conditions *after* update_model, -# as update_model updates model.concatenated_initial_conditions, by concatenting -# the (unmodified) initial conditions for each variable -model.concatenated_initial_conditions = solution1.y[:, -1][:, np.newaxis] - -# simulate 1 hour of rest -t_start = solution1.t[-1] -t_end = t_start + 3600 / tau.evaluate(0) -t_eval2 = np.linspace(t_start, t_end, 120) -solution2 = model.default_solver.solve(model, t_eval2) - -# process variables for later plotting -time2 = pybamm.ProcessedVariable(model.variables["Time [h]"], solution2.t, solution2.y) -voltage2 = pybamm.ProcessedVariable( - model.variables["Terminal voltage [V]"], solution2.t, solution2.y, mesh=mesh -) -current2 = pybamm.ProcessedVariable( - model.variables["Current [A]"], solution2.t, solution2.y, mesh=mesh -) - -# plot -plt.subplot(121) -plt.plot(time1(t_eval1), voltage1(t_eval1), time2(t_eval2), voltage2(t_eval2)) -plt.xlabel("Time [h]") -plt.ylabel("Voltage [V]") -plt.subplot(122) -z = np.linspace(0, 1, 10) -plt.plot(time1(t_eval1), current1(t_eval1), time2(t_eval2), current2(t_eval2)) -plt.xlabel("Time [h]") -plt.ylabel("Current [A]") -plt.show() diff --git a/results/drive_cycles/user_sin_current_simulation.py b/results/drive_cycles/user_sin_current_simulation.py deleted file mode 100644 index c3f56d39b9..0000000000 --- a/results/drive_cycles/user_sin_current_simulation.py +++ /dev/null @@ -1,63 +0,0 @@ -# -# Simulate user-defined current profile which takes parameters -# -import pybamm -import numpy as np - - -# create user-defined function -def my_fun(t, A, omega): - return A * np.sin(2 * np.pi * omega * t) - - -# choose amplitude and frequencies -A = pybamm.electrical_parameters.I_typ -frequencies = [0.1, 1] - -# load models (need to create new instances of model, not copies) -pybamm.set_logging_level("INFO") -models = [None] * len(frequencies) -for i in range(len(frequencies)): - models[i] = pybamm.lithium_ion.SPM() - -# load parameter values and process models -param = models[0].default_parameter_values -for i, frequency in enumerate(frequencies): - - def current(t): - return my_fun(t, A, frequency) - - param.update({"Current function": current}) - param.process_model(models[i]) - -# discretise models -for model in models: - # create geometry - geometry = model.default_geometry - param.process_geometry(geometry) - mesh = pybamm.Mesh( - geometry, models[-1].default_submesh_types, model.default_var_pts - ) - disc = pybamm.Discretisation(mesh, model.default_spatial_methods) - disc.process_model(model) - -# Example: simulate for 30 seconds -simulation_time = 30 # end time in seconds -tau = param.process_symbol( - pybamm.standard_parameters_lithium_ion.tau_discharge -).evaluate(0) - -# loop over frequencies -solutions = [None] * len(frequencies) -labels = [None] * len(frequencies) -for i, frequency in enumerate(frequencies): - # need enough timesteps to resolve output - npts = 20 * simulation_time * frequency - t_eval = np.linspace(0, simulation_time / tau, npts) - solutions[i] = model.default_solver.solve(model, t_eval) - labels[i] = "Frequency: {} Hz".format(frequency) - -# plot -output_variables = ["Current [A]", "Terminal voltage [V]"] -plot = pybamm.QuickPlot(models, mesh, solutions, output_variables, labels) -plot.dynamic_plot() diff --git a/tests/unit/test_parameters/test_standard_parameters_lead_acid.py b/tests/unit/test_parameters/test_standard_parameters_lead_acid.py index 0521776d59..fed248f02f 100644 --- a/tests/unit/test_parameters/test_standard_parameters_lead_acid.py +++ b/tests/unit/test_parameters/test_standard_parameters_lead_acid.py @@ -17,7 +17,7 @@ def test_scipy_constants(self): def test_all_defined(self): parameters = pybamm.standard_parameters_lead_acid parameter_values = pybamm.lead_acid.BaseModel().default_parameter_values - output_file = "results/2019_08_sulzer_thesis/parameters.txt" + output_file = "lead_acid_parameters.txt" pybamm.print_parameters(parameters, parameter_values, output_file) # test print_parameters with dict and without C-rate del parameter_values["Cell capacity [A.h]"]