List of all functions

advect_interface_regrid!(Ri, V_ip, dt, x_left, x_right, C_left, C_right, nr)

Update the interface position and calculate new grids based on the advection velocity. Units may differ from SI units if non-dimensionalisation has been performed.

Arguments

• Ri::Float64: Radii [interface total length in [m].
• V_ip::Float64: Advection velocity in [m/s].
• dt::Float64: Time step in [s].
• x_left::Vector{Float64}: Left distance nodes in [m].
• x_right::Vector{Float64}: Right distance nodes in [m].
• C_left::Vector{Float64}: Composition values on the left nodes in [-].
• C_right::Vector{Float64}: Composition values on the right nodes in [-].
• nr::Vector{Int}: Resolution of the left and the right grid.

Returns

• Fl_regrid::Int: Flag indicating if regridding was performed (1) or not (0).
• x_left::Vector{Float64}: Updated left distance nodes.
• x_right::Vector{Float64}: Updated right distance nodes.
• C_left::Vector{Float64}: Updated composition values on the left nodes.
• C_right::Vector{Float64}: Updated composition values on the right nodes.
• nr::Vector{Int}: Updated resolution.
• Ri::Float64: Updated radii [interface total length].

blocktest(L1, R1, L2, R2)

Constructs a block matrix and a block vector from given input matrices and vectors.

Arguments

• L1::Matrix: The first input matrix.
• R1::Vector: The first input vector.
• L2::Matrix: The second input matrix.
• R2::Vector: The second input vector.

Returns

• Lblock::SparseMatrixCSC: The block matrix constructed from L1 and L2.
• Rblock::Vector: The block vector constructed from R1 and R2.

calc_mass_err(Mass, Mass0)

Calculate the mass error between the final mass Mass[end] and the initial mass Mass0.

Arguments

• Mass::Vector: A vector containing the mass values in [-].
• Mass0::Number: The initial mass value in [-].

Output

• ErrM::Number: The calculated mass error.

calc_mass_vol(x_left, x_right, C_left, C_right, n, rho)

Calculate the total mass based on the volume of the phase.

Arguments

• x_left::Float64: Left nodes in [m].
• x_right::Float64: Right nodes in [m].
• C_left::Vector{Float64}: Composition values of the left phase in [-].
• C_right::Vector{Float64}: Composition values of the right phase in [-].
• n::Int: Number which defines the geometry.
• rho::Vector{Float64}: Densities of the left and right phase [kg/m³].

Returns

• Mtot::Float64: The total mass.

calc_volume(x1, x2, ndim)

Calculation of all volumes. The density in both phases is constant. Subsequently, shrinking and expanding volumes are not considered.

Arguments

• x1: Nodes of the left phase.
• x2: Nodes of the right phase.
• ndim: Geometry factor.

Returns

• V1: Array of volumes for the left phase.
• V2: Array of volumes for the right phase.
• dVC: Array of total volume changes.

calculate_density(X_A, Y_A, rho_left, rho_right, C_leftB, C_rightB, T)

Calculate the density of the phases at a given temperature T using interpolation in 2D.

Arguments

• X_A: X-axis values for interpolation (composition X in [-])
• Y_A: Y-axis values for interpolation (temperature T in [K])
• rho_left: Left density values for interpolation in [kg/m^3]
• rho_right: Right density values for interpolation in [kg/m^3]
• C_leftB: Left composition values for interpolation in [-]
• C_rightB: Right composition values for interpolation in [-]
• T: Temperature at which to calculate the density in [K]

Returns

• rho: Array of normalized densities

coeff_trans_line(eq_values)

Extracts coefficients for linear least squares from the input eq_values.

Arguments

• eq_values: A 2x3 matrix containing the coefficients for the upper and lower transition lines.

Returns

• coeff_up: A 1x3 vector containing the coefficients for the upper transition line.
• coeff_do: A 1x3 vector containing the coefficients for the lower transition line.

composition(coeff_up, coeff_do, T)

Compute the composition of two components A and B at a given temperature.

Arguments

• coeff_up::Vector{Float64}: Coefficients for the composition of component B.
• coeff_do::Vector{Float64}: Coefficients for the composition of component A.
• T::Float64: Temperature at which to compute the composition in [K].

Returns

• C_left::Float64: Composition of component A in [-].
• C_right::Float64: Composition of component B in [-].

construct_matrix_fem(x_left, x_right, C_left, C_right, D_l, D_r, dt, n, res)

Constructs the global matrix for the FEM solver in a diffusion couple problem. Units may differ from SI units if non-dimensionalisation has been performed.

Arguments

• x_left::Vector{Float64}: Left grid spatial points in [m].
• x_right::Vector{Float64}: Right grid spatial points in [m].
• C_left::Vector{Float64}: Composition values of the left phase in [-].
• C_right::Vector{Float64}: Composition values of the right phase in [-].
• D_l::Float64: Diffusion coefficient on the left side in [m²/s].
• D_r::Float64: Diffusion coefficient on the right side in [m²/s].
• dt::Float64: Time step in [s].
• n::Int: Geometry definition.
• res::Vector{Float64}: Resolution.

Returns

• L_g::SparseMatrixCSC{Float64, Int}: Global stiffness matrix.
• R_g::Vector{Float64}: Global RHS vector.
• Co_l::Vector{Float64}: Stores left side composition values before the update.
• Co_r::Vector{Float64}: Stores right side composition values before the update.

create_grid!(Ri, nr, MRefin,verbose)

Create grid with or without variable spacing.Units may differ from SI units if non-dimensionalisation has been performed.

Arguments

• Ri::Vector{Float64}: Initial radii [interface total length] in [m].
• nr::Vector{Int}: Resolution vector.
• MRefin::Int: The refinement factor for the grid.
• verbose::Bool: A boolean indicating whether to print additional information.

Returns

• x_left::Matrix{Float64}: Left nodes in [m].
• x_right::Matrix{Float64}: Right nodes in [m].
• dx1::Float64: The grid spacing on the left side in [m].
• dx2::Float64: The grid spacing on the right side in [m].
• x0::Matrix{Float64}: Initial grid spacing for the whole domain in [m].

define_new_grid(Ri, nr, Rfact, verbose)

This function defines a new grid based on the given parameters. Units may differ from SI units if non-dimensionalisation has been performed.

Arguments

• Ri: Radii [interface total length] in [m].
• nr: Resolution of nodes on the left and right side.
• Rfact: Grid refinement factor.
• verbose: A boolean indicating whether to print additional information.

Returns

• Ri: Radii [interface total length] in [m].
• nr: Resolution of nodes on the left and right sides.
• x_left:Left nodes in [m].
• x_right: Right nodes in [m].
• dx_left: Left grid spacing in [m].
• dx_right: Right grid spacing in [m].
• Sc_left: The scaling factor for the left side.
• Sc_right: The scaling factor for the right side.

fill_matrix!(C, x, D, dt, ndim, nels)

fillmatrix! function fills the global matrices Lg and R_g with the corresponding local matrices and vectors. Units may differ from SI units if non-dimensionalisation has been performed.

Arguments

• C: Composition matrix in [-].
• x: Spatial grid points.
• D: Diffusion coefficient in [m²/s].
• dt: Time step in [s].
• ndim: Geometry factor.
• nels: Number of elements.

Returns

• L_g: Global LHS matrix.
• R_g: Global RHS vector.

find_dt(dx1, dx2, V_ip, D_l, D_r, CFL)

Find the important time step dt based on the given parameters. The function calculates the time step dt based on the advection and diffusion properties of the system. Usually, advection time scale dtV are more dominat than diffusion time scale dtD. However, we included a dumping of dt, if dt>dtDto ensure the visibility of diffusion processes within the code. If the advection velocityV_ipis zero,dtD` is used instead. Units may differ from SI units if non-dimensionalisation has been performed.

Arguments

• dx1: Left spatial step size next to the interface.
• dx2: Right spatial step size next to the interface.
• V_ip: The advection velocity in [m/s].
• D_l: The diffusion coefficient on the left side in [m²/s].
• D_r: The diffusion coefficient on the right side in [m²/s].
• CFL: The Courant-Friedrichs-Lewy number.

Returns

• dt: The calculated time step.

interp2(X1d, Y1d, Z2d, xi, yi)

Perform 2D bilinear interpolation within bounds.

Arguments

• X1d::Vector: 1D monotonic (increasing) array representing the x-coordinates.
• Y1d::Vector: 1D monotonic (increasing) array representing the y-coordinates.
• Z2d::Matrix: 2D array representing the values to be interpolated.
• xi::Vector: 1D array of x-values to interpolate.
• yi::Vector: 1D array of y-values to interpolate.

Returns

• zi::Vector: 1D array of interpolated values corresponding to xi and yi.

Notes

• The lengths of xi and yi must be the same.
• The values in xi and yi should be within the bounds of X1d and Y1d respectively.
• The function uses 2D bilinear interpolation to compute the interpolated values.

linear_interpolation_1D(x, y, x_interp)

Perform linear interpolation in 1D.

Arguments

• x::AbstractArray: The x coordinates of data points.
• y::AbstractArray: The y coordinates of data points.
• x_interp::Real: The x coordinate at which interpolation is desired.

Returns

• y_interp::Real: The interpolated y value at x_interp.

linspace_interface(L1, L2, LIP, nx1, nx2, dX1_dXN)

This function calculates the adaptive grid depending on the position of the interface. Note: This function assumes that L1 < L2 and nx1 <= nx2. Units may differ from SI units if non-dimensionalisation has been performed.

Arguments

• L1: Length of the left side in [m].
• L2: Length of the right side in [m].
• LIP: Interface position in [m].
• nx1: Number of nodes on the left side.
• nx2: Number of nodes on the right side.
• dX1_dXN: Ratio of the first grid spacing to the last grid spacing.

Returns

• x_left: Array of nodes on the left side in [m].
• x_right: Array of nodes on the right side in [m].

make_dx_right(R, d1, n)

Constructs an array containing the right side dx.

Arguments

• R::Number: The common ratio to scale all dx.
• d1::Number: The initial value of dx.
• n::Integer: Number of dx elements.

Returns

• dx::Array{Float64,1}: Spatial distancenes for the right grid.

newton_solver(S, d, n, tol, max_iter, verbose)

Solves the non-linear equation F(x) = (1 - x^n) / (1-x) - S / d using the Newton method.

Arguments

• S::Float64: The value of S in the equation (scaling factor).
• d::Float64: The value of d in the equation (e.g. spatial distance).
• n::Int: The value of n in the equation (exponential factor).
• tol::Float64: The tolerance for convergence.
• max_iter::Int: The maximum number of iterations.
• verbose::Bool: Whether to print iteration information.

Returns

• x::Float64: The solution to the equation.

pchip(x, y, X)

Shape-preserving piecewise Cubic Hermite Interpolating Polynomial.

This function computes the shape-preserving piecewise cubic Hermite interpolating polynomial for the given data points (x, y).

Arguments

• x::Vector: The x-coordinates of the data points.
• y::Vector: The y-coordinates of the data points.
• X::Vector: The x-coordinates at which to evaluate the interpolating polynomial.

Returns

• P::Vector: The interpolated values at the points X.

regrid!(Fl_regrid, x_left, x_right, C_left, C_right, Ri, V_ip, nr, nmin, MRefin, verbose)

Regrid the grid and interpolate the composition profiles. Units may differ from SI units if non-dimensionalisation has been performed.

Arguments

• Fl_regrid::Int: Flag indicating whether to regrid or not.
• x_left::Vector: Vector of left spatial nodes in [m].
• x_right::Vector: Vector of right spatial nodes in [m].
• C_left::Vector: Vector of left composition values in [-].
• C_right::Vector: Vector of right composition values in [-].
• Ri::Vector: Radii [interface total length] in [m].
• V_ip::Float64: Velocity of the interface in [m/s].
• nr::Vector: Resolution.
• nmin::Int: Minimum grid size.
• MRefin::Int: Refinement factor.
• verbose::Bool: Whether to print additional information.

Returns

• x_left::Matrix: Matrix of left nodes in [m].
• x_right::Matrix: Matrix of right nodes in [m].
• C_left::Vector: Vector of left composition values in [-].
• C_right::Vector: Vector of right composition values in [-].
• dx1::Float64: Grid spacing at the left interface in [m].
• dx2::Float64: Grid spacing at the right interface in [m].
• nr::Vector: Updated resolution.

set_inner_bc_Lasaga!(Cl_i, beta, t, KD, D_r, D_l, D0, C_left, C_right, dx1, dx2, rho, L_g, R_g, nr)

Set inner boundary conditions for the special case of major element diffusion in a diffusion couple (Lasaga, 1983). Units may differ from SI units if non-dimensionalisation has been performed.

Arguments

• Cl_i::Float64: Initial composition on the left side in [-].
• beta::Float64: Variable from Lasagas semi-analytical solution.
• t::Float64: Time in [s].
• KD::Float64: Distribution coefficient.
• D_r::Float64: Diffusion coefficient on the right side in [m²/s].
• D_l::Float64: Diffusion coefficient on the left side in [m²/s].
• D0::Array{Float64}: Pre-exponential factor within the equation for the diffusion coefficient (D` at T0).
• C_left::Array{Float64}: Array of concentrations on the left side in [-].
• C_right::Array{Float64}: Array of concentrations on the right side in [-].
• dx1::Float64: Grid spacing on the left side in [m].
• dx2::Float64: Grid spacing on the right side in [m].
• rho::Array{Float64}: Array of densities in [kg/m³].
• L_g::Array{Float64}: Global left-hand side matrix.
• R_g::Array{Float64}: Global right-hand side vector.
• nr::Array{Int64}: Resolution.

Returns

• L_g::Array{Float64}: Updated global left-hand side matrix.
• R_g::Array{Float64}: Updated global right-hand side vector.
• ScF::Float64: Scaling factor.
• BC_left::Float64: Modelled left inner boundary condition (interface) in [-].
• BC_right::Float64: Modelled right inner boundary condition (interface) in [-].
• BC_left_Las::Float64: Left inner boundary condition following Lasaga (1983) in [-].
• BC_right_Las::Float64: Right inner boundary condition following Lasaga (1983) in [-].

set_inner_bc_flux!(L_g, R_g, KD, D_l, D_r, x_left, x_right, V_ip, rho, nr)

Set the inner boundary conditions at the interface using fluxes. Units may differ from SI units if non-dimensionalisation has been performed.

Arguments

• L_g::Matrix: The global matrix representing the system of equations (LHS).
• R_g::Vector: The global vector representing the right-hand side (RHS) of the system of equations.
• KD::Float64: The distribution coefficient.
• D_l::Float64: The diffusion coefficient on the left side in [m²/s].
• D_r::Float64: The diffusion coefficient on the right side in [m²/s].
• x_left::Vector: Left spatial nodes in [m].
• x_right::Vector: Right spatial nodes in [m].
• V_ip::Float64: The interface velocity in [m/s].
• rho::Vector: The density value in [kg/m³].
• nr::Vector: Resolution.

Returns

• L_g::Matrix: Updated LHS matrix.
• R_g::Vector: Updated RHS vector.
• ScF::Float64: The scale factor used to reduce the condition number.

set_inner_bc_mb!(L_g, R_g, dVolC, Mtot, KD, nr)

Set the inner boundary conditions at the interface using mass balance (MB). Units may differ from SI units if non-dimensionalisation has been performed.

Arguments

• L_g::Matrix: The global matrix representing the system of equations (LHS).
• R_g::Vector: The global vector representing the right-hand side (RHS) of the system of equations.
• dVolC::Vector: The volume change vector in [m³].
• Mtot::Float64: The total mass in [-].
• KD::Float64: The distribution coefficient.
• nr::Vector: Resolution.

Returns

• L_g::Matrix: Updated LHS matrix.
• R_g::Vector: Updated RHS vector.
• ScF::Float64: The scale factor used to reduce the condition number.

set_inner_bc_stefan!(L_g, R_g, C_left, C_right, nr)

Uses the Stefan boundary conditions to set the inner boundary conditions of the system of equations.

Arguments

• L_g::Matrix: The matrix representing the left-hand side of the system of equations.
• R_g::Vector: The vector representing the right-hand side of the system of equations.
• C_left::Vector: Composition of the left phase in [-].
• C_right::Vector: Composition of the right phase in [-].
• nr::Vector: Resolution of the model.

Returns

• L_g::Matrix: The modified matrix L_g with Stefan conditions applied.
• R_g::Vector: The modified vector R_g with Stefan conditions applied.
• ScF::Float64: Scaling factor for reducing the condition number of the matrix.

set_outer_bc!(BCout, L_g, R_g, C_left, C_right, ScF)

Set the outer boundary conditions (Dirichlet or Neumann) for the diffusion-advection problem. Units may differ from SI units if non-dimensionalisation has been performed.

Arguments

• BCout: An array indicating the type of boundary condition at the outer boundaries.
• L_g: The global left-hand side matrix of the diffusion-advection problem.
• R_g: The global right-hand side vector of the diffusion-advection problem.
• C_left: Composition vector of the left side in [-].
• C_right: Composition vector of the right side in [-].
• ScF: A scaling factor.

Details

• If BCout[1] is equal to 1, the left outer boundary condition is Dirichlet, otherwise it is Neumann.
• If BCout[2] is equal to 1, the right outer boundary condition is Dirichlet, otherwise it is Neumann.

Returns

• L_g: The updated global left-hand side matrix.
• R_g: The updated global right-hand side vector.

sinusoid_profile(C0, n, L, D, t, G)

Calculates a sinusoidal composition profile.

This function takes the initial composition C0, the number of sinusoidal modes n, the length of the system L, the diffusion coefficient D, the time t, and the amplitude G as input parameters. It calculates the composition profile at a given time t using the sinusoidal equation. Units may differ from SI units if non-dimensionalisation has been performed.

Arguments

• C0: Initial composition at position x at `t = 0.0 in [-].
• n: Number of sinusoidal modes.
• L: Length of the modelling domain in [m].
• D: Diffusion coefficient in [m²/s].
• t: Time in [s].
• G: Amplitude.

Returns

• C: Composition profile at time t in [-].

solve_soe(L_g, R_g, res)

Solves a system of equations.

Arguments

• L_g: The left-hand side matrix of the system of equations.
• R_g: The right-hand side vector of the system of equations.
• res: Resolution.

Returns

• C_left: The updated composition vector of the left side in [-].
• C_right: The updated composition vector of the right side in [-].

trapezoidal_integration(x, fx)

Trapezoidal integration of fx over x.

Arguments

• x: Array of x values.
• fx: Array of f(x) values.

Returns

• Int: The integrated value.

update_t_dependent_param!(D0, Di, Ea1, Ea2, KD_ar, R, T_ar, t_ar, t, t_tot)

Update the time dependent parameters D_l, D_r, KD, and T based on the given inputs. If Di = [-1.0 -1.0], the diffusion coefficient will be calculated based on the Arrhenius relation. Units may differ from SI units if non-dimensionalisation has been performed.

Arguments

• D0::Vector{Float64}: Pre-exponential factor within the calculation of the diffusion coefficient in [m²/s].
• Di::Vector{Float64}: Diffusion coefficients in [m²/s].
• Ea1::Float64: Activation energy for the left phase in [J/mol].
• Ea2::Float64: Activation energy for the right phase in [J/mol].
• KD_ar::Vector{Float64}: Array of distributioncoefficients.
• R::Float64: Gas constant in [J/(mol*K)].
• T_ar::Vector{Float64}: Array of temperatures in [K].
• t_ar::Vector{Float64}: Array of time in [s].
• t::Float64: Current time in [s].
• t_tot::Float64: Total time in [s].

Returns

• D_l::Float64: Updated diffusion coefficient for the left side in [m²/s].
• D_r::Float64: Updated diffusion coefficient for the right side in [m²/s].
• KD::Float64: Updated distribution coefficient.
• T::Float64: Updated temperature in [K].

update_t_dependent_param_simple!(D0, Di, Ea1, R, T_ar, t_ar, t, t_tot)

Update the dependent parameters D and T based on the given inputs. If Di = [-1.0], the diffusion coefficient will be calculated based on the Arrhenius relation. Units may differ from SI units if non-dimensionalisation has been performed.

Arguments

• D0: Pre-exponential factor within the calculation of the diffusion coefficient in [m²/s].
• Di: Initial diffusion coefficient in [m²/s].
• Ea1: Activation energy in [J/mol].
• R: Gas constant in [J/(mol*K)].
• T_ar: Array of temperatures in [K].
• t_ar: Array of time in [s].
• t: Current time in [s].
• t_tot: Total time in [s].

Returns

• D_l: Updated diffusion coefficient in [m²/s].
• T: Updated temperature in [K].

update_time!(t, dt, it, t_tot)

Update the time and time related variables t, dt, and it for a given total time t_tot. Units may differ from SI units if non-dimensionalisation has been performed.

Arguments

• t::Number: Current time in [s].
• dt::Number: Time step in [s].
• it::Integer: Time iterations.
• t_tot::Number: Total time in [s].

Returns

• t::Number: Updated time in [s].
• dt::Number: Updated time step in [s].
• it::Integer: Updated time iterations.