%% Protonic Membrane model
%% Load and parse input from given json files
% The source of the json files can be seen in :battmofile:`protonicMembrane.json<ProtonicMembrane/jsonfiles/protonicMembrane.json>` and
% :battmofile:`1d-PM-geometry.json<ProtonicMembrane/jsonfiles/1d-PM-geometry.json>`
filename = fullfile(battmoDir(), 'ProtonicMembrane', 'jsonfiles', 'protonicMembrane.json');
jsonstruct_material = parseBattmoJson(filename);
filename = fullfile(battmoDir(), 'ProtonicMembrane', 'jsonfiles', '1d-PM-geometry.json');
jsonstruct_geometry = parseBattmoJson(filename);
jsonstruct = mergeJsonStructs({jsonstruct_material, jsonstruct_geometry});
%% Input structure setup
% We setup the input parameter structure which will we be used to instantiate the model
inputparams = ProtonicMembraneInputParams(jsonstruct);
%%
% We setup the grid, which is done by calling the function :battmo:`setupProtonicMembraneGrid`
[inputparams, gen] = setupProtonicMembraneGrid(inputparams, jsonstruct);
%% Model setup
% We instantiate the model for the protonic membrane cell
model = ProtonicMembrane(inputparams);
%%
% The model is equipped for simulation using the following command (this step may become unnecessary in future versions)
model = model.setupForSimulation();
%% Initial state setup
% We setup the initial state using a default setup included in the model
state0 = model.setupInitialState();
%% Schedule
% We setup the schedule, which means the timesteps and also the control we want to use. In this case we use current
% control and the current equal to zero (see :battmofile:`here <ProtonicMembrane/jsonfiles/protonicMembrane.json#118>`).
%
% We compute the steady-state solution and the time stepping here does not correspond to time values but should be seen
% as step-wise increase of the effect of the non-linearity (in particular in the expression of the conductivity which
% includes highly nonlineaer effect with the exponential terms. We do not detail here the method).
schedule = model.Control.setupSchedule(jsonstruct);
%% Simulation
% We run the simulation
[~, states, report] = simulateScheduleAD(state0, model, schedule);
%% Plotting
%
%%
% We setup som shortcuts for convenience and introduce plotting options
an = 'Anode';
ct = 'Cathode';
elyte = 'Electrolyte';
ctrl = 'Control';
set(0, 'defaultlinelinewidth', 3);
set(0, 'defaultaxesfontsize', 15);
%%
% We recover the position of the mesh cell of the discretization grid. This is used when plotting the spatial
% distribution of some of the variables.
xc = model.(elyte).grid.cells.centroids(:, 1);
%%
% We consider the solution obtained at the last time step, which corresponds to the solution at steady-state. The second
% line adds to the state variable all the variables that are derived from our primary unknowns.
state = states{end};
state = model.addVariables(state, schedule.control);
%%
% Plot of electromotive potential
figure
plot(xc, state.(elyte).pi)
title('Electromotive potential (\pi)')
xlabel('x / m')
ylabel('\pi / V')
%%
% Plot of electronic chemical potential
figure
plot(xc, state.(elyte).pi - state.(elyte).phi)
title('Electronic chemical potential (E)')
xlabel('x / m')
ylabel('E / V')
%%
% Plot of electrostatic potential
figure
plot(xc, state.(elyte).phi)
title('Electrostatic potential (\phi)')
xlabel('x / m')
ylabel('\phi / V')
%%
% Plot of the conductivity
figure
plot(xc, log(state.(elyte).sigmaEl))
title('logarithm of conductivity (\sigma)')
xlabel('x / m')
xlabel('log(\sigma/Siemens)')
%% Evolution of the Faradic efficiency
% We increase the current density from 0 to 1 A/cm^2 and plot the faraday efficiency.
%
% We sample the current value from 0 to 1 A/cm^2.
%
Is = linspace(0, 1*ampere/((centi*meter)^2), 20);
%%
% We run the simulation for each current value and collect the results in the :code:`endstates`.
%
endstates = {};
for iI = 1 : numel(Is)
model.Control.I = Is(iI);
[~, states, report] = simulateScheduleAD(state0, model, schedule);
state = states{end};
state = model.addVariables(state, schedule.control);
endstates{iI} = state;
end
%%
% We plot the profile of the electromotive potential for the mininum and maximum current values.
%
figure
hold on
unit = ampere/((centi*meter)^2); % shortcut
state = endstates{1};
plot(xc, state.(elyte).pi, 'displayname', sprintf('I=%g A/cm^2', Is(1)/unit));
state = endstates{end};
plot(xc, state.(elyte).pi, 'displayname', sprintf('I=%g A/cm^2', Is(end)/unit));
title('Electromotive potential (\pi)')
xlabel('x / m')
ylabel('\pi / V')
legend
%%
% We retrieve and plot the cell potential
%
E = cellfun(@(state) state.(an).pi - state.(ct).pi, endstates);
figure
plot(Is/unit, E, '*-');
xlabel('Current density / A/cm^2')
ylabel('E / V')
title('Cell potential');
%%
% We retrieve and plot the Faradic efficiency
%
feff = cellfun(@(state) state.(an).iHp/state.(an).i, endstates);
figure
plot(Is/unit, feff, '*-');
xlabel('Current density / A/cm^2')
ylabel('Faradic efficiency / -')
title('Faradic efficiency');