solvers

Solve the diffusion equation.

solvers.solver(z_init, geometry, cvar, MU_fun, L_fun, yVa_fun, nt, dt, ideal_lattice, k_dislo, k_pores, rho_dislo, rho_pores, DVa_fun, Va_method, saved_steps, BC, show_completion, verbose, stencil, logger)[source]

Solve diffusion equation.

Explicit (forward Euler) solver with non-equilibrium vacancies and pores, finite sink strengths that represent dislocation climb and pore growth.

Parameters:
  • z_init (1D array) – Initial node positions, shape (nz,).

  • geometry (str) – Domain geometry (planar, cylindrical or spherical).

  • cvar (composition_variables.CompositionVariables) – Composition variables (x, y, c, Vm, …).

  • MU_fun (function) – Compute chemical potential from site fractions (see MU_funy in simu.System).

  • L_fun (func) – Same for phenomenological coefficients.

  • yVa_fun (function) – Compute equilibrium vacancy site fraction, see simu.System.add_Va_model().

  • nt (int) – Number of time steps, including time 0.

  • dt (float) – Time step (s).

  • ideal_lattice (bool) – Whether lattice is ideal, i.e., vacancy fraction is kept at equilibrium.

  • k_dislo (float or str) – Sink strength associated with dislocation climb (str should be a valid file name in current job folder)

  • k_pores (float or str) – Sink strength associated with pore growth (str should be a valid file name in current job folder)

  • rho_dislo (float or str) – Dislocation density, used to compute k_dislo from local DVa.

  • rho_pores (float) – Density used to compute k_pores from local DVa.

  • DVa_fun (function) – Compute vacancy diffusion coefficient from atom fractions and site fractions, see thermo_functions.make_DVa_fun().

  • Va_method (str) – Type of function for equilibrium vacancy site fraction. Determined based on thermo_method in simu.System.add_thermo_functions().

  • saved_steps (list) – Steps for which results are returned.

  • BC (dict) – Types and functions of left and right BC, see simu.Simulation.add_BC().

  • show_completion (bool) – Print completion rate while solver is running.

  • verbose (str) – Verbosity level.

  • stencil (str) – Name of discretization stencil. See compute_resistance().

Returns:

  • res (dict of dicts) – Dict of stored variables:

    {n: {var: array for var in variables} for n in saved_steps}

    where variables are:

    • z : node positions, 1D array

    • c : concentrations, 2D array

    • fp : pore fractions, 1D array

    • mu : chemical potentials, 2D array

    • Jlat : fluxes in lattice-fixed frame, 2D array

    • yVa_eq : equilibrium vacancy fraction, 1D array

    • v : velocity field of lattice in laboratory frame, 1D array.

    • gamma_V : relative volume variation rate associated with molar volume variation, 1D array.

    • gamma_p : relative volume variation rate associated with pore growth, 1D array.

    • alpha_d : sink term associated with dislocation climb, 1D array.

    • alpha_p : sink term associated with pore growth, 1D array.

    • L : Onsager coefficients, 3D array

    • deformation : cumulative deformation in diffusion direction (eps_zz), 1D array.

    The first 3 variables (z, c, fp) are stored at the start of the time step, while the others are stored at the end of the step.

  • dur (float) – Solver run time in s.

solvers.solver_core(z, cnod, c, y_Va, dt, Jlat, Vk, V0, Vp, V, Vm, alpha_d, alpha_p, geometry)[source]

Compute fluxes in laboratory frame and then new concentrations.

See argument definitions in solver().

solvers.compute_alpha_ideal(z, cnod, c, y_Va, dt, Jlat, Vk, V0, Vp, V, Vm, yVa_fun, Va_method, geometry)[source]

Compute ideal sink terms.

See argument definitions in solver().

The term related to pore growth, alpha_p, is zero.

The term related to dislocation climb, alpha_d, is computed from analytical expression if the equilibrium vacancy fraction is composition-fixed. If not, an iterative method is used instead, with the analytical version as a starting point. See Gheno 2022 [1] for details.

The error increases when partial molar volumes are different. The error is reduced by increasing the number of loops. Two loops is found to be sufficient in practice (y0 - y0_eq = 1e-13 - 1e-14).

solvers.compute_boundary_fluxes(t, MU_diff, dz, BC, MU_fun, L, L_fun, Vk)[source]

Compute fluxes on domain boundaries.

Parameters:
  • t (float) – Time (s).

  • MU_diff (2D array) – Diffusion potentials (MU_k - MU_0), initial shape (ninds + 1, nz - 1).

  • dz (1D array) – Space step (m).

  • BC (dict) – Types and functions of left and right BC, see simu.Simulation.add_BC().

  • MU_fun (function) – Compute chemical potentials from site fractions.

  • L (3D array) – Onsager coefficients, initial shape (ninds + 1, ninds + 1, nz - 1).

  • Vk (1D array) – Partial molar volumes, shape (ninds + 2).

Returns:

  • J[‘left’] (1D array) – Flux on left-hand boundary, shape (ninds + 1,).

  • J[‘right’] (1D array) – Flux on right-hand boundary, shape (ninds + 1,).

solvers.compute_resistance(comps, L, dz, stencil)[source]

Compute diffusion resistance.

Defined as dz/L, with L the Onsager coefficients.

Parameters:
  • comps (list of str) – System constituents.

  • L (3D array) – Onsager coefficients, initial shape (ninds + 1, ninds + 1, nz - 1).

  • dz (1D array) – Space step, initial shape (nz - 1,).

  • stencil (str) – Name of discretization stencil. Possible values: ‘H’, ‘A’, ‘G’.

Raises:

Exception – If stencil is not supported.

Returns:

RL – Diffusion resistance, initial shape (ninds + 1, nz - 2).

Return type:

2D array

solvers.remesh(n, z, dz_min, dz_max, cmid, deformation, verbose, logger)[source]

Add or delete nodes.

Parameters:
  • n (int) – Time step.

  • z (1D array) – Node positions, initial shape (nz,).

  • dz_min (float) – Minimum of initial dz.

  • dz_max (float) – Maximum of initial dz.

  • cmid (2D array) – Concentrations, initial shape (ninds + 1, nz - 1).

  • deformation (1D array) – Cumulated deformation, initial shape (nz - 1,).

  • verbose (str) – Verbosity level.

Returns:

  • new_z (1D array) – New node positions.

  • new_cmid (2D array) – New concentrations.

  • new_deformation (1D array) – New deformation.