Experiment_3_reactive_astrocyte_Part3

"""

Solving the metabolic system in a reactive astrocytic domain

"""

import numpy as np
from dolfin import *
from cutfem import *
from timeit import default_timer as timer
from Dimensionless_parameters import *

parameters["form_compiler"]["no-evaluate_basis_derivatives"] = False
parameters['allow_extrapolation'] = True
# Start timer
startime = timer()

date = time.strftime("%d%m%Y")

############################

#   READ FROM HDF5 File  #

############################

f = HDF5File(mpi_comm_world(),  './AD_results_'+ date + '/hdfFile', 'r')
bg_mesh = Mesh()
f.read(bg_mesh, "bg_mesh", False)
W = FunctionSpace(bg_mesh, "CG", 1)
level_set = Function(W)
f.read(level_set, 'level_set')
f.close()


# read the volume from part 1
volume = np.load('./AD_results_' + date + '/Volume.npy')


mesh = CutFEMTools_fictitious_domain_mesh(bg_mesh,level_set,0,0)

L = np.load('./AD_results_' + date + '/L.npy')

# Define the sum of ATP+ADP inside the cell
sum_atp_adp = 3.2
# define the dimensionless parameter alpha
alpha = 0.16


# define the dimensionless diffusion concentration for each species
D_a_bar, D_b_bar, D_c_bar, D_d_bar, D_e_bar, D_f_bar = dimensionless_diffusion_coeff()

# Dimensionless parameters
t_c, f_bar, lac_bar, beta_a, beta_b, beta_c, beta_d, gamma_b, gamma_c, gamma_d, gamma_e, mu_e, mu_g, csi_b, csi_c, csi_e, tau_b, tau_c = dim_param(L, alpha, sum_atp_adp, 1.0)

# store the dimensionless time
np.save('./AD_results_' + date + '/t_c.npy', t_c)

# Define final dimensionless time
T = 90/t_c # final time
print('final time', T)
num_step = 250 # number of time step
dt = T / num_step
k = Constant(dt)

# cutfem parameters
beta = 10.0
gamma =0.1

# compute normals
n = FacetNormal(mesh)
h = CellSize(mesh)

# Compute Mesh -Levelset intersection and corresponding marker
mesh_cutter = MeshCutter(mesh,level_set)

# Diminishing quadrature points
q_degree = 3

# Define new measures associated with the interior domains and facets
dx = Measure("dx")[mesh_cutter.domain_marker()]
dS = Measure("dS")[mesh_cutter.interior_facet_marker(0)]
dxq = dc(0, metadata={"num_cells": 1, 'quadrature_degree': q_degree})
dsq = dc(1, metadata={"num_cells": 1})

dxc = dx(0, metadata={'quadrature_degree': q_degree}) + dxq

# Finite Element space for the concentration
T = FunctionSpace(mesh,"P",1)
V0 = MixedFunctionSpace([T, T, T, T, T, T])

# Define test functions
(v_1, v_2, v_3, v_4, v_5, v_6) = TestFunctions(V0)

# Define Trial functions which must be Functions instead of Trial Functions cause the pb is non linear
u = Function (V0)

# Define the initial condition of concetrations
a_0 = 0.0
b_0 = 1.6 / sum_atp_adp
c_0 = 1.6 / sum_atp_adp
d_0 = 0.0
e_0 = 0.0
f_0 = 0.0


u_0 = Expression(('a_0', 'b_0', 'c_0','d_0', 'e_0', 'f_0'), a_0=a_0, b_0=b_0, c_0=c_0, d_0=d_0, e_0=e_0, f_0=f_0, degree=1)

u_n = interpolate(u_0, V0)

a, b, c, d, e, f = split(u)
a_n, b_n, c_n, d_n, e_n, f_n = split(u_n)

# Reaction sites 

from Enzymes import Gauss_expression_3d


# From Part 2 read the coordinates of the reaction sites
list_of_enzymes = np.load('./AD_results_' + date + '/enzymes_coordinates.npy')

M = len(list_of_enzymes[0])


coordinate_enzymes_hxk = list_of_enzymes[0:3]
coordinate_enzymes_pyrk = list_of_enzymes[3:6]
coordinate_enzymes_ldh = list_of_enzymes[6:9]

# Read coordinates of mitochondria
mito = np.load('Images/AD_mito4cutfem_scaled_sigma.npy')
mito_points = mito.T

# Read the adaptive parameter eta
eta_hxk, eta_pyrk, eta_ldh, eta_mito = np.load('./AD_results_' + date + '/eta.npy')

# Variance dimensionless
sigma = 1./L

def Sum_Gaussian_Mito(M, coordinate_enzymes_test, sigma):
    Gaussian = 0
    for i in range(M):
        Gaussian += Gauss_expression_3d(coordinate_enzymes_test[0,i], coordinate_enzymes_test[1,i],coordinate_enzymes_test[2,i], sigma[i])
    return(Gaussian)

def Sum_Gaussian(M, coordinate_enzymes_test, sigma):
    if M == 1:
        Gaussian = Gauss_expression_3d(coordinate_enzymes_test[0], coordinate_enzymes_test[1],coordinate_enzymes_test[2], sigma)
    else:
        Gaussian = 0
        for i in range(M):
            Gaussian += Gauss_expression_3d(coordinate_enzymes_test[0,i], coordinate_enzymes_test[1,i],coordinate_enzymes_test[2,i], sigma)
    return(Gaussian)

Gaussian_hxk = Sum_Gaussian(M, coordinate_enzymes_hxk, sigma)
Gaussian_pyrk = Sum_Gaussian(M, coordinate_enzymes_pyrk, sigma)
Gaussian_ldh = Sum_Gaussian(M, coordinate_enzymes_ldh, sigma)
Gaussian_mito = Sum_Gaussian_Mito(len(mito.T[0]), mito_points[:-1], mito_points[-1])



# Reaction rates 
k_hxk = Constant(0.0619)
k_pyrk = Constant(1.92)
k_ldh = Constant(0.719)
k_mito = Constant(0.0813)

K_act = Constant(0.169)

# Spatial reaction rates
K_hxk = Gaussian_hxk/eta_hxk  * k_hxk *  Constant(volume)
K_pyrk = Gaussian_pyrk/eta_pyrk  * k_pyrk *  Constant(volume)
K_ldh = Gaussian_ldh/eta_ldh  * k_ldh *  Constant(volume)
K_mito = Gaussian_mito/eta_mito  * k_mito *  Constant(volume)


################################################################

# Define the source term

# read the subdomain volume computed in Part 2
subdomain_volume = np.load('./AD_results_' + date + '/subdomain_volume_glc.npy')

influx =  f_bar * volume /subdomain_volume

radius_influx = 0.0307

x_infl, y_infl, z_infl =  0.313278, 0.520833, 0.0302059
x_infl2, y_infl2, z_infl2 = 0.104167, 0.0, 0.0815185
x_infl3, y_infl3, z_infl3 = 0.364583, 0.572917, 0.0402746

f_1 = Expression('(pow(x[0] - x_0,2) + pow(x[1] - y_0,2) + pow(x[2]- z_0,2)) < (r * r) ? influx : 0', influx=influx, r=radius_influx, x_0 = x_infl , y_0 = y_infl, z_0 = z_infl, degree=1)
f_1 += Expression('(pow(x[0] - x_0,2) + pow(x[1] - y_0,2) + pow(x[2]- z_0,2)) < (r * r) ? influx : 0', influx=influx, r=radius_influx, x_0 = x_infl2 , y_0 = y_infl2, z_0 = z_infl2, degree=1)
f_1 += Expression('(pow(x[0] - x_0,2) + pow(x[1] - y_0,2) + pow(x[2]- z_0,2)) < (r * r) ? influx : 0', influx=influx, r=radius_influx, x_0 = x_infl3, y_0 = y_infl3, z_0 = z_infl3, degree=1)

# Define Lac

subdomain_outflux_volume = np.load('./AD_results_' + date + '/subdomain_volume_lac.npy')

outflux = lac_bar * volume/ subdomain_outflux_volume

radius_outflux = 0.035

x_out, y_out, z_out = 0.923611, 1., 0.0203796
x_out2, y_out2, z_out2 =  1., 0.909722, 0.0135864
x_out3, y_out3, z_out3 = 0.506944, 0.993056, 0.0543457
x_out4, y_out4, z_out4 =0.989583, 0.510142, 0.0402746

eta_f = Expression('(pow(x[0] - x_0,2) + pow(x[1] - y_0,2) + pow(x[2]- z_0,2)) < (r * r) ? outflux : 0', outflux = outflux, r=radius_outflux, x_0 = x_out , y_0 = y_out, z_0 = z_out, degree=1)
eta_f += Expression('(pow(x[0] - x_0,2) + pow(x[1] - y_0,2) + pow(x[2]- z_0,2)) < (r * r) ? outflux : 0', outflux = outflux, r=radius_outflux, x_0 = x_out2 , y_0 = y_out2, z_0 = z_out2, degree=1)
eta_f += Expression('(pow(x[0] - x_0,2) + pow(x[1] - y_0,2) + pow(x[2]- z_0,2)) < (r * r) ? outflux : 0', outflux = outflux, r=radius_outflux, x_0 = x_out3, y_0 = y_out3, z_0 = z_out3, degree=1)
eta_f += Expression('(pow(x[0] - x_0,2) + pow(x[1] - y_0,2) + pow(x[2]- z_0,2)) < (r * r) ? outflux : 0', outflux = outflux, r=radius_outflux, x_0 = x_out4, y_0 = y_out4, z_0 = z_out4, degree=1)

#  Weak formulation dimensionless system

F = ((a - a_n) / k) * v_1 * dxc \
    + D_a_bar * dot(grad(a), grad(v_1)) * dxc + K_hxk * beta_a * a * b**2 * v_1 * dxc \
    + ((b - b_n) / k) * v_2 * dxc  \
    + D_b_bar * dot(grad(b), grad(v_2)) * dxc + 2 * K_hxk * beta_b * a * b**2 * v_2 * dxc - 2 * K_pyrk * gamma_b * d * c**2 * v_2 * dxc - 28 * K_mito * csi_b * e * c**28 * v_2 * dxc + K_act * tau_b * b * v_2 * dxc\
    + ((c - c_n) / k)*v_3 * dxc \
    + D_c_bar * dot(grad(c), grad(v_3)) * dxc - 2 * K_hxk * beta_c * a * b**2 * v_3 * dxc  + 2 * K_pyrk * gamma_c * d * c**2 * v_3 * dxc - K_act *tau_c * b * v_3 * dxc + 28 * K_mito * csi_c * e * c**28 * v_3 * dxc\
    + ((d - d_n) / k)*v_4 * dxc\
    + D_d_bar * dot(grad(d),grad(v_4)) * dxc - 2 * K_hxk * beta_d * a * b**2 * v_4 * dxc + K_pyrk * gamma_d * d * c**2 * v_4 * dxc\
    + ((e - e_n) / k)*v_5 * dxc\
    + D_e_bar * dot(grad(e),grad(v_5)) * dxc - K_pyrk * gamma_e * d * c**2 * v_5 * dxc + K_ldh * mu_e * e * v_5 * dxc + K_mito * csi_e * e * c**28 * v_5 * dxc\
    + ((f - f_n) / k)*v_6 * dxc\
    + D_f_bar * dot(grad(f),grad(v_6)) * dxc - K_ldh * mu_g * e * v_6 * dxc + eta_f * f * v_6 * dxc\
    - f_1 * v_1 * dxc


F += avg(gamma) * avg(h) * D_a_bar * dot(jump(grad(a), n), jump(grad(v_1), n)) * dS(1) + avg(gamma) * avg(h) * D_b_bar *dot(jump(grad(b), n), jump(grad(v_2), n)) * dS(1)\
    + avg(gamma) * avg(h) * D_c_bar * dot(jump(grad(c), n), jump(grad(v_3), n)) * dS(1) + avg(gamma) * avg(h) * D_d_bar * dot(jump(grad(d), n), jump(grad(v_4), n)) * dS(1)\
    + avg(gamma) * avg(h) * D_e_bar * dot(jump(grad(e), n), jump(grad(v_5), n)) * dS(1) + avg(gamma) * avg(h) * D_f_bar * dot(jump(grad(f), n), jump(grad(v_6), n)) * dS(1)



#Create VTK files for visualization output
vtkfile_a = File('AD_results_' + date + '/a/a.pvd')
vtkfile_b = File('AD_results_' + date + '/b/b.pvd')
vtkfile_c = File('AD_results_' + date + '/c/c.pvd')
vtkfile_d = File('AD_results_' + date + '/d/d.pvd')
vtkfile_e = File('AD_results_' + date + '/e/e.pvd')
vtkfile_f = File('AD_results_' + date + '/f/f.pvd')

# Compute Jacobian
J = derivative(F, u)

# Fictitious domain
composite_mesh = CompositeMesh()
composite_mesh.add(mesh)

V = CompositeFunctionSpace(composite_mesh)
V.add(V0);
V.build();

# Constrain dofs outside
FidoTools_compute_constrained_dofs(V, mesh_cutter)

a = FidoForm(V, V)
form_a = create_dolfin_form(J)
a.add(form_a, mesh_cutter)

L = FidoForm(V)
form_L = create_dolfin_form(F)
L.add(form_L, mesh_cutter);

# Space for the solution
cutmesh0 = CutFEMTools_physical_domain_mesh(mesh, mesh_cutter.cut_cells(0), mesh_cutter.domain_marker(), 0);
V0Phys = FunctionSpace(cutmesh0, "CG", 1);

a_inter = Function(V0Phys);
b_inter = Function(V0Phys);
c_inter = Function(V0Phys);
d_inter = Function(V0Phys);
e_inter = Function(V0Phys);
f_inter = Function(V0Phys);

# Time stepping
t = [0.0]

# Function that compute the average concentration inside the astrocyte
def compute_average_concentration(conc, V):

    int_conc = conc * dxc

    form_conc_mean = create_dolfin_form(int_conc)

    composite_form_conc_mean = CompositeForm(V)
    composite_form_conc_mean.add(form_conc_mean)

    cut_cells = mesh_cutter.cut_cells(0)

    quadrature = Quadrature(cut_cells.type().cell_type(),cut_cells.geometry().dim(),order=2)

    composite_form_conc_mean.cut_form(0).set_quadrature(0, quadrature);
    composite_form_conc_mean.cut_form(0).set_cut_mesh(0, cut_cells);
    composite_form_conc_mean.cut_form(0).add_single_parent_mesh_id(0, 0);

    mean_conc = composite_assemble(composite_form_conc_mean)

    return(mean_conc)

def compute_LAC_efflux(bg_mesh, level_set, subdom2comp):

    r = radius_outflux

    class Sub_Efflux(SubDomain):
        def inside(self, x, on_boundary):
            return ((pow(x[0] - x_out, 2) + pow(x[1] - y_out, 2) + pow(x[2] - z_out, 2)) <= (r * r))

    class Sub_Efflux2(SubDomain):
        def inside(self, x, on_boundary):
                return ( (pow(x[0] - x_out2, 2) + pow(x[1] - y_out2, 2) + pow(x[2] - z_out2, 2)) < (r * r))

    class Sub_Efflux3(SubDomain):
        def inside(self, x, on_boundary):
                return ( (pow(x[0] - x_out3, 2) + pow(x[1] - y_out3, 2) + pow(x[2] - z_out3, 2)) < (r * r) )

    class Sub_Efflux4(SubDomain):
        def inside(self, x, on_boundary):
                return ( (pow(x[0] - x_out4, 2) + pow(x[1] - y_out4, 2) + pow(x[2] - z_out4, 2)) < (r * r) )

    # Create mesh

    mesh = CutFEMTools_fictitious_domain_mesh(bg_mesh, level_set, 0, 0)

    # Compute Mesh -Levelset intersection and corresponding marker
    mesh_cutter = MeshCutter(mesh, level_set)

    marker_test = mesh_cutter.domain_marker()

    Sub_Efflux().mark(marker_test, 50)
    Sub_Efflux2().mark(marker_test, 50)
    Sub_Efflux3().mark(marker_test, 50)
    Sub_Efflux4().mark(marker_test, 50)

    # file = File("./test_domain_marker.pvd")
    # file << marker_test

    # Define new measures associated with the interior domains and facets
    dx = Measure("dx", domain=mesh)[marker_test]
    dxq = dc(0, metadata={"num_cells": 1})  # cut cells inside the Circle

    dx_eff = dx(50) + dxq

    V0 = FunctionSpace(mesh, "Lagrange", 1)

    # Fictitious domain
    composite_mesh = CompositeMesh()
    composite_mesh.add(mesh)

    V = CompositeFunctionSpace(composite_mesh)
    V.add(V0);
    V.build();

    # Constrain dofs outside
    FidoTools_compute_constrained_dofs(V, mesh_cutter)

    u_e_V0 = project(subdom2comp, V0)

    psi0 = u_e_V0 * dx_eff

    form_psi0 = create_dolfin_form(psi0)

    composite_form_psi0 = CompositeForm(V)
    composite_form_psi0.add(form_psi0)

    area_ls_cutfem = composite_assemble(composite_form_psi0)

    return(area_ls_cutfem)

# Empty list to save the values

list_a =[]
list_b =[]
list_c =[]
list_d =[]
list_e = []
list_f = []

time_list = []


# add the time = 0
time_list.append(t[0])

mean_a_ = compute_average_concentration(a_n, V)
mean_b_ = compute_average_concentration(b_n, V)
mean_c_ = compute_average_concentration(c_n, V)
mean_d_ = compute_average_concentration(d_n, V)
mean_e_ = compute_average_concentration(e_n, V)
mean_f_ = compute_average_concentration(f_n, V)


list_a.append(mean_a_/ volume)
list_b.append(mean_b_/ volume)
list_c.append(mean_c_/ volume)
list_d.append(mean_d_/ volume)
list_e.append(mean_e_/ volume)
list_f.append(mean_f_/ volume)

efflux_lac = []

lac_eff = compute_LAC_efflux(bg_mesh, level_set, f_n)
efflux_lac.append(lac_eff)

# Parameter for solver
Nmax = 50
abs_tol = 1.0e-12
rel_tol = 1.0e-07


# Initial residual
initial_residual = composite_assemble(L)
absolute0 = initial_residual.norm('l2')

for i in range(num_step):
    print("Timestep", i)
    print("Time", t)

    n = 1
    while n < Nmax:

        A = composite_assemble(a)
        b = composite_assemble(L)

        uc = CompositeFunction(V)

        solve(A, uc.vector(), -b, 'mumps')

        u.vector().axpy(1.0, uc.part(0).vector())

        residual = composite_assemble(L)

        # absolute residual
        absolute = residual.norm('l2')

        # relative residual
        relative = absolute / absolute0

        if absolute < abs_tol or relative < rel_tol:
            break
        else:
            n += 1

    _a, _b, _c, _d, _e, _f = u.split()

    ## To save the vtk uncomment this section.
    # 
    ## save  GLC
    #a_inter.interpolate(_a);
    ## save  ATP
    #b_inter.interpolate(_b);
    ## save  ADP
    #c_inter.interpolate(_c);
    ## save  GLY
    #d_inter.interpolate(_d);
    ## save  PYR
    #e_inter.interpolate(_e);
    ## save  LAC
    #f_inter.interpolate(_f);
    #
    #vtkfile_a << (a_inter, t[0])
    #vtkfile_b << (b_inter, t[0])
    #vtkfile_c << (c_inter, t[0])
    #vtkfile_d << (d_inter, t[0])
    #vtkfile_e << (e_inter, t[0])
    #vtkfile_f << (f_inter, t[0])
    
    u_n._assign(u)

    t[0] = t[0] + dt

    time_list.append(t[0])
    
    #Compute average concentration inside astrocyte of each species
    mean_a_ = compute_average_concentration(_a, V)
    mean_b_ = compute_average_concentration(_b, V)
    mean_c_ = compute_average_concentration(_c, V)
    mean_d_ = compute_average_concentration(_d, V)
    mean_e_ = compute_average_concentration(_e, V)
    mean_f_ = compute_average_concentration(_f, V)


    list_a.append(mean_a_/ volume)
    list_b.append(mean_b_/ volume)
    list_c.append(mean_c_/ volume)
    list_d.append(mean_d_/ volume)
    list_e.append(mean_e_/ volume)
    list_f.append(mean_f_/ volume)

    # Compute the LAC concentration in the subregions where the efflux is defined
    lac_eff = compute_LAC_efflux(bg_mesh, level_set, _f)
    efflux_lac.append(lac_eff)

# stop time
aftersolve = timer()
tottime = aftersolve-startime
print('Final time', tottime)

# Create a single list with all the solutions
list_of_list = [list_a, list_b, list_c, list_d, list_e, list_f, time_list]

# save using numpy
print(list_a, list_b, list_c, list_d, list_e, list_f, time_list)
np.save('./AD_results_' + date + '/AD_3Dastro_list' + pathology, np.asarray(list_of_list))

np.save('./AD_results_' + date + '/lac_efflux' + pathology, np.asarray(efflux_lac))