<reponame>dmitry-kabanov/fickettmodel<gh_stars>0
import numpy as np
import math
import matplotlib.pyplot as plt
from scipy.integrate import ode
from scipy.integrate import quad

def compute_u_for_given_lamda(lamda):
    # As we compute steady solution we start at x = 0 and move to the right,
    # that is, it's more convenient to consider detonation moving to the left.
    # Therefore, when we compute u, it should be negative.
    # u = -D - np.sqrt(D**2 - Q * lamda)
    # assert u <= 0.0, \
    #     'Steady u should be negative for detonation moving to the left'
    # return u
    return -D - math.sqrt(D_SQ - Q * lamda)

def compute_omega_for_given_lamda(lamda):
    u = compute_u_for_given_lamda(lamda)
    return K * (1 - lamda) * math.exp(THETA * (-u + Q * lamda))

def rhs(t, x):
    u, lamda = compute_state_at_x(x)
    
    return u + D

def compute_state_at_x(x):
    lamda = compute_lamda(x)
    lamda = np.minimum(lamda, 1.0)
    #assert 0.0 <= lamda <= 1.0, 'Lamda must be within range [0, 1]'
    u = compute_u_for_given_lamda(lamda)
    
    return u, lamda

def compute_lamda(x):
    dt = x/100.0
    y0, t0 = [0], 0
    r = ode(func_lamda).set_integrator('dopri5')
    r.set_initial_value(y0, t0)
    
    while r.successful() and r.t < x:
        r.integrate(r.t + dt)
        
    return r.y[0]

def func_lamda(x, lamda):
    omega = compute_omega_for_given_lamda(lamda)
    
    return omega / D

def compute_k():
    """ Compute $k = \int_0^{1/2} D / \omega d\lambda$ """
    rhs = lambda lamda: D / compute_omega_for_given_lamda(lamda)
    k, _ = quad(rhs, 0, 0.5)
    assert k > 0, 'Preexponent must be positive'
    return k

def compute_reaction_length():
    rhs = lambda lamda: D / compute_omega_for_given_lamda(lamda)
    length, _ = quad(rhs, 0, 1 - 1.0e-6)
    assert length > 0, 'Length must be positive'
    return length

Q = 2
D = np.sqrt(Q)
D_SQ = D**2
THETA = 2
K = 1
K = compute_k()
L = compute_reaction_length()

tfinal = 100.0
dt = 0.5

time = []
x = []

y0, t0 = [L], 0.0
r = ode(rhs).set_integrator('dopri5')
r.set_initial_value(y0, t0)

while r.successful() and r.t < tfinal:
    time.append(r.t + dt)
    r.integrate(r.t + dt)
    x.append(r.y)
    
plt.plot(x, time)
plt.show()

import numpy as np
import math
import matplotlib.pyplot as plt
from scipy.integrate import ode
from scipy.integrate import quad

def compute_u_for_given_lamda(lamda):
    # As we compute steady solution we start at x = 0 and move to the right,
    # that is, it's more convenient to consider detonation moving to the left.
    # Therefore, when we compute u, it should be negative.
    # u = -D - np.sqrt(D**2 - Q * lamda)
    # assert u <= 0.0, \
    #     'Steady u should be negative for detonation moving to the left'
    # return u
    return -D - math.sqrt(D_SQ - Q * lamda)

def compute_omega_for_given_lamda(lamda):
    u = compute_u_for_given_lamda(lamda)
    return K * (1 - lamda) * math.exp(THETA * (-u + Q * lamda))

def rhs(t, x):
    u, lamda = compute_state_at_x(x)
    
    return u + D

def compute_state_at_x(x):
    lamda = compute_lamda(x)
    lamda = np.minimum(lamda, 1.0)
    #assert 0.0 <= lamda <= 1.0, 'Lamda must be within range [0, 1]'
    u = compute_u_for_given_lamda(lamda)
    
    return u, lamda

def compute_lamda(x):
    dt = x/100.0
    y0, t0 = [0], 0
    r = ode(func_lamda).set_integrator('dopri5')
    r.set_initial_value(y0, t0)
    
    while r.successful() and r.t < x:
        r.integrate(r.t + dt)
        
    return r.y[0]

def func_lamda(x, lamda):
    omega = compute_omega_for_given_lamda(lamda)
    
    return omega / D

def compute_k():
    """ Compute $k = \int_0^{1/2} D / \omega d\lambda$ """
    rhs = lambda lamda: D / compute_omega_for_given_lamda(lamda)
    k, _ = quad(rhs, 0, 0.5)
    assert k > 0, 'Preexponent must be positive'
    return k

def compute_reaction_length():
    rhs = lambda lamda: D / compute_omega_for_given_lamda(lamda)
    length, _ = quad(rhs, 0, 1 - 1.0e-6)
    assert length > 0, 'Length must be positive'
    return length

Q = 2
D = np.sqrt(Q)
D_SQ = D**2
THETA = 2
K = 1
K = compute_k()
L = compute_reaction_length()

tfinal = 100.0
dt = 0.5

time = []
x = []

y0, t0 = [L], 0.0
r = ode(rhs).set_integrator('dopri5')
r.set_initial_value(y0, t0)

while r.successful() and r.t < tfinal:
    time.append(r.t + dt)
    r.integrate(r.t + dt)
    x.append(r.y)
    
plt.plot(x, time)
plt.show()
# As we compute steady solution we start at x = 0 and move to the right, # that is, it's more convenient to consider detonation moving to the left. # Therefore, when we compute u, it should be negative. # u = -D - np.sqrt(D**2 - Q * lamda) # assert u <= 0.0, \ #     'Steady u should be negative for detonation moving to the left' # return u #assert 0.0 <= lamda <= 1.0, 'Lamda must be within range [0, 1]' Compute $k = \int_0^{1/2} D / \omega d\lambda$