Modelling fiber lasers in Julia as a BVDAE problem

Specific Domains Modelling & Simulations
1

Hello,
I am trying to model a fiber laser using Julia. I’ve managed to simulate an amplifier which is an initial value DAE problem, now i need to extend the code to a BVDAE problem to handle different laser resonator architectures. I’ve read the relevant documentation and based on it tried to use the Ascher method however the solver does not converge. Code fails before that, when setting up boundary conditions (I think). Below is a working example, since the system is stiff I attach all variables. Thank you in advance for any help :slight_smile:

import OrdinaryDiffEq as ODE
import BoundaryValueDiffEq as bvp

c = 3e8  # [m/s]
h = 6.626e-34  # [Js]
wavl_s = 2036e-9  # [m]
wavl_p = 793e-9  # [m]
sigma_es = 2.5e-25  # 2.09e-25  # [m ^ 2]
sigma_ep = 6.02e-25  # 8.0e-25  # [m^2]
sigma_as = 1.25e-27  # 3.42e-27 [m^2]
sigma_ap = 8.41e-25  # 8.7e-25 #9.45e-25 [m^2]
k3101 = 1.0e-22  # [m^3/s]
k1310 = 0.084 * k3101  # [m^3/s]
a = 25.0e-6 / 2
b = 400e-6 / 2
A_core = pi * (a) ^ 2
A_clad = pi * (b) ^2
NA_cladding = 0.46
NA_core = 0.09
Vs = 2 * pi * (a / wavl_s) * NA_core
overlap_p = A_core/A_clad
overlap_s = 1 - exp((-2 * a^2) / (A_core^2))
const_s = wavl_s / (h * c * A_core)
const_p = wavl_p / (h * c * A_core)
pump_abs_dbm = 4.2
pump_abs = log(10) * pump_abs_dbm / 10
Ntot = (pump_abs) / (overlap_p * sigma_ap)
tau1 = 460.0e-6
tau3 = 60.0e-6
b31 = 0.7
b30 = 1 - b31

function DAElaser!(du,u,p,t)
    n1,n3,pf,pr,sf,sr = u
    n0 = Ntot - n1 - n3
    w03 = const_p * sigma_ap * n0 * (pf + pr) * overlap_p
    w30 = const_p * sigma_ep * n3 * (pf + pr) * overlap_p
    w10 = const_s * sigma_es * n1 * (sf + sr) * overlap_s
    w01 = const_s * sigma_as * n0 * (sf + sr) * overlap_s
    cr1 = k3101 * n3 * n0 - k1310 * n1^2
    #du[1] =  (n1 / tau1) + ( (b30 * n3) / tau3) + w10 + w30 - w03 - w01 - cr1 
    du[1] = ( (b31 * n3) / tau3) + 2 * cr1 + w01 - w10 - (n1 / tau1)
    du[2] = w03 - w30 - (n3 / tau3) - cr1
    du[3] = -overlap_p * ((sigma_ap*n0) - (sigma_ep*n3)) * pf
    du[4] = overlap_p * ((sigma_ap*n0) - (sigma_ep*n3)) * pr
    du[5] = overlap_s * ((sigma_es*n1) - (sigma_as*n0)) * sf
    du[6] = -overlap_s * ((sigma_es*n1) - (sigma_as*n0)) * sr
end

u0 = [0.05*Ntot, 0.9*Ntot, 30, 0.0, 1e-3, 0.0] 
zspan=(0.0, 6.5)
m_matrix = [0 0 0 0 0 0
            0 0 0 0 0 0
            0 0 1.0 0 0 0
            0 0 0 1.0 0 0
            0 0 0 0 1.0 0
            0 0 0 0 0 1.0]
hr = 0.99
oc=0.04

function bc!(res, u, p, t)
    res[1] = u[1] - u0[3]
    res[2] = u[3] - hr*u[4]
    res[3] = u[4] - oc*u[3]
    res[4] = u[1]
end

f = bvp.BVPFunction(DAElaser!,bc!; mass_matrix=m_matrix)
prob_bvp = bvp.BVProblem(f, u0, zspan)
z = [0.0, 0.0, zspan[end], zspan[end]]
solver = bvp.Ascher4(; zeta=z, jac_alg = bvp.BVPJacobianAlgorithm(bvp.AutoForwardDiff()))
sol = bvp.solve(prob_bvp, solver; dt=0.01)