## Supplimentaty material for the Paper ##

The code below solves the following equation set. 

$$
\frac{2}{\pi} A \bigg[\text{erf}\bigg(\frac{y+W}{\sqrt{4\alpha_{Th}\,x}}\bigg) -  \text{erf}\bigg(\frac{y-W}{\sqrt{4\alpha_{Th}\,x}}\bigg)\bigg ]\\
\cdot e^{(- B x)}- D = 0
$$

$$
-B D +\frac{2 A }{\pi}\left(\frac{\left(- \frac{W}{2} + \frac{y}{2}\right) e^{- \frac{\left(- W + y\right)^{2}}{4 \alpha_{Th} \,x}}}{\sqrt{\pi} x \sqrt{\alpha_{Th}\,x}} - \frac{\left(\frac{W}{2} + \frac{y}{2}\right) e^{- \frac{\left(W + y\right)^{2}}{4 \alpha_{Th} x}}}{\sqrt{\pi} x \sqrt{\alpha_{Th} \,x}}\right) e^{(- B x)}=0
$$

with

 $ A=  (\gamma C_D^\circ + C_A^\circ) $, $ B = \frac{\pi^2}{4M^2}\alpha_{Tv}$, and $D =(C_A^\circ + C^\ast)$

(The symbols name are in input cell)
 

**The solution steps:**

1. Provide model input
2. Obtain estimate of $L_{max}$ using Liedl et al. (2011)
3. Adjust the starting value for $W_{max}$ and $X_{wmax}$
4. Iteratively solve the equation set using Newton-Krylov (LGMRES method, see below for details)


**First we import the Python library** 


```python
import numpy as np
from sympy import diff, symbols
import matplotlib.pyplot as plt
import scipy.optimize as so
import scipy.special as ss
from scipy import optimize
from scipy.special import erf
```

### Step 1: the model input: 


```python
M = 2
W = 5
ath = 0.036
atv = 0.0018
Ca = 7.68
Cd = 17.3
Ct = 0.01
ga = 3.5
```


### Step II: Obtaining $L_{max}$ using Liedl et al. (2011) ###

The expression provided in Liedl et al. (2011) is:

$$
\text{erf}\Bigg(\frac{W}{\sqrt{4\alpha_{Th}L_{max}}}\Bigg)\,\text{e}^{-\alpha_{Tv}L_{max} \big(\frac{\pi}{2M}\big)^2} =\frac{\pi}{4}\frac{\gamma C^*+C_{A}^\circ}{\gamma C_{D}^\circ+C_{A}^\circ}
$$

$L_{max}$ is implicit and requires iterative solution. We use Newton-Raphson (NR) method to obtain $L_{max}$.

The code below solves NR using estimate for differential and using exact FD


**The main NR routine**


```python
def f_lm(x): # Liedl 2011- f(x) = 0
    return erf(W/(np.sqrt(4*ath*x)))*np.exp(-atv*x*(np.pi/(2*M))**2)-(np.pi/4)*((ga*Ct+Ca)/(ga*Cd + Ca)) 

def df_lm(x): # FD estimate df
    h = 1e-4
    return (f_lm(x+h) - f_lm(x-h))/(2*h)

def df_lm1(x): # exact df
    return -W/(2*x*np.sqrt(np.pi)*np.sqrt(atv*x))*np.exp(-(W**2)/(4*atv*x))*np.exp(-(np.pi**2*atv*x)/(4*M**2)) \
        - (1/(4*M**2))*np.pi**2*atv*np.exp((-np.pi**2*atv*x)/(4*M**2))*erf(W/(2*np.sqrt(atv*x)))  

def NR(x): # Newton Raphson simulation using FD
    iterat = 0
    tol = 1e-06
    h = f_lm(x)/df_lm(x)
    print("Iter. Nr ","", "Lmax","       " , " Residual")
    print('-'*35)   
    
    while abs(h)>= tol:
        h = f_lm(x)/df_lm(x) 
        x = x - h
        iterat+= 1 
        #print(" ",iterat,  "    " , "%.4f "% x, "   " , "%.2e " %  h) 
    return x
        
def NR1(x): # Newton Raphson simulation using exact df
    iterat = 0
    tol = 1e-06
    h = f_lm(x)/df_lm1(x)
    
    while abs(h)>= tol:
        h = f_lm(x)/df_lm1(x)
        x = x - h
        iterat+= 1
    return x 
```

**Starting value for NR method**

The starting value is based on Liedl et al. (2005), which provides explicit $L_{max}$ from 2D set-up. 


```python
ma_1 = -1/(np.pi*(ath/W**2)*np.log(1-(0.25*np.pi*((ga*Ct+Ca)/(ga*Cd+Ca)))))
ma_2 = -2/(np.pi**2*(atv/M**2))*np.log(0.25*np.pi*((ga*Ct+Ca)/(ga*Cd+Ca))) 
ma_3 = 4/np.pi**2*(M**2/atv)*np.log((4/np.pi)*((ga*Cd+ Ca)/(ga*Ct+ Ca)))

ma_x0 = np.minimum(np.maximum(ma_1, ma_2), ma_3)

min_x0 = np.minimum(-1/(np.pi*(ath/W**2) *np.log(1-(0.25*np.pi*((ga*Ct+Ca)/(ga*Cd+Ca))))), \
                    -2/(np.pi**2*(atv/M**2))*np.log(0.25*np.pi*((ga*Ct+Ca)/(ga*Cd+Ca))))

if ma_x0 == ma_3:
    x0 = ma_x0
else:
    x0 = min_x0 
```


```python
lmax = NR(x0) # simulations using exact differential
print(lmax) 
```

    Iter. Nr   Lmax          Residual
    -----------------------------------
    1333.098505604007


### Step III: Starting value of $W_{max}$ and $X_{wmax}$ ###

- We use $L_{max}$ from above to obtain the starting value of $X_{wmax}$
- For $W_{max}$ the starting value is arbitarily considered (See accompanied table to choose appropriate starting value)


```python
A = (ga*Cd+Ca)
B = (np.pi**2*atv)/(4*M**2)
D = Ca+Ct
x_s = lmax/ 4
print(x_s)
y_s = W * 1.3
print(y_s)
st = [x_s, y_s] 
```

    333.27462640100174
    6.5


### Step IV: Solution of $W_{max}$ and $X_{wmax}$ using Newton-Krylov (NK) non-linear system solver using LGMRES method. ###

The code used below is described in detailed here: 
https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.newton_krylov.html

First we set-up our function


```python
def f(x): 
    return [(2*A)/(np.pi)*(ss.erf((x[1]+W)/(np.sqrt(4*ath*x[0])))-ss.erf((x[1]-W)/(np.sqrt(4*ath*x[0]))))*np.exp(-B*x[0])-D, 
            -B*D + (2*A)/np.pi*np.exp(-B*x[0])*((((-W/2 + x[1]/2)*np.exp(-(x[1]-W)**2/(4*ath*x[0])))/(x[0]*np.sqrt(np.pi*ath*x[0]))) - (((W/2 + x[1]/2)*np.exp(-(x[1]+W)**2/(4*ath*x[0])))/(x[0]*np.sqrt(np.pi*ath*x[0]))))]
```

**The main NK routine**

Few solver conditions:

tolerance: 6e-6
Preconditioner for Jacobian: None
Maximum iteration: Not set.

The solution format (see below) is ($X_{wmax}$, $W_{max}$) and $L_{max}$


```python
sol = so.newton_krylov(f, st, method= 'lgmres', inner_maxiter= 10 ,outer_k= 3, verbose = 1, line_search= 'wolfe')
sol, lmax
#print(lmax)
```

    0:  |F(x)| = 4.38903; step 1
    1:  |F(x)| = 0.784995; step 1
    2:  |F(x)| = 0.0390243; step 1
    3:  |F(x)| = 0.00483961; step 1
    4:  |F(x)| = 0.00483965; step 1
    5:  |F(x)| = 0.00483965; step 1
    6:  |F(x)| = 0.00483965; step 1
    7:  |F(x)| = 0.00483965; step 1
    8:  |F(x)| = 0.00483965; step 1
    9:  |F(x)| = 7.02837; step 1
    10:  |F(x)| = 1.60188; step 1
    11:  |F(x)| = 0.144272; step 1
    12:  |F(x)| = 0.00483912; step 1
    13:  |F(x)| = 0.00483965; step 1
    14:  |F(x)| = 8.40867; step 1
    15:  |F(x)| = 2.08032; step 1
    16:  |F(x)| = 0.228314; step 1
    17:  |F(x)| = 0.00330359; step 1
    18:  |F(x)| = 0.951433; step 1
    19:  |F(x)| = 0.0557367; step 1
    20:  |F(x)| = 0.00238468; step 1
    21:  |F(x)| = 0.00227957; step 0.0437571
    22:  |F(x)| = 0.00218419; step 0.0415658
    23:  |F(x)| = 0.00209682; step 0.0397667
    24:  |F(x)| = 0.00201605; step 0.038311
    25:  |F(x)| = 0.00194074; step 0.0371648
    26:  |F(x)| = 0.00186995; step 0.0363023
    27:  |F(x)| = 0.00180289; step 0.0357017
    28:  |F(x)| = 0.0017389; step 0.0353461
    29:  |F(x)| = 0.00167741; step 0.0352211
    30:  |F(x)| = 0.00161794; step 0.0353125
    31:  |F(x)| = 0.00156011; step 0.035613
    32:  |F(x)| = 0.00150356; step 0.036116
    33:  |F(x)| = 0.001448; step 0.0368208
    34:  |F(x)| = 0.00139317; step 0.0377298
    35:  |F(x)| = 0.00133886; step 0.0388509
    36:  |F(x)| = 0.00128486; step 0.0401969
    37:  |F(x)| = 0.00123098; step 0.0417865
    38:  |F(x)| = 0.00117707; step 0.0436476
    39:  |F(x)| = 0.00112296; step 0.0458156
    40:  |F(x)| = 0.0010685; step 0.048338
    41:  |F(x)| = 0.00102453; step 0.0512743
    42:  |F(x)| = 0.00102118; step 0.0547092
    43:  |F(x)| = 0.00101499; step 0.0587499
    44:  |F(x)| = 0.00100574; step 0.0635482
    45:  |F(x)| = 0.000993113; step 0.0693056
    46:  |F(x)| = 0.000976687; step 0.076312
    47:  |F(x)| = 0.000955872; step 0.0850005
    48:  |F(x)| = 0.0010017; step 0.1711
    49:  |F(x)| = 0.000977867; step 0.180974
    50:  |F(x)| = 0.000911354; step 0.207964
    51:  |F(x)| = 0.000807632; step 0.271372
    52:  |F(x)| = 0.000647892; step 0.426889
    53:  |F(x)| = 0.000324548; step 1
    54:  |F(x)| = 3.14415e-08; step 1





    (array([553.79859438,  11.04398976]), 1333.098505604007)



6db27f04145043da426a8996ae6868c65d48691e
ipynb
Jupyter Notebook
Wmax_computation.ipynb
prabhasyadav/Wmax
206233c25f4245a345859abb364117a43d5b5a50
CC-BY-4.0
Wmax_computation.ipynb
prabhasyadav/Wmax
206233c25f4245a345859abb364117a43d5b5a50
CC-BY-4.0
Wmax_computation.ipynb
prabhasyadav/Wmax
206233c25f4245a345859abb364117a43d5b5a50
CC-BY-4.0
Qwen/Qwen-72B
1. YES
2. YES