Skip to content
Permalink
master
Switch branches/tags

Name already in use

A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?

DrugReleaseSystemModel/implant_model.py

Go to file
 
 
Cannot retrieve contributors at this time
1929 lines (1552 sloc) 84.9 KB
Raw Blame
Open in GitHub Desktop
# =============================================================================
# Import Libraries
# =============================================================================
import numpy as np
import matplotlib.pyplot as plt
import math as ma
from scipy.sparse import SparseEfficiencyWarning
from scipy.sparse.linalg import spsolve
from scipy import sparse
from scipy.integrate import odeint
import time
from time import process_time
import matplotlib.tri as tri
import warnings
# =============================================================================
# Set Model Conditions
# =============================================================================
warnings.simplefilter('ignore', SparseEfficiencyWarning)
# Set Model to in Vitro or in Vivo
# in Vivo requires a defined rb (smallest radius), pH, and Diff. of Tissue
# in Vivo considering setting Swell = 0 to turn off swelling
Model = 'in Vitro' # 'in Vivo'
MW = 29000
Lt = 23#23
pH_bath = 7.8#7.4
# in vivo parameters
rb_inVivo = 2.285/1000
pH_bath_inVivo = 7.4
Diff_Tissue_Drug_inVivo = 1e-10 #m^2/day
# Toggles to set degredation, swelling, and drug diff. coef. opt.
Degr = 1
Opt = 0
ke1 = 0.0039#42
ke2 = 6e-5#0
kr1 = 1* 6.65e-5#4453125e-05
kr2 = 0.075#17
if Degr == 0: ke1 = ke2 = kr1 = kr2 = 0
Beta = 8.1
# =============================================================================
# Material Properties
# =============================================================================
rho_p = 1340
rho_w = 997
rho_s = 1030#854
rho_d = 1600
Mw_w = 18
Mw_s = 99.133 #89.13 #NMP
Mw_d = 332.31
Mw_col = (90.08 + 76.05 - 1.008*2)*4
Mw_cm = (90.08 + 76.05 - 1.008*2)
CH0 = 1000*(10**(-pH_bath))
M0 = 72;
m = 4;
Col0 = 0
Cm0 = 0
Pmass_rat = 0.4
Smass_rat = 0.6
Vol = 5.01507e-8
Mass = Vol/(0.39/rho_p + 0.6/rho_s)
#rb0 = (Vol/4*3/np.pi)**(1/3)
Polymass = 0.39 * Mass
Solmass = 0.6 * Mass
DrugTotal = 0.01 * Mass
Pmass_rat = 0.4
Smass_rat = 0.6
Vol = 5e-8#5.01507e-8
Mass = Vol/(0.39/rho_p + 0.6/rho_s + 0.01/rho_d)
rb0 = (Vol/4*3/np.pi)**(1/3)
Polymass = 0.39 * Mass
Solmass = 0.6 * Mass
DrugTotal = 0.01 * Mass
fDrug = 0
MWpolyunit = (72.068 + 90.1008)/2
MWcol = (MWpolyunit)*4
mm = MW/(MWpolyunit)
SwellM = np.array([
[0,1/3,1,2,3,4,5,6,7,8,9,10]
,[1.14,1.113839116,1.125772077,1.319071372,1.402670617,1.456901822,1.519157263,1.550285064,1.530332265,1.531844287,1.538057059,1.547226309]
,[1.03719,1.087972141,1.094369578,1.157885346,1.205755558,1.229106885,1.273335458,1.321601701,1.343241483,1.470949449,1.470949449,1.470949449]
,[0.96715,0.970342272,0.988546752,1.040669341,1.075871936,1.127231584,1.166579349,1.218221573,1.232711366,1.317084591,1.275656165,0.997112764]
,[1.0,1.123140564,1.126275839,1.209823608,1.27311905,1.331267788,1.36737552,1.427182005,1.482569932,1.531785663,1.519902324,1.4893521]
,[1.0,0.981655651,0.98410626,0.972763299,1.008192926,1.037573234,1.042276866,1.092698255,1.106364359,1.182551823,1.412915932,1.299325116]
,[1.0,1.040536374,1.051502505,1.069818144,1.094057022,1.134256645,1.135462148,1.206042205,1.197495562,1.22378213,1.303643527,1.327323056]
])
expp = np.array([-13.9051,0.168151,0.0289776])
def Diff(Mw_drug,Mw_poly):
return np.exp(expp[0] * Mw_drug**expp[1] * Mw_poly**expp[2])
def StokeEin(Mw,rho):
Rad = ((3*Mw)/(4*np.pi*(6.022e23)*rho))**(1/3)
return (8.314*300)/((6e23)*3*np.pi*(1e-3)*Rad)
Dpolymer_H2O =(3e-10)*(24*3600)
Dpolymer_Sol = Diff(Mw_s,MW/1000) *60*60*24
Dpolymer_Col = Diff(Mw_col,MW/1000) *60*60*24
Dpolymer_Cm = Diff(Mw_cm,MW/1000) *60*60*24
Dpolymer_Drug = Diff(Mw_d,MW/1000) *(24*3600)
if Opt == 1:
if MW == 53000: Dpolymer_Drug = (2.5e-18) *(24*3600)
if MW == 29000: Dpolymer_Drug = (3.5e-18) *(24*3600)
if MW == 15000: Dpolymer_Drug = (16e-18) *(24*3600)
Dpore_H2O = (2e-9)*(24*3600)
Dpore_Sol = StokeEin(Mw_s,rho_s) *(24*60*60)
Dpore_Col = StokeEin(Mw_col,rho_p) *(24*3600)
Dpore_Cm = StokeEin(Mw_cm,rho_p) *(24*3600)
Dpore_Drug = StokeEin(Mw_d,rho_d) *(24*3600)
Dpolymer_CH = Dpolymer_H2O
Dpore_CH = Dpore_H2O
Dpore_Drug_Tissue = Dpore_Drug
if Model == 'in Vitro':
if MW == 53000:
Si = 3
Ti = 5
rp = 0.019031/1000/2
rs = 0.390/1000
thick = 0.02/1000
rps = 0.1*rs
rb = SwellM[Si,0]* 2.285/1000
rb0 = SwellM[Si,0]* 2.285/1000
Data = np.array([[0,0.000694444,0.020833333,0.041666667,0.083333333,0.166666667,0.25,0.333333333,1,2,3,4,5,6,7,8,9,10,11,12,13,14,16,18,21]
,[0,1.119142804,9.262339789,11.98741229,15.53883212,20.24333537,22.5071373,24.45746377,30.61285371,31.01585248,32.50954267,33.77538118,34.88172184,35.29867581,36.50942933,38.51772902,40.07357018,41.98452285,45.2975403,47.94281328,51.54607499,55.08537035,63.08351632,67.7997305,70.4767212]
,[0,0.133326478,0.967230761,0.93595068,1.083056477,1.95273396,2.077156874,1.555949026,2.307256901,3.807566178,1.811207162,1.663381142,1.482157325,1.420050663,0.9641997,1.719178327,1.960206207,2.456818305,3.308653552,5.374767878,5.308919469,6.23759649,6.897256065,6.502648253, 5.622977768]])
DataMn = np.array([[1,3,7,10,14,17,21]
,[1,1.007763134,0.640007222,0.624480953,0.355569597,0.312908467,0.309044954]
,[0.017810236,0.144697465,0.060462983,0.053449922,0.010278687,0.060640608,0.150238972]])
DataML = np.array([[1,3,5,7,8,10,12,14,16,18,21]
,[46.8134715,47.20207254,46.94300518,47.4611399,41.24352332,40.33678756,34.89637306,35.02590674,32.8238342,27.12435233,26.60621762]
,[1.683937824,1.295336788,1.943005181,1.813471503,2.720207254,1.424870466,1.554404145,1.03626943,1.424870466,1.813471503,2.461139896]])
elif MW == 29000:
Si = 2
Ti = 1
rp = 0.010/1000/2
rs =0.390/1000
rps = 0.1*rs
thick = 0.01/1000
rb = SwellM[Si,0]* 2.285/1000
rb0 = SwellM[Si,0]* 2.285/1000
Data = np.array([[0,0.000694444,0.020833333,0.041666667,0.083333333,0.166666667,0.25,0.333333333,1,2,3,4,5,6,7,8,9,10,11,12,13,14,16,18,21]
,[0,0.802561858,7.226856607,8.440377965,9.989583769,12.85666748,14.50869886,16.13060466,22.67370573,26.16872192,29.4456466,32.51134381,34.71584224,37.06137629,39.60130539,44.04048026,47.83582988,53.32927363,57.31533667,60.66826522,66.10551599,69.33759084,76.24546257,76.85323346,75.03531962]
,[0,0.455660094,0.853294764,0.994756217,1.2508843,1.429279326,1.421216848,1.559346725,2.864824775,1.550817231,1.345094475,0.716603102,1.185592383,1.695622715,2.514804205,4.716818952,5.798447622,6.724925016,7.176424076,9.078293132,7.139356076,6.451263526,4.1746822,6.425503411,6.870846248]])
DataMn = np.array([[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0]])
DataML = np.array([[1,3,5,7,8,10,12,14,16,18,21]
,[47.79220779,44.93506494,41.55844156,42.72727273,37.66233766,36.62337662,35.97402597,32.85714286,28.05194805,25.32467532,22.98701299]
,[0.909090909,1.558441558,1.948051948,2.857142857,2.727272727,1.558441558,2.337662338,1.038961039,0.909090909,0.779220779,2.467532468]])
elif MW == 15000:
Si = 1
Ti = 1
rp = 0.1/1000/2
rs = 0.147/1000
rps = 0
thick = 0.04/1000
rb = SwellM[Si,0]* 2.285/1000
rb0 = SwellM[Si,0]* 2.285/1000
Data = np.array([[0,0.000694444,0.020833333,0.041666667,0.083333333,0.166666667,0.25,0.333333333,1,2,3,4,5,6,7,8,9,10,11,12,13,14,16,18,21]
,[0,0.462476336,5.123648484,7.433042562,10.02455512,13.51578389,15.74610437,17.48820752,24.79401366,31.76581687,36.68680098,43.08394241,48.22577652,53.4191935,58.2543482,62.43889017,66.95240196,73.19940025,78.12550524,81.06608078,83.23052045,83.60509647,81.03100642,79.59390574,79.76481325]
,[0,0.078657889,0.663636899,1.115984807,1.546849241,1.82003899,1.988922704,1.678892329,1.724597036,2.615812986,1.709123185,2.75052848,3.110626217,4.041749624,3.637894973,3.511931844,3.287644934,4.29795183,4.180872898,2.562052363,3.368844522,3.122991375,4.130334619,4.815340276,4.380351231]])
DataMn = np.array([[1,3,7,10,14,17,21]
,[1,0.751796346,0.461691644,0.414268117,0.351221515,0.304342024,0.318723055]
,[0.054732257,0.06828421,0.004970172,0.007388712,0.033398182,0.016940949,0.090968699]])
DataML = np.array([[1,3,5,7,8,10,12,14,16,18,21]
,[45.51813472,40.20725389,38.13471503,33.86010363,33.47150259,29.19689119,25.69948187,21.29533679,19.22279793,15.85492228,13.52331606]
,[2.202072539,1.683937824,2.07253886,1.683937824,1.813471503,0.647668394,2.720207254,1.813471503,0.777202073,0.647668394,0.777202073]])
elif MW == 22000:
Si = 4
Ti = 3
rp = 0.1/1000/2
rb = 2.21/1000
rs = 0.285/1000
thick = 0.04/1000
rps = 0
rb = 2.567/1000
rb0 = 2.567/1000
elif Model == 'in Vivo':
Si = 1
Ti = 1
rp = 0.1/1000/2
rs = 0.147/1000
rps = 0
thick = 0.04/1000
rb = SwellM[Si,0]* rb_inVivo
rb0 = rb_inVivo
pH_bath = pH_bath_inVivo
Dpore_Drug_Tissue = Diff_Tissue_Drug_inVivo
nx = 201
nttt = 21
ntt = 101
Nt = 2
Lx = 6/2/1000
dx = Lx/nx
dx0 = Lx/nx
dt = dx**2/(1e-8)
Nt = 2
LT = Lt
TitraM = np.array([
[0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0] #%WaterContent
,[0.017202142,0.052222222,0.59832664,0.691097724,0.921659973,0.974431058,1.147215529,1.044431058,0.901659973,0.958875502]
,[0.005612,0.016823,0.08972,0.120561,0.072897,0.075708,0.145802,0.131775,0.120561,0.092523]
,[0.003262,0.215883,0.871488,1.028035,1.013547,1.032717,1.037855,1.014956,0.997656,1.042054] #3A
,[0.011212,0.056068,0.039252,0.056075,0.028038,0.028037,0.016822,0.025234,0.014018,0.019626]
,[0.014469,0.580369,0.882703,1.061687,1.033173,0.976636,0.987381,0.998126,0.989252,0.941112] #4A
,[0.011215,0.098131,0.047663,0.044867,0.019626,0.042057,0,0.036441,0.019626,0.019626]
,[0.008876,0.642051,1.011675,0.927108,1.094848,1.021481,0.964958,0.972892,1.025701,1.005605] #4.5A
,[0.005607,0.064486,0.095327,0.014019,0.058871,0.042049,0.019626,0.011215,0.030841,0.056075]
,[0.01168,0.54392,0.9556,0.991586,1.047192,1.052336,1.00421,0.992519,1.042523,1.033642] #7E
,[0.011215,0.064485,0.061682,0.022429,0.011215,0.044866,0.030841,0.100935,0.070093,0.036448]])
Titr_fits = np.zeros((3,9))
Titr_fits = np.array([[0,8,0,33,0,50,0,60,0,116],
[0, 0.17923357,0,0.140819,0,0.113849,0,0.118412,0,0.117123],
[0,4.85877588,0,8.15931,0,5.52442,0,43.8569,0,7.25775]])
TitraMV = np.copy(TitraM)
TitraMV[0,:] = 0.6*TitraM[0,:]
WperV = np.linspace(0,0.6,1000)
TitraFV = np.zeros(len(WperV))
TitraFV[:] = 1 + (0-1)/(1+np.power(WperV[:]/Titr_fits[1,Ti],Titr_fits[2,Ti]))
def TFunc(phi_W):
return 1 + (0-1)/(1+np.power(phi_W[:]/Titr_fits[1,Ti],Titr_fits[2,Ti]))
def TFunc8_33(phi_W,MW):
return (abs(33-MW/1000)/(33-8) * (1+(0-1)/(1+np.power(phi_W[:]/Titr_fits[1,1],Titr_fits[2,1])))
+ abs(8-MW/1000)/(33-8) * (1+(0-1)/(1+np.power(phi_W[:]/Titr_fits[1,3],Titr_fits[2,3]))))
def TFunc33_50(phi_W,MW):
return (abs(50-MW/1000)/(50-33) * (1+(0-1)/(1+np.power(phi_W[:]/Titr_fits[1,3],Titr_fits[2,3])))
+ abs(33-MW/1000)/(50-33) * (1+(0-1)/(1+np.power(phi_W[:]/Titr_fits[1,5],Titr_fits[2,5]))))
def TFunc50_60(phi_W,MW):
return (abs(60-MW/1000)/(60-50) * (1+(0-1)/(1+np.power(phi_W[:]/Titr_fits[1,5],Titr_fits[2,5])))
+ abs(50-MW/1000)/(60-50) * (1+(0-1)/(1+np.power(phi_W[:]/Titr_fits[1,7],Titr_fits[2,7]))))
def TFunc15(phi_W):
return (abs(33-15)/(33-8) * (1+(0-1)/(1+np.power(phi_W[:]/Titr_fits[1,1],Titr_fits[2,1])))
+ abs(8-15)/(33-8) * (1+(0-1)/(1+np.power(phi_W[:]/Titr_fits[1,3],Titr_fits[2,3]))))
def TFunc29(phi_W):
return (abs(33-29)/(33-8) * (1+(0-1)/(1+np.power(phi_W[:]/Titr_fits[1,1],Titr_fits[2,1])))
+ abs(8-29)/(33-8) * (1+(0-1)/(1+np.power(phi_W[:]/Titr_fits[1,3],Titr_fits[2,3]))))
def TFunc53(phi_W):
return (abs(60-52)/(60-50) * (1+(0-1)/(1+np.power(phi_W[:]/Titr_fits[1,5],Titr_fits[2,5])))
+ abs(50-52)/(60-50) * (1+(0-1)/(1+np.power(phi_W[:]/Titr_fits[1,7],Titr_fits[2,7]))))
def PlotW():
plt.figure(14)
plt.plot(TitraM[0,:],TitraM[Ti,:],'ok',label='Experimental Data')
plt.plot(WperV/0.6,TFunc(WperV))
Xext_1D = np.zeros((nx,Nt))
Xc_1D = np.zeros((nx,Nt))
Vc = 4.19e-24
p = 1e-2
etaA = 6.022e23
lam = 1
Xmax = 0.655
Xc0 = 0
#Estimating on and off rates
# ~10^7 M^-1s^-1 -- 10^8 M^-1s^-1 for a small molecule colliding with a protein
# ~10^6 M^-1s^-1 -- 10^7 M^-1s^-1 for protein-protein
# Kd * kon == koff
# 10^-6.4 * 10^7 = koff
# kon =10^7
kon = 200 #1e7 # can be increased back to 1e7
Kdpol = (1.4*10**-4)*(60*60*24)
Kdol = (1.4*10**-4)*(60*60*24)
Kdm = (1.4*10**-4)*(60*60*24)
Kddrug = (2.6*10**-4)*(60*60*24)
Kdw = (1*10**-14)*(60*60*24)
kondrug = kon
koffdrug = kondrug*Kddrug
konpol = kon
koffpol = konpol*Kdpol
konol = kon
koffol = konol*Kdol
konm = kon
koffm = konm*Kdm
konw = kon
koffw = konw*Kdw
Xcend_1D = np.zeros((nx,Nt)); Maskfit_1D = np.zeros((nx,Nt)); invMaskfit_1D = np.zeros((nx,Nt))
X = np.zeros((nx));Y = np.zeros((nx));CH_1D = np.zeros((nx,Nt));CH2O_1D = np.zeros((nx,Nt))
COH_1D = np.zeros((nx,Nt));Rs_1D = np.zeros((nx,Nt));Rrs_1D = np.zeros((nx,Nt));Res_1D = np.zeros((nx,Nt))
Col_1D = np.zeros((nx,Nt));Colcoo_1D = np.zeros((nx,Nt));Cm_1D = np.zeros((nx,Nt));Cmcoo_1D = np.zeros((nx,Nt));
Cpol_1D = np.zeros((nx,Nt));Cpolcoo_1D = np.zeros((nx,Nt));DDrug_1D = np.zeros((nx,Nt));
DCH_1D = np.zeros((nx,Nt));DCol_1D = np.zeros((nx,Nt));DCm_1D = np.zeros((nx,Nt));
Drug_1D = np.zeros((nx,Nt));Drugcoo_1D = np.zeros((nx,Nt))
pinchCol_1D = np.zeros((nx,ntt));pinchColcoo_1D = np.zeros((nx,ntt));pinchCm_1D = np.zeros((nx,ntt))
pinchCmcoo_1D = np.zeros((nx,ntt));pinchCpol_1D = np.zeros((nx,ntt));pinchCpolcoo_1D = np.zeros((nx,ntt))
pinchDrug_1D = np.zeros((nx,ntt));pinchDrugcoo_1D = np.zeros((nx,ntt));pinchCH2O_1D = np.zeros((nx,ntt))
pinchCOH_1D = np.zeros((nx,ntt));pinchCH_1D = np.zeros((nx,ntt))
DrugRelease_1D = np.zeros((nx,Nt));VporeDrug_1D = np.zeros((nx,Nt));VporeCol_1D = np.zeros((nx,Nt));
VporeCm_1D = np.zeros((nx,Nt));Vpore_1D = np.zeros((nx,Nt));Mn_1D = np.zeros((nx,Nt));
Ce_1D = np.zeros((nx,Nt));Cend_1D = np.zeros((nx,Nt));Rol_1D = np.zeros((nx,Nt));
Rm_1D = np.zeros((nx,Nt));ProfileDrug = np.zeros(Nt);ProfileXc = np.zeros(Nt);ProfileMn = np.zeros(Nt);
Mask_1D = np.ones(nx);invMask_1D = np.ones(nx);tisMask_1D = np.ones(nx);invtisMask_1D = np.ones(nx)
CH2O_1D = np.zeros((nx,Nt));DCH2O_1D = np.zeros((nx,Nt));CSol_1D = np.zeros((nx,Nt))
DCSol_1D = np.zeros((nx,Nt));CeSwell_1D = np.ones(Nt);TimeT = np.zeros(Nt)
DegWeight = np.zeros((nx,Nt));DegWeightAcid = np.zeros(nx);DegWeightSciss = np.zeros(nx)
x = np.linspace(0,Lx,nx)
dX = np.zeros(Nt)
dX[:] = dx0 #*SwellM[Si,0]
phi_SM = np.zeros((nx,Nt))
phi_PM = np.zeros((nx,Nt))
phi_WM = np.zeros((nx,Nt))
Rb = np.zeros(Nt)
dX = np.zeros(Nt)
Rb[:] = rb0 *SwellM[Si,0]
dX[:] = dx
TimeR = np.zeros(Nt)
TimeT = np.zeros(Nt)
Variables = (
Mask_1D,invMask_1D,tisMask_1D,invtisMask_1D,Maskfit_1D,invMaskfit_1D,
Drug_1D,Drugcoo_1D,Col_1D,Colcoo_1D,CH_1D,COH_1D,CH2O_1D,Cpol_1D,Cpolcoo_1D,Cmcoo_1D,Cm_1D,
Ce_1D,Cend_1D,
DDrug_1D,DCol_1D,DCH_1D,DCm_1D,
Res_1D,Rrs_1D,Xc_1D,Xext_1D,Xcend_1D,Rol_1D,Rm_1D,
Mn_1D,VporeCol_1D,VporeDrug_1D,Vpore_1D,CSol_1D,dt,DegWeightSciss,DegWeightAcid,DegWeight,rb,dx,rb0,CeSwell_1D,
phi_SM,phi_PM,phi_WM,Rb,dX,TimeR,TimeT,DrugRelease_1D)
def ImplantModel(Variables):
# =============================================================================
# Initialize Variables
# =============================================================================
(Mask_1D,invMask_1D,tisMask_1D,invtisMask_1D,Maskfit_1D,invMaskfit_1D,
Drug_1D,Drugcoo_1D,Col_1D,Colcoo_1D,CH_1D,COH_1D,CH2O_1D,Cpol_1D,Cpolcoo_1D,Cmcoo_1D,Cm_1D,
Ce_1D,Cend_1D,
DDrug_1D,DCol_1D,DCH_1D,DCm_1D,
Res_1D,Rrs_1D,Xc_1D,Xext_1D,Xcend_1D,Rol_1D,Rm_1D,
Mn_1D,VporeCol_1D,VporeDrug_1D,Vpore_1D,CSol_1D,dt,DegWeightSciss,DegWeightAcid,DegWeight,rb,dx,rb0,CeSwell_1D,
phi_SM,phi_PM,phi_WM,Rb,dX,TimeR,TimeT,DrugRelease_1D) = Variables
# =============================================================================
# Create Mask weighting function
# =============================================================================
(Mask_1D,invMask_1D,tisMask_1D,invtisMask_1D
) = Maskfunc(
rb,rp,rs,x,Lx,Mask_1D,invMask_1D,tisMask_1D,invtisMask_1D,dx)
# =============================================================================
# Initialize Functions
# =============================================================================
(Drug_1D,Cpol_1D,Col_1D,Cm_1D,CH_1D,DDrug_1D,DCol_1D,DCm_1D,fDrug,Maskfit_1D,invMaskfit_1D,Cpol0,Ce0,CSol_1D,CH2O_1D,DCSol_1D,DCH2O_1D,Sol_conc,Ce0, phi_SM,phi_PM,phi_WM
)= Initialize(
Mask_1D,invMask_1D,Maskfit_1D,invMaskfit_1D,invtisMask_1D,tisMask_1D,CSol_1D,CH2O_1D,rb,dx, phi_SM,phi_PM,phi_WM)
A = np.zeros((nx,nx))
TimeT_mem = np.copy(TimeT)
TimeR_mem = np.copy(TimeR)
Ce_1D_mem = np.copy(Ce_1D)
Cend_1D_mem = np.copy(Cend_1D)
Res_1D_mem = np.copy(Res_1D)
Rrs_1D_mem = np.copy(Rrs_1D)
Rol_1D_mem = np.copy(Rol_1D)
Rm_1D_mem = np.copy(Rm_1D)
Cm_1D_mem = np.copy(Cm_1D)
Cmcoo_1D_mem = np.copy(Cmcoo_1D)
Col_1D_mem = np.copy(Col_1D)
Colcoo_1D_mem = np.copy(Colcoo_1D)
Cpol_1D_mem = np.copy(Cpol_1D)
Cpolcoo_1D_mem = np.copy(Cpolcoo_1D)
Drug_1D_mem = np.copy(Drug_1D)
Drugcoo_1D_mem = np.copy(Drugcoo_1D)
CH_1D_mem = np.copy(CH_1D)
DDrug_1D_mem = np.copy(DDrug_1D)
DCSol_1D_mem = np.copy(DCSol_1D)
DCm_1D_mem = np.copy(DCm_1D)
DCol_1D_mem = np.copy(DCol_1D)
DCH2O_1D_mem = np.copy(DCH2O_1D)
DCH_1D_mem = np.copy(DCH_1D)
VporeCol_1D_mem = np.copy(VporeCol_1D)
VporeDrug_1D_mem = np.copy(VporeDrug_1D)
Vpore_1D_mem = np.copy(Vpore_1D)
Xc_1D_mem = np.copy(Xc_1D)
Xcend_1D_mem = np.copy(Xcend_1D)
Xext_1D_mem = np.copy(Xext_1D)
Maskfit_1D_mem = np.copy(Maskfit_1D)
invMaskfit_1D_mem = np.copy(invMaskfit_1D)
CH2O_1D_mem = np.copy(CH2O_1D)
COH_1D_mem = np.copy(COH_1D)
CSol_1D_mem = np.copy(CSol_1D)
DegWeight_mem = np.copy(DegWeight)
Mn_1D_mem = np.copy(Mn_1D)
phi_SM_mem = np.copy(phi_SM)
phi_WM_mem = np.copy(phi_WM)
phi_PM_mem = np.copy(phi_PM)
DrugRelease_1D_mem = np.copy(DrugRelease_1D)
# =============================================================================
# =============================================================================
# =============================================================================
# # # Begin Time Iteration
# =============================================================================
# # =============================================================================
# =============================================================================
skp = 1
t = 0; n = 1; nstart = n#int(ntdeg/skpe)
cc = 0; c = 0
DrugStart = np.sum(invtisMask_1D*(Drug_1D_mem[:,-1]+Drugcoo_1D_mem[:,-1]))
rb = SwellM[Si,cc] * rb0
while TimeT[t] <= (Lt ):
drugleft = np.sum(invtisMask_1D*(Drug_1D[:,t]+Drugcoo_1D[:,t]))/DrugStart
if TimeT[t] > c:
print(f'Day = {int(TimeT[t]):<4} \tDrugLeft = {ma.trunc(drugleft*1000)/1000:0.3f}')
#print('Day =',int(TimeT[t]),'\tDrug Left =',ma.trunc(drugleft*1000)/1000)
c += 1
dt = 0.5* dx**2 / (np.sum(Mask_1D*DDrug_1D[:,t])/np.sum(Mask_1D))
if dt <= 1e-3: dt = 1e-3
# if drugleft < 0.1: dt = 0.1
count = n%skp
if count == 0:
TimeR_mem = np.append(TimeR_mem,TimeR[t+1])
TimeT_mem = np.append(TimeT_mem,TimeT[t+1])
Ce_1D_mem = np.append(Ce_1D_mem,Ce_1D[:,:t+1],axis=1)
Cend_1D_mem = np.append(Cend_1D_mem,Cend_1D[:,:t+1],axis=1)
Res_1D_mem = np.append(Res_1D_mem,Res_1D[:,:t+1],axis=1)
Rrs_1D_mem = np.append(Rrs_1D_mem,Rrs_1D[:,:t+1],axis=1)
Rol_1D_mem = np.append(Rol_1D_mem,Rol_1D[:,:t+1],axis=1)
Rm_1D_mem = np.append(Rm_1D_mem,Rm_1D[:,:t+1],axis=1)
Cm_1D_mem = np.append(Cm_1D_mem,Cm_1D[:,:t+1],axis=1)
Cmcoo_1D_mem = np.append(Cmcoo_1D_mem,Cmcoo_1D[:,:t+1],axis=1)
Col_1D_mem = np.append(Col_1D_mem,Col_1D[:,:t+1],axis=1)
Colcoo_1D_mem = np.append(Colcoo_1D_mem,Colcoo_1D[:,:t+1],axis=1)
Cpol_1D_mem = np.append(Cpol_1D_mem,Cpol_1D[:,:t+1],axis=1)
Cpolcoo_1D_mem = np.append(Cpolcoo_1D_mem,Cpolcoo_1D[:,:t+1],axis=1)
Drug_1D_mem = np.append(Drug_1D_mem,Drug_1D[:,:t+1],axis=1)
Drugcoo_1D_mem = np.append(Drugcoo_1D_mem,Drugcoo_1D[:,:t+1],axis=1)
CH_1D_mem = np.append(CH_1D_mem,CH_1D[:,:t+1],axis=1)
DDrug_1D_mem = np.append(DDrug_1D_mem,DDrug_1D[:,:t+1],axis=1)
DCSol_1D_mem = np.append(DCSol_1D_mem,DCSol_1D[:,:t+1],axis=1)
DCm_1D_mem = np.append(DCm_1D_mem,DCm_1D[:,:t+1],axis=1)
DCol_1D_mem = np.append(DCol_1D_mem,DCol_1D[:,:t+1],axis=1)
DCH2O_1D_mem = np.append(DCH2O_1D_mem,DCH2O_1D[:,:t+1],axis=1)
DCH_1D_mem = np.append(DCH_1D_mem,DCH_1D[:,:t+1],axis=1)
VporeCol_1D_mem = np.append(VporeCol_1D_mem,VporeCol_1D[:,:t+1],axis=1)
VporeDrug_1D_mem = np.append(VporeDrug_1D_mem,VporeDrug_1D[:,:t+1],axis=1)
Vpore_1D_mem = np.append(Vpore_1D_mem,Vpore_1D[:,:t+1],axis=1)
Xc_1D_mem = np.append(Xc_1D_mem,Xc_1D[:,:t+1],axis=1)
Xcend_1D_mem = np.append(Xcend_1D_mem,Xcend_1D[:,:t+1],axis=1)
Xext_1D_mem = np.append(Xext_1D_mem,Xext_1D[:,:t+1],axis=1)
Maskfit_1D_mem = np.append(Maskfit_1D_mem,Maskfit_1D[:,:t+1],axis=1)
invMaskfit_1D_mem = np.append(invMaskfit_1D_mem,invMaskfit_1D[:,:t+1],axis=1)
CH2O_1D_mem = np.append(CH2O_1D_mem,CH2O_1D[:,:t+1],axis=1)
COH_1D_mem = np.append(COH_1D_mem,COH_1D[:,:t+1],axis=1)
CSol_1D_mem = np.append(CSol_1D_mem,CSol_1D[:,:t+1],axis=1)
DegWeight_mem = np.append(DegWeight_mem,DegWeight[:,:t+1],axis=1)
Mn_1D_mem = np.append(Mn_1D_mem,Mn_1D[:,:t+1],axis=1)
phi_SM_mem = np.append(phi_SM_mem,phi_SM[:,:t+1],axis=1)
phi_WM_mem = np.append(phi_WM_mem,phi_WM[:,:t+1],axis=1)
phi_PM_mem = np.append(phi_PM_mem,phi_PM[:,:t+1],axis=1)
DrugRelease_1D_mem = np.append(DrugRelease_1D_mem,DrugRelease_1D[:,:t+1],axis=1)
# =============================================================================
# What was old will be new again
# =============================================================================
TimeT[t] = TimeT[t+1]; TimeT[t+1] = 0
TimeR[t] = TimeR[t+1]; TimeR[t+1] = 0
Ce_1D[:,t] = Ce_1D[:,t+1]; Ce_1D[:,t+1] = np.zeros(nx)
Cend_1D[:,t] = Cend_1D[:,t+1]; Cend_1D[:,t+1] = np.zeros(nx)
Res_1D[:,t] = Res_1D[:,t+1]; Res_1D[:,t+1] = np.zeros(nx)
Rrs_1D[:,t] = Rrs_1D[:,t+1]; Rrs_1D[:,t+1] = np.zeros(nx)
Rol_1D[:,t] = Rol_1D[:,t+1]; Rol_1D[:,t+1] = np.zeros(nx)
Rm_1D[:,t] = Rm_1D[:,t+1]; Rm_1D[:,t+1] = np.zeros(nx)
Cm_1D[:,t] = Cm_1D[:,t+1]; Cm_1D[:,t+1] = np.zeros(nx)
Cmcoo_1D[:,t] = Cmcoo_1D[:,t+1]; Cmcoo_1D[:,t+1] = np.zeros(nx)
Col_1D[:,t] = Col_1D[:,t+1]; Col_1D[:,t+1] = np.zeros(nx)
Colcoo_1D[:,t] = Colcoo_1D[:,t+1]; Colcoo_1D[:,t+1] = np.zeros(nx)
Cpol_1D[:,t] = Cpol_1D[:,t+1]; Cpol_1D[:,t+1] = np.zeros(nx)
Cpolcoo_1D[:,t] = Cpolcoo_1D[:,t+1]; Cpolcoo_1D[:,t+1] = np.zeros(nx)
Drug_1D[:,t] = Drug_1D[:,t+1]; Drug_1D[:,t+1] = np.zeros(nx)
Drugcoo_1D[:,t] = Drugcoo_1D[:,t+1]; Drugcoo_1D[:,t+1] = np.zeros(nx)
CH_1D[:,t] = CH_1D[:,t+1]; CH_1D[:,t+1] = np.zeros(nx)
DDrug_1D[:,t] = DDrug_1D[:,t+1]; DDrug_1D[:,t+1] = np.zeros(nx)
DCSol_1D[:,t] = DCSol_1D[:,t+1]; DCSol_1D[:,t+1] = np.zeros(nx)
DCm_1D[:,t] = DCm_1D[:,t+1]; DCm_1D[:,t+1] = np.zeros(nx)
DCol_1D[:,t] = DCol_1D[:,t+1]; DCol_1D[:,t+1] = np.zeros(nx)
DCH2O_1D[:,t] = DCH2O_1D[:,t+1]; DCH2O_1D[:,t+1] = np.zeros(nx)
DCH_1D[:,t] = DCH_1D[:,t+1]; DCH_1D[:,t+1] = np.zeros(nx)
VporeCol_1D[:,t] = VporeCol_1D[:,t+1]; VporeCol_1D[:,t+1] = np.zeros(nx)
VporeDrug_1D[:,t] = VporeDrug_1D[:,t+1]; VporeDrug_1D[:,t+1] = np.zeros(nx)
Vpore_1D[:,t] = Vpore_1D[:,t+1]; Vpore_1D[:,t+1] = np.zeros(nx)
Xc_1D[:,t] = Xc_1D[:,t+1]; Xc_1D[:,t+1] = np.zeros(nx)
Xcend_1D[:,t] = Xcend_1D[:,t+1]; Xcend_1D[:,t+1] = np.zeros(nx)
Xext_1D[:,t] = Xext_1D[:,t+1]; Xext_1D[:,t+1] = np.zeros(nx)
Maskfit_1D[:,t] = Maskfit_1D[:,t+1]; Maskfit_1D[:,t+1] = np.zeros(nx)
invMaskfit_1D[:,t] = invMaskfit_1D[:,t+1]; invMaskfit_1D[:,t+1] = np.zeros(nx)
CH2O_1D[:,t] = CH2O_1D[:,t+1]; CH2O_1D[:,t+1] = np.zeros(nx)
COH_1D[:,t] = COH_1D[:,t+1]; COH_1D[:,t+1] = np.zeros(nx)
CSol_1D[:,t] = CSol_1D[:,t+1]; CSol_1D[:,t+1] = np.zeros(nx)
DegWeight[:,t] = DegWeight[:,t+1]; DegWeight[:,t+1] = np.zeros(nx)
Mn_1D[:,t] = Mn_1D[:,t+1]; Mn_1D[:,t+1] = np.zeros(nx)
phi_SM[:,t] = phi_SM[:,t+1]; phi_SM[:,t+1] = np.zeros(nx)
phi_WM[:,t] = phi_WM[:,t+1]; phi_WM[:,t+1] = np.zeros(nx)
phi_PM[:,t] = phi_PM[:,t+1]; phi_PM[:,t+1] = np.zeros(nx)
DrugRelease_1D[:,t] = DrugRelease_1D[:,t+1]; DrugRelease_1D[:,t+1] = np.zeros(nx)
CeSwell_1D[t] = CeSwell_1D[t+1]; CeSwell_1D[t+1] = 0
Rb[t] = Rb[t+1]; Rb[t+1] = 0
dX[t] = dX[t+1]; dX[t+1] = 0
TimeT[t+1] = TimeT[t] + dt
TimeR[t+1] = TimeR[t] + dt
CeSwell_1D[t+1] = CeSwell_1D[t]
if (TimeT[t+1]) <= SwellM[0,1]:
rb = SwellM[Si,0] * rb0
CeSwell_1D[t+1] = ( (4/3*np.pi)* ((rb**3) - (rb - rs)**3) )/( (4/3*np.pi)* ((rb0**3) - (rb0-rs)**3) )
if (TimeT[t+1]) >= SwellM[0,cc+1]:
if cc < 10: cc += 1
rb = SwellM[Si,cc] * rb0
CeSwell_1D[t+1] = ( (4/3*np.pi)*((rb**3) - (rb - rs)**3) ) /( (4/3*np.pi)*((rb0**3) - (rb0-rs)**3) )
Rb[t] = rb
dX[t] = dx0 *(rb/rb0)
(Maskfit_1D,invMaskfit_1D,CSol_1D,CH2O_1D,COH_1D,DCSol_1D,DCH2O_1D,phi_SM,phi_PM,phi_WM
) = Solidification(
Maskfit_1D,invMaskfit_1D,CSol_1D,CH2O_1D,COH_1D,DCSol_1D,DCH2O_1D,Sol_conc,t,dt ,rb,dx ,Vpore_1D,phi_SM,phi_PM,phi_WM,Ce0,invtisMask_1D,CeSwell_1D,Cend_1D )
# Solve for H+ concentration from acid dissociation
(Drug_1D,Drugcoo_1D,Col_1D,Colcoo_1D,CH_1D,COH_1D,CH2O_1D,Cpol_1D,Cpolcoo_1D,Cm_1D,Cmcoo_1D
) = AcidDiss(
Drug_1D,Drugcoo_1D,Col_1D,Colcoo_1D,CH_1D,COH_1D,CH2O_1D,Cpol_1D,Cpolcoo_1D,Cm_1D,Cmcoo_1D,pinchDrug_1D,pinchDrugcoo_1D,pinchCol_1D,pinchColcoo_1D,pinchCH_1D,pinchCOH_1D,pinchCH2O_1D,pinchCpol_1D,pinchCpolcoo_1D,pinchCm_1D,pinchCmcoo_1D,t,dt,CSol_1D,Sol_conc ,Cend_1D)
# Solve for Res, Rrs, Ce, Cend, Xc, and Mn
(Ce_1D,Cend_1D,Res_1D,Rrs_1D,Rol_1D,Rm_1D,Col_1D,Cm_1D,Xc_1D,Xext_1D,Xcend_1D,Mn_1D,Mask_1D,invMask_1D,Cpol0,t,DegWeightSciss,omega
) = ScissionRate(
Ce_1D,Cend_1D,Res_1D,Rrs_1D,Rol_1D,Rm_1D,Col_1D,Cm_1D,Xc_1D,Xext_1D, Xcend_1D,Mn_1D,Mask_1D,invMask_1D,Cpol0,Ce0,t,dt ,CH2O_1D,CSol_1D,Sol_conc,DegWeightSciss,rb,dx ,CeSwell_1D,CH_1D )
# Solve for diffusion of constituents
(Drug_1D,Drugcoo_1D,Col_1D,Colcoo_1D,CH_1D,COH_1D,CH2O_1D,Cmcoo_1D,Cm_1D,DDrug_1D,DCol_1D,DCH_1D,DCm_1D,t,A
) = Diffusion(
Drug_1D,Drugcoo_1D,Col_1D,Colcoo_1D,CH_1D,COH_1D,CH2O_1D,Cmcoo_1D,Cm_1D,DDrug_1D,DCol_1D,DCH_1D,DCm_1D,t,dt,A ,rb,dx,Rol_1D,Rm_1D )
# Solve for porosity of polymer
(Maskfit_1D,invMaskfit_1D,VporeCol_1D,VporeDrug_1D,Vpore_1D,Rol_1D,Rm_1D,Col_1D,Colcoo_1D,Cm_1D,Cmcoo_1D,Drug_1D,Drugcoo_1D,DCol_1D,DCm_1D,DDrug_1D,DCH_1D,t,fDrug
) = Porosity(
Maskfit_1D,invMaskfit_1D,VporeCol_1D,VporeDrug_1D,Vpore_1D,Rol_1D,Rm_1D,Col_1D,Colcoo_1D,Cm_1D,Cmcoo_1D,Drug_1D,Drugcoo_1D,DCol_1D,DCm_1D,DDrug_1D,DCH_1D,t,fDrug,Ce0,dt,rb,dx,CeSwell_1D,Mask_1D,Xc_1D,Xcend_1D,Xext_1D )
# Update Diffusivity
(Drug_1D,Drugcoo_1D,Col_1D,Colcoo_1D,CH_1D,COH_1D,CH2O_1D,Cmcoo_1D,Cm_1D,DDrug_1D,DCol_1D,DCH_1D,DCm_1D,t
) = UpdateDiffusion(
Drug_1D,Drugcoo_1D,Col_1D,Colcoo_1D,CH_1D,COH_1D,CH2O_1D,Cmcoo_1D,Cm_1D,DDrug_1D,DCol_1D,DCH_1D,DCm_1D,t,dt ,rb,dx,invtisMask_1D,tisMask_1D,Rol_1D,Rm_1D )
#Outputs
DrugRelease_1D[:,t] =invtisMask_1D[:] * (Drug_1D[:,t] +Drugcoo_1D[:,t])
n+=1
# =============================================================================
# Output Variables
# =============================================================================
# SumDrug = np.sum(DrugRelease_1D_mem*MWdrug*dx ,axis = 0)
# SumMn = np.sum(Mn_1D_mem,axis = 0)
# SumXc = np.sum(Xc_1D_mem,axis = 0)
#
# ProfileMn = np.zeros(len(Mn_1D_mem))
# ProfileDrug = np.zeros(len(Drug_1D_mem))
# ProfileXc = np.zeros(len(Xc_1D_mem))
#
# ProfileMn[:-nstart] = (SumMn[(nstart):])/SumMn[nstart]
# ProfileDrug[:-nstart] =100* (SumDrug[nstart:])/SumDrug[nstart]
# ProfileXc[:-nstart] = (SumXc[nstart:])/SumXc[nstart]
nend = n
return(DDrug_1D_mem,Drug_1D_mem,Drugcoo_1D_mem,Maskfit_1D_mem,CSol_1D_mem,CH2O_1D_mem,DDrug_1D_mem,A,
omega,Ce0,invtisMask_1D,phi_SM_mem,phi_PM_mem,phi_WM_mem,CH_1D_mem,nstart,nend,Ce_1D_mem,Cend_1D_mem,
Xc_1D_mem,Rb,dX,TimeT_mem,TimeR_mem,DrugRelease_1D_mem,CeSwell_1D,
Rrs_1D_mem,Res_1D_mem,Vpore_1D_mem,Mn_1D_mem, invMaskfit_1D_mem, DCH_1D_mem,
Cpol_1D_mem, Cpolcoo_1D_mem, Col_1D_mem, Colcoo_1D_mem, Cm_1D_mem, Cmcoo_1D_mem, Rol_1D_mem,Rm_1D_mem)
Copy_Ce_1D=np.zeros(nx)
Copy_Col_1D=np.zeros(nx)
Copy_Colcoo_1D=np.zeros(nx)
Copy_Cm_1D=np.zeros(nx)
Copy_Cmcoo_1D=np.zeros(nx)
Copy_Cpol_1D=np.zeros(nx)
Copy_Cpolcoo_1D=np.zeros(nx)
Copy_Res_1D=np.zeros(nx)
Copy_Rrs_1D=np.zeros(nx)
Copy_Rol_1D=np.zeros(nx)
Copy_Rm_1D=np.zeros(nx)
Copy_Cend_1D=np.zeros(nx)
Copy_Drug_1D=np.zeros(nx)
Copy_Drugcoo_1D=np.zeros(nx)
Copy_CSol_1D=np.zeros(nx)
Copy_CH2O_1D=np.zeros(nx)
Copy_CH_1D=np.zeros(nx)
Copy_COH_1D=np.zeros(nx)
Copy_CSol_1D = np.zeros(nx)
OCopy_Ce_1D=np.zeros(nx)
OCopy_Col_1D=np.zeros(nx)
OCopy_Colcoo_1D=np.zeros(nx)
OCopy_Cm_1D=np.zeros(nx)
OCopy_Cmcoo_1D=np.zeros(nx)
OCopy_Cpol_1D=np.zeros(nx)
OCopy_Cpolcoo_1D=np.zeros(nx)
OCopy_Res_1D=np.zeros(nx)
OCopy_Rrs_1D=np.zeros(nx)
OCopy_Rol_1D=np.zeros(nx)
OCopy_Rm_1D=np.zeros(nx)
OCopy_Cend_1D=np.zeros(nx)
OCopy_Drug_1D=np.zeros(nx)
OCopy_Drugcoo_1D=np.zeros(nx)
OCopy_CSol_1D=np.zeros(nx)
OCopy_CH2O_1D=np.zeros(nx)
OCopy_CH_1D=np.zeros(nx)
OCopy_COH_1D=np.zeros(nx)
OCopy_CSol_1D = np.zeros(nx)
def Maskfunc(
rb,rp,rs,x,Lx,Mask_1D,invMask_1D,tisMask_1D,invtisMask_1D,dx):
ns = 100;
shifty = thick + 2*rp
shiftx = 2*np.sqrt((shifty**2) - (0.5*shifty)**2)
start = 0
# sy = np.linspace(start,ns-1+start,ns)*shifty
sx = np.linspace(start,ns-1+start,ns)*(shiftx)
I = 0;
while I < ns:
#Creates Interior Pores over entire Domain
tempA = ((x[:]-sx[I])**2) > (rp)**2
Mask_1D[:] = Mask_1D[:] * tempA*1
tempB = ((x[:]-sx[I]-0.5*shiftx)**2) > (rp)**2
Mask_1D[:] = Mask_1D[:] * tempB*1
#Removes everything except core region, where interior pores exist
tempc =((x[:]-0*0.5*Lx)**2) < (rb-rs)**2
Mask_1D[:] = Mask_1D[:] * tempc*1
#Adds everything outside of the interior core, sets up for making solid shell
tempc2 =((x[:]-0*0.5*Lx)**2) > (rb-rs)**2
Mask_1D[:] = Mask_1D[:] + tempc2*1
#Crops shell
tempd =((x[:]-0*0.5*Lx)**2) < (rb)**2
Mask_1D[:] = Mask_1D[:] * tempd*1
# =============================================================================
#Shell pore Circle base
tempE = ((x[:]-0*0.5*Lx-((rb-rs)+0.5*rs))**2) > (rps)**2
Mask_1D[:] = Mask_1D[:] * tempE*1
# tempE = ((x[:]-0*0.5*Lx-0*((rb-rs)+0.5*rs))**2) > (4*rp)**2
# Mask_1D[:] = Mask_1D[:] * tempE*1
# tempE = ((x[:]-0*0.5*Lx+((rb-rs)+0.5*rs))**2) > (4*rp)**2
# Mask_1D[:] = Mask_1D[:] * tempE*1
# tempE = ((x[:]-0*0.5*Lx-0*((rb-rs)+0.5*rs))**2) > (4*rp)**2
# Mask_1D[:] = Mask_1D[:] * tempE*1
I+=1
M = (Mask_1D[:] < 1)
invMask_1D = M*1
tis = ((x[:]-0*0.5*Lx)**2) > (rb)**2
tisMask_1D[:] = tis*1
invtisMask_1D = 1-tisMask_1D
return(Mask_1D,invMask_1D,tisMask_1D,invtisMask_1D)
def Initialize( Mask_1D,invMask_1D,Maskfit_1D,invMaskfit_1D,invtisMask_1D,tisMask_1D,CSol_1D,CH2O_1D,rb,dx, phi_SM,phi_PM,phi_WM):
Vol = 5e-8#5.01507e-8
Mass = Vol/(0.39/rho_p + 0.6/rho_s + 0.01/rho_d)
Polymass = 0.39 * Mass
Solmass = 0.6 * Mass
DrugTotal = 0.01 * Mass
DrugTotalmol = DrugTotal/Mw_d
# PoreV = (4/3*np.pi*(rb-rs)**3)*0.74048 * (rp)**3/(rp+0.5*thick)**3
# PolymerV = (4/3*np.pi*(rb-rs)**3) - PoreV + 4/3*np.pi*(rb**3 - (rb-rs)**3)
Drug0 = (DrugTotalmol)/(Vol)
Drug_1D[:,0] = Drug0 * invtisMask_1D[:]
Drug_1D[:,1] = Drug0 * invtisMask_1D[:]
fDrug = 0
##############################################################################
Ce0 = (1-Xc0)* 1000*(Polymass/MWpolyunit)/(Vol) #mol/m^3
omega = Ce0
Cpol0 = Ce0/mm
Mn_1D[:,0] = Mask_1D[:] * ((Ce0+omega*Xc0) * M0)/(Cpol0 + invMask_1D[:])
Mn_1D[:,1] = Mask_1D[:] * ((Ce0+omega*Xc0) * M0)/(Cpol0 + invMask_1D[:])
Sol_conc = 1000*(Solmass/99.133)/(Vol)
phi_S = invtisMask_1D*((Solmass/1/99.133)/(Vol) *Mw_s/rho_s)
phi_P =invtisMask_1D* (MW*(Polymass/MW/Vol)/rho_p)
phi_W = (1-phi_P-phi_S)
CH2O0 = phi_W * 1000*rho_w/Mw_w
CSol_1D[:,0] = 1000* phi_S*rho_s/Mw_s
CSol_1D[:,1] = 1000* phi_S*rho_s/Mw_s
phi_SM[:,1] = phi_S; phi_SM[:,0] = phi_S
phi_PM[:,1] = phi_P; phi_PM[:,0] = phi_P
phi_WM[:,1] = phi_W; phi_WM[:,0] = phi_W
##############################################################################
Ce_1D[:,0] = Ce0*Mask_1D[:]
Ce_1D[:,1] = Ce0*Mask_1D[:]
Cend_1D[:,0] = Mask_1D[:]* 2*Cpol0
Cend_1D[:,1] = Mask_1D[:]* 2*Cpol0
Cpol_1D[:,0] = Cpol0*Mask_1D[:]
Cpol_1D[:,1] = Cpol0*Mask_1D[:]
Col_1D[:,0] = Col0*Mask_1D[:]
Col_1D[:,1] = Col0*Mask_1D[:]
Cm_1D[:,0] = Cm0*Mask_1D[:]
Cm_1D[:,1] = Cm0*Mask_1D[:]
Rs_1D[:,0] = 0*Mask_1D[:]
CH_1D[:,0] = CH0 #0.1*(10**-7.4)*Mask[0:,0:] + (10**-7.4)*invMask[0:,0:];
CH2O_1D[:,0] = CH2O0
CH2O_1D[:,1] = CH2O0
Maskfit_1D[:,0] = 0 * Mask_1D[:]
Maskfit_1D[:,1] = 0 * Mask_1D[:]
invMaskfit_1D[:,0] = 1 - Maskfit_1D[:,0]
invMaskfit_1D[:,1] = 1 - Maskfit_1D[:,1]
DCol_1D[:,0] = invMaskfit_1D[:,0]*Dpore_Col + Maskfit_1D[:,0]*(Dpolymer_Col)
DCm_1D[:,0] = invMaskfit_1D[:,0]*Dpore_Cm + Maskfit_1D[:,0]*(Dpolymer_Cm)
DDrug_1D[:,0] = invMaskfit_1D[:,0]*Dpore_Drug + Maskfit_1D[:,0]*(Dpolymer_Drug)
DCH_1D[:,0] = invMaskfit_1D[:,0]*Dpore_CH + Maskfit_1D[:,0]*(Dpolymer_CH)
DCol_1D[:,1] = invMaskfit_1D[:,1]*Dpore_Col + Maskfit_1D[:,1]*(Dpolymer_Col)
DCm_1D[:,1] = invMaskfit_1D[:,1]*Dpore_Cm + Maskfit_1D[:,1]*(Dpolymer_Cm)
DDrug_1D[:,1] = invMaskfit_1D[:,1]*Dpore_Drug + Maskfit_1D[:,1]*(Dpolymer_Drug)
DCH_1D[:,1] = invMaskfit_1D[:,1]*Dpore_CH + Maskfit_1D[:,1]*(Dpolymer_CH)
DCH2O_1D[0:,0] = invMaskfit_1D[:,0]*Dpore_H2O + Maskfit_1D[:,0]*(Dpolymer_H2O)
DCH2O_1D[0:,1] = invMaskfit_1D[:,1]*Dpore_H2O + Maskfit_1D[:,1]*(Dpolymer_H2O)
DCSol_1D[0:,0] = invMaskfit_1D[:,0]*Dpore_Sol + Maskfit_1D[:,0]*(Dpolymer_Sol)
DCSol_1D[0:,1] = invMaskfit_1D[:,1]*Dpore_Sol + Maskfit_1D[:,1]*(Dpolymer_Sol)
Drug_1D[-1,0] = 0
Col_1D[-1,0] = 0
Drugcoo_1D[-1,0] = 0
Colcoo_1D[-1,0] = 0
CH_1D[-1,0] = CH0
COH_1D[-1,0] = 0
CSol_1D[-1,0] = 0
Cpol_1D[-1,0] = 0
Cpolcoo_1D[-1,0] = 0
Drug_1D[-1,1] = 0
Col_1D[-1,1] = 0
Drugcoo_1D[-1,1] = 0
Colcoo_1D[-1,1] = 0
CH_1D[-1,1] = 0
COH_1D[-1,1] = 0
CSol_1D[-1,1] = 0
Cpol_1D[-1,1] = 0
Cpolcoo_1D[-1,1] = 0
Drug_1D[0,0] = Drug_1D[1,0]
Col_1D[0,0] = 0
Drugcoo_1D[0,0] = 0
Colcoo_1D[0,0] = 0
CH_1D[0,0] = 0
COH_1D[0,0] = 0
CSol_1D[0,0] = 0
Cpol_1D[0,0] = 0
Cpolcoo_1D[0,0] = 0
Drug_1D[0,1] = Drug_1D[1,0]
Col_1D[0,1] = 0
Drugcoo_1D[0,1] = 0
Colcoo_1D[0,1] = 0
CH_1D[0,1] = CH0
COH_1D[0,1] = 0
CSol_1D[0,1] = 0
Cpol_1D[0,1] = 0
Cpolcoo_1D[0,1] = 0
return(Drug_1D,Cpol_1D,Col_1D,Cm_1D,CH_1D,DDrug_1D,DCol_1D,DCm_1D,fDrug,Maskfit_1D,invMaskfit_1D,Cpol0,Ce0,CSol_1D,CH2O_1D,DCSol_1D,DCH2O_1D,Sol_conc,Ce0, phi_SM,phi_PM,phi_WM)
def Solidification(Maskfit_1D,invMaskfit_1D,CSol_1D,CH2O_1D,COH_1D,DCSol_1D,DCH2O_1D,Sol_conc,t,dt ,rb,dx,Vpore_1D,phi_SM,phi_PM,phi_WM ,Ce0,invtisMask_1D,CeSwell_1D,Cend_1D):
r = np.linspace(0,Lx,nx)
CH2O_1D[0,t] = CH2O_1D[1,t]
CH2O_1D[1:,t] = invtisMask_1D[1:]*CH2O_1D[1:,t] + (1-invtisMask_1D[1:]) *CH2O_1D[1:,0]
CSol_1D[0,t] = CSol_1D[1,t]
phi_SM[:,t] = invtisMask_1D*(Mw_s*CSol_1D[:,t] /1000/rho_s)
phi_PM[:,t] = invtisMask_1D* (MW*(Polymass/MW/(Vol*CeSwell_1D[t+1]))/rho_p)#(( Polymass/1/MWpolyunit)/(Vol*CeSwell_1D[t+1]) - Col0)/mm*MW/rho_p
phi_WM[:,t] = invtisMask_1D[:]*(1 - phi_SM[:,t] - phi_PM[:,t])
phi_WM[phi_WM[:,t]<=0] = 0
RT = 8.314*300
a1 = Mw_w*RT/(Dpore_H2O*rho_w)
A1 = Mw_w*RT/(Dpolymer_H2O*rho_w)
b1 = rho_s*Mw_w/(rho_w*Mw_s)
a2 = Mw_s*RT/(Dpore_Sol*rho_s)
A2 = Mw_s*RT/(Dpolymer_Sol*rho_s)
b2 = rho_w*Mw_s/(rho_s*Mw_w)
from scipy.optimize import fsolve # Move this outside of the loop move variables into definitions
def equations(p):
x, y = p
with np.errstate(invalid='ignore'):
return (a1-b1*np.sqrt(x*y)-x, a2-b2*np.sqrt(x*y)-y)
def equations2(p):
x, y = p
with np.errstate(invalid='ignore'):
return (A1-b1*np.sqrt(x*y)-x, A2-b2*np.sqrt(x*y)-y)
D11Bar_theo, D22Bar_theo = fsolve(equations, (10000, 10000))
D11Bar_poly_theo, D22Bar_poly_theo = fsolve(equations2, (10000, 10000))
D12Bar_theo= np.sqrt(D11Bar_theo*D22Bar_theo)
D12Bar_poly_theo= np.sqrt(D11Bar_poly_theo*D22Bar_poly_theo)
D12Bar_poly_Swell = (1/D12Bar_poly_theo) * np.exp( -Beta*(1- (CeSwell_1D[t+1]))) #*(2e4)
D12bar = invMaskfit_1D[:,t]*1/D12Bar_theo + Maskfit_1D[:,t]*(D12Bar_poly_Swell+(1.3*Vpore_1D[:,t-1]**2-0.3*Vpore_1D[:,t-1]**3)*(1/D12Bar_theo-D12Bar_poly_Swell))
w1 = 1#2/3
w2 = 1#4/3
w3 = 0#1/3
r = np.linspace(0,Lx,nx)#/Lx
dr = r[1]-r[0]
I = np.linspace(1,nx,nx)
# Water Diffusion
sig_xRC = np.zeros(nx)
sig_xLC = np.zeros(nx)
sig_xyC = np.zeros(nx)
sig_xRC[2:nx] = -(w1)*(dt/(dr)**2)*(D12bar[1:nx-1]/2/I[1:nx-1]*I[2:nx]+D12bar[2:nx]/2 )
sig_xLC[0:nx-2] = -(w1)*(dt/(dr)**2)*(D12bar[1:nx-1]/2/I[1:nx-1]*I[:nx-2]+D12bar[0:nx-2]/2 )
sig_xyC[1:nx-1] = 1 + (w1)*(dt/(dr)**2)*(D12bar[1:nx-1] + D12bar[2:nx]/2+ D12bar[:nx-2]/2)
sig_xyC[0] = sig_xyC[-1] = 1
A5 = sparse.spdiags([sig_xyC[:],sig_xRC[:],sig_xLC[:]],[0,1,-1],nx,nx)
B5 = sparse.spdiags([(w2)*np.ones(nx)],[0],nx,nx)
C5 = sparse.spdiags([(w3)*np.ones(nx)],[0],nx,nx)
# Solvent Diffusion
sig_xRC = np.zeros(nx)
sig_xLC = np.zeros(nx)
sig_xyC = np.zeros(nx)
sig_xRC[2:nx] = -(w1)*(dt/(dr)**2)*(D12bar[1:nx-1]/2/I[1:nx-1]*I[2:nx]+D12bar[2:nx]/2 )
sig_xLC[0:nx-2] = -(w1)*(dt/(dr)**2)*(D12bar[1:nx-1]/2/I[1:nx-1]*I[:nx-2]+D12bar[0:nx-2]/2 )
sig_xyC[1:nx-1] = 1 + (w1)*(dt/(dr)**2)*(D12bar[1:nx-1] + D12bar[2:nx]/2+ D12bar[:nx-2]/2)
sig_xyC[0] = sig_xyC[-1] = 1
A6 = sparse.spdiags([sig_xyC[:],sig_xRC[:],sig_xLC[:]],[0,1,-1],nx,nx)
B6 = sparse.spdiags([(w2)*np.ones(nx)],[0],nx,nx)
C6 = sparse.spdiags([(w3)*np.ones(nx)],[0],nx,nx)
bh2o = (B5.dot(Mw_w/rho_w/1000 * CH2O_1D[:,t] + phi_PM[:,t]) )
ch2o = (C5.dot(Mw_w/rho_w/1000 * CH2O_1D[:,t-1] + phi_PM[:,t-1]) )
Bh2o = bh2o-ch2o
bsol = (B6.dot(Mw_s/rho_s/1000 * CSol_1D[:,t] ) )
csol = (C6.dot(Mw_s/rho_s/1000 * CSol_1D[:,t-1] ) )
Bsol = bsol-csol
temph = rho_w/Mw_w*1000* (spsolve(A5, Bh2o)- phi_PM[:,t])
temps = rho_s/Mw_s*1000* spsolve(A6, Bsol)
CH2O_1D[:,t+1] = temph[:]
CSol_1D[:,t+1] = temps[:]
CH2O_1D[-1,t+1] = CH2O_1D[-1,0]
# =============================================================================
# Solidification
# =============================================================================
Maskfit_1D[:,t+1] = TFunc(phi_WM[:,t]) * Mask_1D[:]
# if MW == 15000: Maskfit_1D[:,t+1] = TFunc15(phi_WM[:,t]) * Mask_1D[:]
# if MW == 29000: Maskfit_1D[:,t+1] = TFunc29(phi_WM[:,t]) * Mask_1D[:]
# if MW == 53000: Maskfit_1D[:,t+1] = TFunc53(phi_WM[:,t]) * Mask_1D[:]
if MW <= 33000: Maskfit_1D[:,t+1] = TFunc8_33(phi_WM[:,t],MW) * Mask_1D[:]
elif MW >= 50000: Maskfit_1D[:,t+1] = TFunc50_60(phi_WM[:,t],MW) * Mask_1D[:]
else: Maskfit_1D[:,t+1] = TFunc33_50(phi_WM[:,t],MW) * Mask_1D[:]
# Maskfit_1D[(Maskfit_1D[:,t+1])>=1,t+1] = 1
# Maskfit_1D[(Maskfit_1D[:,t+1])>=0.9999,t+1] = 1
# Maskfit_1D[0,t+1] = Maskfit_1D[1,t+1]
invMaskfit_1D[:,t+1] = 1 - Maskfit_1D[:,t+1]
# =============================================================================
# Update Diffusion
# =============================================================================
D12Bar_poly_Swell = (1/D12Bar_poly_theo) * np.exp( -Beta*(1- (CeSwell_1D[t+1])))
D12bar = invMaskfit_1D[:,t+1]*1/D12Bar_theo + Maskfit_1D[:,t+1]*(D12Bar_poly_Swell + (1.3*Vpore_1D[:,t]**2-0.3*Vpore_1D[:,t]**3)*(1/D12Bar_theo - D12Bar_poly_Swell))#Dpolymer_Sol_Swell))
w1 = 1#2/3
w2 = 1#4/3
w3 = 0#1/3
r = np.linspace(0,Lx,nx)#/Lx
dr = r[1]-r[0]
I = np.linspace(1,nx,nx)
# Water Diffusion
sig_xRC = np.zeros(nx)
sig_xLC = np.zeros(nx)
sig_xyC = np.zeros(nx)
sig_xRC[2:nx] = -(w1)*(dt/(dr)**2)*(D12bar[1:nx-1]/2/I[1:nx-1]*I[2:nx]+D12bar[2:nx]/2 )
sig_xLC[0:nx-2] = -(w1)*(dt/(dr)**2)*(D12bar[1:nx-1]/2/I[1:nx-1]*I[:nx-2]+D12bar[0:nx-2]/2 )
sig_xyC[1:nx-1] = 1 + (w1)*(dt/(dr)**2)*(D12bar[1:nx-1] + D12bar[2:nx]/2+ D12bar[:nx-2]/2)
sig_xyC[0] = sig_xyC[-1] = 1
A5 = sparse.spdiags([sig_xyC[:],sig_xRC[:],sig_xLC[:]],[0,1,-1],nx,nx)
B5 = sparse.spdiags([(w2)*np.ones(nx)],[0],nx,nx)
C5 = sparse.spdiags([(w3)*np.ones(nx)],[0],nx,nx)
# Solvent Diffusion
sig_xRC = np.zeros(nx)
sig_xLC = np.zeros(nx)
sig_xyC = np.zeros(nx)
sig_xRC[2:nx] = -(w1)*(dt/(dr)**2)*(D12bar[1:nx-1]/2/I[1:nx-1]*I[2:nx]+D12bar[2:nx]/2 )
sig_xLC[0:nx-2] = -(w1)*(dt/(dr)**2)*(D12bar[1:nx-1]/2/I[1:nx-1]*I[:nx-2]+D12bar[0:nx-2]/2 )
sig_xyC[1:nx-1] = 1 + (w1)*(dt/(dr)**2)*(D12bar[1:nx-1] + D12bar[2:nx]/2+ D12bar[:nx-2]/2)
sig_xyC[0] = sig_xyC[-1] = 1
A6 = sparse.spdiags([sig_xyC[:],sig_xRC[:],sig_xLC[:]],[0,1,-1],nx,nx)
B6 = sparse.spdiags([(w2)*np.ones(nx)],[0],nx,nx)
C6 = sparse.spdiags([(w3)*np.ones(nx)],[0],nx,nx)
bh2o = (B5.dot(Mw_w/rho_w/1000 * CH2O_1D[:,t] + phi_PM[:,t]) )
ch2o = (C5.dot(Mw_w/rho_w/1000 * CH2O_1D[:,t-1] + phi_PM[:,t-1]) )
Bh2o = bh2o-ch2o
bsol = (B6.dot(Mw_s/rho_s/1000 * CSol_1D[:,t] ) )
csol = (C6.dot(Mw_s/rho_s/1000 * CSol_1D[:,t-1] ) )
Bsol = bsol-csol
temph = rho_w/Mw_w*1000* (spsolve(A5, Bh2o)- phi_PM[:,t])
temps = rho_s/Mw_s*1000* spsolve(A6, Bsol)
CH2O_1D[:,t+1] = temph[:]
CSol_1D[:,t+1] = temps[:]
CH2O_1D[-1,t+1] = CH2O_1D[-1,0]
COH_1D[:,t+1] = COH_1D[:,t]
return Maskfit_1D,invMaskfit_1D,CSol_1D,CH2O_1D,COH_1D,DCSol_1D,DCH2O_1D,phi_SM,phi_PM,phi_WM
def RK(fCol,fColcoo,fCm,fCmcoo,fCpol,fCpolcoo,fDrug,fDrugcoo,fCH,Col0,Colcoo0,Cm0,Cmcoo0,Cpol0,Cpolcoo0,Drug0,Drugcoo0,CH0,dt,ntt,t):
col = np.zeros((nx,ntt + 1))
colcoo = np.zeros((nx,ntt + 1))
cm = np.zeros((nx,ntt + 1))
cmcoo = np.zeros((nx,ntt + 1))
cpol = np.zeros((nx,ntt + 1))
cpolcoo = np.zeros((nx,ntt + 1))
drug = np.zeros((nx,ntt + 1))
drugcoo = np.zeros((nx,ntt + 1))
ch = np.zeros((nx,ntt + 1))
ut = np.zeros(ntt + 1)
h = (dt) / float(ntt)
col[:,0] = Col0
colcoo[:,0] = Colcoo0
cm[:,0] = Cm0
cmcoo[:,0] = Cmcoo0
cpol[:,0] = Cpol0
cpolcoo[:,0] = Cpolcoo0
drug[:,0] = Drug0
drugcoo[:,0] = Drugcoo0
ch[:,0] = CH0
ut[0] = TimeT[t]
for i in range(0, ntt):
Col = col[:,i]
Colcoo = colcoo[:,i]
Cm = cm[:,i]
Cmcoo = cmcoo[:,i]
Cpol = cpol[:,i]
Cpolcoo = cpolcoo[:,i]
Drug = drug[:,i]
Drugcoo = drugcoo[:,i]
CH = ch[:,i]
Ut = ut[i]
k1 = h* fCol(Ut,Col,Colcoo,Cm,Cmcoo,Cpol,Cpolcoo,Drug,Drugcoo,CH)
l1 = h* fColcoo(Ut,Col,Colcoo,Cm,Cmcoo,Cpol,Cpolcoo,Drug,Drugcoo,CH)
m1 = h* fCm(Ut,Col,Colcoo,Cm,Cmcoo,Cpol,Cpolcoo,Drug,Drugcoo,CH)
n1 = h* fCmcoo(Ut,Col,Colcoo,Cm,Cmcoo,Cpol,Cpolcoo,Drug,Drugcoo,CH)
o1 = h* fCpol(Ut,Col,Colcoo,Cm,Cmcoo,Cpol,Cpolcoo,Drug,Drugcoo,CH)
p1 = h* fCpolcoo(Ut,Col,Colcoo,Cm,Cmcoo,Cpol,Cpolcoo,Drug,Drugcoo,CH)
q1 = h* fDrug(Ut,Col,Colcoo,Cm,Cmcoo,Cpol,Cpolcoo,Drug,Drugcoo,CH)
r1 = h* fDrugcoo(Ut,Col,Colcoo,Cm,Cmcoo,Cpol,Cpolcoo,Drug,Drugcoo,CH)
s1 = h* fCH(Ut,Col,Colcoo,Cm,Cmcoo,Cpol,Cpolcoo,Drug,Drugcoo,CH)
k2 = h* fCol(Ut+0.5*h,Col+0.5*k1,Colcoo+0.5*l1,Cm+0.5*m1,Cmcoo+0.5*n1,Cpol+0.5*o1,Cpolcoo+0.5*p1,Drug+0.5*q1,Drugcoo+0.5*r1,CH+0.5*s1)
l2 = h* fColcoo(Ut+0.5*h,Col+0.5*k1,Colcoo+0.5*l1,Cm+0.5*m1,Cmcoo+0.5*n1,Cpol+0.5*o1,Cpolcoo+0.5*p1,Drug+0.5*q1,Drugcoo+0.5*r1,CH+0.5*s1)
m2 = h* fCm(Ut+0.5*h,Col+0.5*k1,Colcoo+0.5*l1,Cm+0.5*m1,Cmcoo+0.5*n1,Cpol+0.5*o1,Cpolcoo+0.5*p1,Drug+0.5*q1,Drugcoo+0.5*r1,CH+0.5*s1)
n2 = h* fCmcoo(Ut+0.5*h,Col+0.5*k1,Colcoo+0.5*l1,Cm+0.5*m1,Cmcoo+0.5*n1,Cpol+0.5*o1,Cpolcoo+0.5*p1,Drug+0.5*q1,Drugcoo+0.5*r1,CH+0.5*s1)
o2 = h* fCpol(Ut+0.5*h,Col+0.5*k1,Colcoo+0.5*l1,Cm+0.5*m1,Cmcoo+0.5*n1,Cpol+0.5*o1,Cpolcoo+0.5*p1,Drug+0.5*q1,Drugcoo+0.5*r1,CH+0.5*s1)
p2 = h* fCpolcoo(Ut+0.5*h,Col+0.5*k1,Colcoo+0.5*l1,Cm+0.5*m1,Cmcoo+0.5*n1,Cpol+0.5*o1,Cpolcoo+0.5*p1,Drug+0.5*q1,Drugcoo+0.5*r1,CH+0.5*s1)
q2 = h* fDrug(Ut+0.5*h,Col+0.5*k1,Colcoo+0.5*l1,Cm+0.5*m1,Cmcoo+0.5*n1,Cpol+0.5*o1,Cpolcoo+0.5*p1,Drug+0.5*q1,Drugcoo+0.5*r1,CH+0.5*s1)
r2 = h* fDrugcoo(Ut+0.5*h,Col+0.5*k1,Colcoo+0.5*l1,Cm+0.5*m1,Cmcoo+0.5*n1,Cpol+0.5*o1,Cpolcoo+0.5*p1,Drug+0.5*q1,Drugcoo+0.5*r1,CH+0.5*s1)
s2 = h* fCH(Ut+0.5*h,Col+0.5*k1,Colcoo+0.5*l1,Cm+0.5*m1,Cmcoo+0.5*n1,Cpol+0.5*o1,Cpolcoo+0.5*p1,Drug+0.5*q1,Drugcoo+0.5*r1,CH+0.5*s1)
k3 = h* fCol(Ut+0.5*h,Col+0.5*k2,Colcoo+0.5*l2,Cm+0.5*m2,Cmcoo+0.5*n2,Cpol+0.5*o2,Cpolcoo+0.5*p2,Drug+0.5*q2,Drugcoo+0.5*r2,CH+0.5*s2)
l3 = h* fColcoo(Ut+0.5*h,Col+0.5*k2,Colcoo+0.5*l2,Cm+0.5*m2,Cmcoo+0.5*n2,Cpol+0.5*o2,Cpolcoo+0.5*p2,Drug+0.5*q2,Drugcoo+0.5*r2,CH+0.5*s2)
m3 = h* fCm(Ut+0.5*h,Col+0.5*k2,Colcoo+0.5*l2,Cm+0.5*m2,Cmcoo+0.5*n2,Cpol+0.5*o2,Cpolcoo+0.5*p2,Drug+0.5*q2,Drugcoo+0.5*r2,CH+0.5*s2)
n3 = h* fCmcoo(Ut+0.5*h,Col+0.5*k2,Colcoo+0.5*l2,Cm+0.5*m2,Cmcoo+0.5*n2,Cpol+0.5*o2,Cpolcoo+0.5*p2,Drug+0.5*q2,Drugcoo+0.5*r2,CH+0.5*s2)
o3 = h* fCpol(Ut+0.5*h,Col+0.5*k2,Colcoo+0.5*l2,Cm+0.5*m2,Cmcoo+0.5*n2,Cpol+0.5*o2,Cpolcoo+0.5*p2,Drug+0.5*q2,Drugcoo+0.5*r2,CH+0.5*s2)
p3 = h* fCpolcoo(Ut+0.5*h,Col+0.5*k2,Colcoo+0.5*l2,Cm+0.5*m2,Cmcoo+0.5*n2,Cpol+0.5*o2,Cpolcoo+0.5*p2,Drug+0.5*q2,Drugcoo+0.5*r2,CH+0.5*s2)
q3 = h* fDrug(Ut+0.5*h,Col+0.5*k2,Colcoo+0.5*l2,Cm+0.5*m2,Cmcoo+0.5*n2,Cpol+0.5*o2,Cpolcoo+0.5*p2,Drug+0.5*q2,Drugcoo+0.5*r2,CH+0.5*s2)
r3 = h* fDrugcoo(Ut+0.5*h,Col+0.5*k2,Colcoo+0.5*l2,Cm+0.5*m2,Cmcoo+0.5*n2,Cpol+0.5*o2,Cpolcoo+0.5*p2,Drug+0.5*q2,Drugcoo+0.5*r2,CH+0.5*s2)
s3 = h* fCH(Ut+0.5*h,Col+0.5*k2,Colcoo+0.5*l2,Cm+0.5*m2,Cmcoo+0.5*n2,Cpol+0.5*o2,Cpolcoo+0.5*p2,Drug+0.5*q2,Drugcoo+0.5*r2,CH+0.5*s2)
k4 = h* fCol(Ut+ h,Col+ k3,Colcoo+ l3,Cm+ m3,Cmcoo+ n3,Cpol+ o3,Cpolcoo+ p3,Drug+ q3,Drugcoo+ r3,CH+ s3)
l4 = h* fColcoo(Ut+ h,Col+ k3,Colcoo+ l3,Cm+ m3,Cmcoo+ n3,Cpol+ o3,Cpolcoo+ p3,Drug+ q3,Drugcoo+ r3,CH+ s3)
m4 = h* fCm(Ut+ h,Col+ k3,Colcoo+ l3,Cm+ m3,Cmcoo+ n3,Cpol+ o3,Cpolcoo+ p3,Drug+ q3,Drugcoo+ r3,CH+ s3)
n4 = h* fCmcoo(Ut+ h,Col+ k3,Colcoo+ l3,Cm+ m3,Cmcoo+ n3,Cpol+ o3,Cpolcoo+ p3,Drug+ q3,Drugcoo+ r3,CH+ s3)
o4 = h* fCpol(Ut+ h,Col+ k3,Colcoo+ l3,Cm+ m3,Cmcoo+ n3,Cpol+ o3,Cpolcoo+ p3,Drug+ q3,Drugcoo+ r3,CH+ s3)
p4 = h* fCpolcoo(Ut+ h,Col+ k3,Colcoo+ l3,Cm+ m3,Cmcoo+ n3,Cpol+ o3,Cpolcoo+ p3,Drug+ q3,Drugcoo+ r3,CH+ s3)
q4 = h* fDrug(Ut+ h,Col+ k3,Colcoo+ l3,Cm+ m3,Cmcoo+ n3,Cpol+ o3,Cpolcoo+ p3,Drug+ q3,Drugcoo+ r3,CH+ s3)
r4 = h* fDrugcoo(Ut+ h,Col+ k3,Colcoo+ l3,Cm+ m3,Cmcoo+ n3,Cpol+ o3,Cpolcoo+ p3,Drug+ q3,Drugcoo+ r3,CH+ s3)
s4 = h* fCH(Ut+ h,Col+ k3,Colcoo+ l3,Cm+ m3,Cmcoo+ n3,Cpol+ o3,Cpolcoo+ p3,Drug+ q3,Drugcoo+ r3,CH+ s3)
ut[i+1] = TimeT[t] + (i+1)*h
col[:,i+1] = col[:,i]+ (k1 + 2*k2 + 2*k3 + k4) / 6
colcoo[:,i+1] = colcoo[:,i]+ (l1 + 2*l2 + 2*l3 + l4) / 6
cm[:,i+1] = cm[:,i]+ (m1 + 2*m2 + 2*m3 + m4) / 6
cmcoo[:,i+1] = cmcoo[:,i]+ (n1 + 2*n2 + 2*n3 + n4) / 6
cpol[:,i+1] = cpol[:,i]+ (o1 + 2*o2 + 2*o3 + o4) / 6
cpolcoo[:,i+1] = cpolcoo[:,i]+ (p1 + 2*p2 + 2*p3 + p4) / 6
drug[:,i+1] = drug[:,i]+ (q1 + 2*q2 + 2*q3 + q4) / 6
drugcoo[:,i+1] = drugcoo[:,i]+ (r1 + 2*r2 + 2*r3 + r4) / 6
ch[:,i+1] = ch[:,i]+ (s1 + 2*s2 + 2*s3 + s4) / 6
return ut,col,colcoo,cm,cmcoo,cpol,cpolcoo,drug,drugcoo,ch
def fCol(Ut,Col,Colcoo,Cm,Cmcoo,Cpol,Cpolcoo,Drug,Drugcoo,CH):
return (konol*Colcoo/m*CH - koffol*Col/m)
def fColcoo(Ut,Col,Colcoo,Cm,Cmcoo,Cpol,Cpolcoo,Drug,Drugcoo,CH):
return (koffol*Col/m - konol*Colcoo/m*CH)
def fCm(Ut,Col,Colcoo,Cm,Cmcoo,Cpol,Cpolcoo,Drug,Drugcoo,CH):
return (konm*Cmcoo*CH - koffm*Cm)
def fCmcoo(Ut,Col,Colcoo,Cm,Cmcoo,Cpol,Cpolcoo,Drug,Drugcoo,CH):
return (koffm*Cm - konm*Cmcoo*CH)
def fCpol(Ut,Col,Colcoo,Cm,Cmcoo,Cpol,Cpolcoo,Drug,Drugcoo,CH):
return (konpol*Cpolcoo*CH - koffpol*Cpol)
def fCpolcoo(Ut,Col,Colcoo,Cm,Cmcoo,Cpol,Cpolcoo,Drug,Drugcoo,CH):
return (koffpol*Cpol - konpol*Cpolcoo*CH)
def fDrug(Ut,Col,Colcoo,Cm,Cmcoo,Cpol,Cpolcoo,Drug,Drugcoo,CH):
return (kondrug*Drugcoo*CH - koffdrug*Drug)
def fDrugcoo(Ut,Col,Colcoo,Cm,Cmcoo,Cpol,Cpolcoo,Drug,Drugcoo,CH):
return (koffdrug*Drug - kondrug*Drugcoo*CH)
def fCH(Ut,Col,Colcoo,Cm,Cmcoo,Cpol,Cpolcoo,Drug,Drugcoo,CH):
return (
koffol*Col/m - konol*Colcoo/m*CH
+ koffm*Cm - konm*Cmcoo*CH
+ koffpol*Cpol - konpol*Cpolcoo*CH
+ koffdrug*Drug - kondrug*Drugcoo*CH
)
def ODEModel(Molecules,tt):
(Col,Colcoo,Cm,Cmcoo,Cpol,Cpolcoo,Drug,Drugcoo,CH)=Molecules
dColdt = (konol*Colcoo/m*CH - koffol*Col/m)
dColcoodt = (koffol*Col/m - konol*Colcoo/m*CH)
dCmdt = (konm*Cmcoo*CH - koffm*Cm)
dCmcoodt = (koffm*Cm - konm*Cmcoo*CH)
dCpoldt = (konpol*Cpolcoo*CH - koffpol*Cpol)
dCpolcoodt = (koffpol*Cpol - konpol*Cpolcoo*CH)
dDrugdt = (kondrug*Drugcoo*CH - koffdrug*Drug)
dDrugcoodt = (koffdrug*Drug - kondrug*Drugcoo*CH)
dCHdt = (koffol*Col/m - konol*Colcoo/m*CH
+ koffm*Cm - konm*Cmcoo*CH
+ koffpol*Cpol - konpol*Cpolcoo*CH
+ koffdrug*Drug - kondrug*Drugcoo*CH)
return dColdt,dColcoodt,dCmdt,dCmcoodt,dCpoldt,dCpolcoodt,dDrugdt,dDrugcoodt,dCHdt
def AcidDiss(
Drug_1D,Drugcoo_1D,Col_1D,Colcoo_1D,CH_1D,COH_1D,CH2O_1D,Cpol_1D,Cpolcoo_1D,Cm_1D,Cmcoo_1D,
pinchDrug_1D,pinchDrugcoo_1D,pinchCol_1D,pinchColcoo_1D,pinchCH_1D,pinchCOH_1D,pinchCH2O_1D,
pinchCpol_1D,pinchCpolcoo_1D,pinchCm_1D,pinchCmcoo_1D,t,dt,CSol_1D,Sol_conc,Cend_1D ):
pinchCol0 = Col_1D[:,t]
pinchColcoo0 = Colcoo_1D[:,t]
pinchCm0 = Cm_1D[:,t]
pinchCmcoo0 = Cmcoo_1D[:,t]
pinchCpol0 = 0.5* Cend_1D[:,t]
pinchCpolcoo0 = Cpolcoo_1D[:,t]
pinchDrug0 = Drug_1D[:,t]
pinchDrugcoo0 = Drugcoo_1D[:,t]
pinchCH0 = CH_1D[:,t]
tt = np.linspace(0,dt,ntt)
init_states = np.array([pinchCol0,pinchColcoo0,pinchCm0,pinchCmcoo0,pinchCpol0,pinchCpolcoo0,pinchDrug0,pinchDrugcoo0,pinchCH0])
for i in range(0,nx):
init_state = init_states[:,i]
Molecules = odeint(ODEModel,init_state,tt)
pinchCol_1D[i,:] = Molecules[:,0]
pinchColcoo_1D[i,:] = Molecules[:,1]
pinchCm_1D[i,:] = Molecules[:,2]
pinchCmcoo_1D[i,:] = Molecules[:,3]
pinchCpol_1D[i,:] = Molecules[:,4]
pinchCpolcoo_1D[i,:] = Molecules[:,5]
pinchDrug_1D[i,:] = Molecules[:,6]
pinchDrugcoo_1D[i,:] = Molecules[:,7]
pinchCH_1D[i,:] = Molecules[:,8]
Col_1D[:,t] = pinchCol_1D[:,ntt-1]
Colcoo_1D[:,t] = pinchColcoo_1D[:,ntt-1]
Cm_1D[:,t] = pinchCm_1D[:,ntt-1]
Cmcoo_1D[:,t] = pinchCmcoo_1D[:,ntt-1]
Cpol_1D[:,t+1] = pinchCpol_1D[:,ntt-1]
Cpolcoo_1D[:,t+1] = pinchCpolcoo_1D[:,ntt-1]
Drug_1D[:,t] = pinchDrug_1D[:,ntt-1]
Drugcoo_1D[:,t] = pinchDrugcoo_1D[:,ntt-1]
CH_1D[:,t] = pinchCH_1D[:,ntt-1]
return(Drug_1D,Drugcoo_1D,Col_1D,Colcoo_1D,CH_1D,COH_1D,CH2O_1D,Cpol_1D,Cpolcoo_1D,Cm_1D,Cmcoo_1D)
def ScissionRate(
Ce_1D,Cend_1D,Res_1D,Rrs_1D,Rol_1D,Rm_1D,Col_1D,Cm_1D,Xc_1D,Xext_1D,
Xcend_1D,Mn_1D,Mask_1D,invMask_1D,Cpol0,Ce0,t,dt ,CH2O_1D,CSol_1D,Sol_conc,DegWeight ,rb,dx,CeSwell_1D,CH_1D):
# =============================================================================
# Recalculate the rate of scission for new time step based on C_H+
# =============================================================================
Xext_1D[:,t+1] = Mask_1D[:]*(p*etaA*(Rrs_1D[:,t]+Res_1D[:,t])*Vc)
if lam == 1:
Xc_1D[:,t+1] = Mask_1D[:]*(Xmax - (Xmax - Xc0)*np.exp(-Xext_1D[:,t+1]))
else:
Xc_1D[:,t+1] = Mask_1D[:]*(Xmax - ((lam-1)*Xext_1D[:,t+1] + (Xmax-Xc0)**(1-lam))**(1/(1-lam)))
omega = Ce0
# Ce_1D[:,t] = Mask_1D[:]*(Ce0 - Rol_1D[:,t] - Rm_1D[:,t] - omega*(Xc_1D[:,t]-Xc0) )
for index in range(0,nx-1):
if Ce_1D[index,t] <= 0.5*Cend_1D[index,t]:
Xcend_1D[index,t] = (Xcend_1D[index,t] + Cend_1D[index,t])
Cend_1D[index,t] = 0
Ce_1D[index,t] = 0
Cm_1D[index,t] = Cm_1D[index,t] + 2 *Mask_1D[index]
Rm_1D[index,t] = Rm_1D[index,t] + 2 *Mask_1D[index]
Mn_1D[index,t+1] = Mask_1D[index] *( M0 + omega*Xc_1D[index,t])/(Xcend_1D[index,t] + invMask_1D[index])
else:
Cend_1D[index,t] = Mask_1D[index]* 2*( Cpol0 + (Rrs_1D[index,t] - Rol_1D[index,t]/m))
Mn_1D[index,t+1] = Mask_1D[index] * ((Ce_1D[index,t]+omega*Xc_1D[index,t]) * M0)/(0.5*Cend_1D[index,t] + invMask_1D[index])
for index in range(0,nx-1):
if np.isnan(Ce_1D[index,t])==1: Ce_1D[index,t] = 0
if np.isnan(Cend_1D[index,t])==1: Cend_1D[index,t] = 0
Rrs_temp = np.zeros((nx,nttt))
Res_temp = np.zeros((nx,nttt))
Ce_temp = np.zeros((nx,nttt))
Cend_temp = np.zeros((nx,nttt))
CH_temp = np.zeros((nx,nttt))
Rrs_temp[:,0] = Rrs_1D[:,t]
Res_temp[:,0] = Res_1D[:,t]
Ce_temp[:,0] = Ce_1D[:,t]
Cend_temp[:,0] = Cend_1D[:,t]
CH_temp[:,0] = CH_1D[:,t]
DegWeightSciss[:] = Maskfit_1D[:,t]
for ttt in range(0,nttt-1):
dttt = dt/nttt
Rrs_temp[:,ttt+1] = Rrs_temp[:,ttt] + DegWeightSciss[:]*(dttt*(Ce_1D[:,t] * (kr1 + kr2*CH_1D[:,t])))
Res_temp[:,ttt+1] = Res_temp[:,ttt] + DegWeightSciss[:]*(dttt*(Cend_1D[:,t]*(ke1 + ke2*CH_1D[:,t])))
Rrs_1D[:,t+1] = Rrs_temp[:,-1]
Res_1D[:,t+1] = Res_temp[:,-1]
Roltemp = np.zeros(nx)
for M in range(1,9):
Roltemp[:] = Roltemp[:] + (Ce0*M*(Rrs_1D[:,t+1]/Ce0)**2*(1-(Rrs_1D[:,t+1]/Ce0))**(M-1) )
Rol_1D[:,t+1] = Roltemp[:]
Rm_1D[:,t+1] = Res_1D[:,t+1]
Ce_1D[:,t+1] = Ce_1D[:,t] + dt*(
-m*(Rol_1D[:,t+1]-Rol_1D[:,t])/dt
-(Rm_1D[:,t+1]-Rm_1D[:,t])/dt
-(Ce_1D[:,t]/(1-Xc_1D[:,t+1]) * (Xc_1D[:,t+1]-Xc_1D[:,t])/dt )
)
Col_1D[:,t+1] = Col_1D[:,t]+ dt*(Rol_1D[:,t+1]-Rol_1D[:,t])/dt
Cm_1D[:,t+1] = Cm_1D[:,t] + dt*(Rm_1D[:,t+1]-Rm_1D[:,t])/dt
return(Ce_1D,Cend_1D,Res_1D,Rrs_1D,Rol_1D,Rm_1D,Col_1D,Cm_1D,Xc_1D,Xext_1D,
Xcend_1D,Mn_1D,Mask_1D,invMask_1D,Cpol0,t,DegWeight,omega)
def Diffusion(
Drug_1D,Drugcoo_1D,Col_1D,Colcoo_1D,CH_1D,COH_1D,CH2O_1D,Cmcoo_1D,Cm_1D,
DDrug_1D,DCol_1D,DCH_1D,DCm_1D,t,dt ,A,rb,dx,Rol_1D,Rm_1D ):
r = np.linspace(0,Lx,nx)
sig_xRC = np.zeros(nx)
sig_xLC = np.zeros(nx)
sig_xyC = np.zeros(nx)
rDCol_1D = np.zeros(nx)
rDDrug_1D = np.zeros(nx)
rDCH_1D = np.zeros(nx)
rDCm_1D = np.zeros(nx)
rDCol_1D[:] = DCol_1D[:,t]
rDDrug_1D[:] = DDrug_1D[:,t]
rDCH_1D[:] = DCH_1D[:,t]
rDCm_1D[:] = DCm_1D[:,t]
# =============================================================================
# Diffusion
# =============================================================================
w1 = 1#2/3
w2 = 1#4/3
w3 = 0#1/3
r = np.linspace(0,Lx,nx)#/Lx
dr = r[1]-r[0]
I = np.linspace(1,nx,nx)
sig_xRC[2:nx] = -(w1)*(dt/(dr)**2)*(rDCol_1D[1:nx-1]/2/I[1:nx-1]*I[2:nx]+rDCol_1D[2:nx]/2 )
sig_xLC[0:nx-2] = -(w1)*(dt/(dr)**2)*(rDCol_1D[1:nx-1]/2/I[1:nx-1]*I[:nx-2]+rDCol_1D[0:nx-2]/2 )
sig_xyC[1:nx-1] = 1 + (w1)*(dt/(dr)**2)*(rDCol_1D[1:nx-1] + rDCol_1D[2:nx]/2+ rDCol_1D[:nx-2]/2)
sig_xyC[0] = sig_xyC[-1] = 1
A1 = sparse.spdiags([sig_xyC[:],sig_xRC[:],sig_xLC[:]],[0,1,-1],nx,nx)
B1 = sparse.spdiags([(w2)*np.ones(nx)],[0],nx,nx)
C1 = sparse.spdiags([(w3)*np.ones(nx)],[0],nx,nx)
sig_xRC[2:nx] = -(w1)*(dt/(dr)**2)*(rDDrug_1D[1:nx-1]/2/I[1:nx-1]*I[2:nx]+rDDrug_1D[2:nx]/2 )
sig_xLC[0:nx-2] = -(w1)*(dt/(dr)**2)*(rDDrug_1D[1:nx-1]/2/I[1:nx-1]*I[:nx-2]+rDDrug_1D[0:nx-2]/2 )
sig_xyC[1:nx-1] = 1 + (w1)*(dt/(dr)**2)*(rDDrug_1D[1:nx-1] + rDDrug_1D[2:nx]/2+ rDDrug_1D[:nx-2]/2)
sig_xyC[0] = sig_xyC[-1] = 1
A2 = sparse.spdiags([sig_xyC[:],sig_xRC[:],sig_xLC[:]],[0,1,-1],nx,nx)
B2 = sparse.spdiags([(w2)*np.ones(nx)],[0],nx,nx)
C2 = sparse.spdiags([(w3)*np.ones(nx)],[0],nx,nx)
sig_xRC[2:nx] = -(w1)*(dt/(dr)**2)*(rDCH_1D[1:nx-1]/2/I[1:nx-1]*I[2:nx]+rDCH_1D[2:nx]/2 )
sig_xLC[0:nx-2] = -(w1)*(dt/(dr)**2)*(rDCH_1D[1:nx-1]/2/I[1:nx-1]*I[:nx-2]+rDCH_1D[0:nx-2]/2 )
sig_xyC[1:nx-1] = 1 + (w1)*(dt/(dr)**2)*(rDCH_1D[1:nx-1] + rDCH_1D[2:nx]/2+ rDCH_1D[:nx-2]/2)
sig_xyC[0] = sig_xyC[-1] = 1
A3 = sparse.spdiags([sig_xyC[:],sig_xRC[:],sig_xLC[:]],[0,1,-1],nx,nx)
B3 = sparse.spdiags([(w2)*np.ones(nx)],[0],nx,nx)
C3 = sparse.spdiags([(w3)*np.ones(nx)],[0],nx,nx)
sig_xRC[2:nx] = -(w1)*(dt/(dr)**2)*(rDCm_1D[1:nx-1]/2/I[1:nx-1]*I[2:nx]+rDCm_1D[2:nx]/2 )
sig_xLC[0:nx-2] = -(w1)*(dt/(dr)**2)*(rDCm_1D[1:nx-1]/2/I[1:nx-1]*I[:nx-2]+rDCm_1D[0:nx-2]/2 )
sig_xyC[1:nx-1] = 1 + (w1)*(dt/(dr)**2)*(rDCm_1D[1:nx-1] + rDCm_1D[2:nx]/2+ rDCm_1D[:nx-2]/2)
sig_xyC[0] = sig_xyC[-1] = 1
A4 = sparse.spdiags([sig_xyC[:],sig_xRC[:],sig_xLC[:]],[0,1,-1],nx,nx)
B4 = sparse.spdiags([(w2)*np.ones(nx)],[0],nx,nx)
C4 = sparse.spdiags([(w3)*np.ones(nx)],[0],nx,nx)
# =============================================================================
# Solving The System of PDEs
# =============================================================================
bdrug = (B2.dot(Drug_1D[:,t]))
cdrug = (C2.dot(Drug_1D[:,t-1]))
Bdrug = bdrug-cdrug
bcol = (B1.dot(Col_1D[:,t]))
ccol = (C1.dot(Col_1D[:,t-1]))
Bcol = bcol-ccol
bcm = (B4.dot(Cm_1D[:,t]))
ccm = (C4.dot(Cm_1D[:,t-1]))
Bcm = bcm-ccm
bdrugcoo = (B2.dot(Drugcoo_1D[:,t]))
cdrugcoo = (C2.dot(Drugcoo_1D[:,t-1]))
Bdrugcoo = bdrugcoo-cdrugcoo
bcolcoo = (B1.dot(Colcoo_1D[:,t]))
ccolcoo = (C1.dot(Colcoo_1D[:,t-1]))
Bcolcoo = bcolcoo-ccolcoo
bcmcoo = (B4.dot(Cmcoo_1D[:,t]))
ccmcoo = (C4.dot(Cmcoo_1D[:,t-1]))
Bcmcoo = bcmcoo-ccmcoo
bch = (B3.dot(CH_1D[:,t]))
cch = (C3.dot(CH_1D[:,t-1]))
Bch = bch-cch
bch2o = (B3.dot(CH2O_1D[:,t]))
cch2o = (C3.dot(CH2O_1D[:,t-1]))
Bch2o = bch2o - cch2o
bcoh = (B3.dot(COH_1D[:,t]))
ccoh = (C3.dot(COH_1D[:,t-1]))
Bcoh = bcoh - ccoh
# =============================================================================
tempd= spsolve(A2, Bdrug)
tempc= spsolve(A1, Bcol)
tempm= spsolve(A4, Bcm)
tempD = spsolve(A2, Bdrugcoo)
tempC = spsolve(A1, Bcolcoo)
tempM = spsolve(A4, Bcmcoo)
temph= spsolve(A3, Bch)
tempH = spsolve(A3, Bcoh)
Drug_1D[:,t+1] = tempd[:]
Col_1D[:,t+1] = tempc[:]
Cm_1D[:,t+1] = tempm[:]
Drugcoo_1D[:,t+1] = tempD[:]
Colcoo_1D[:,t+1] = tempC[:]
Cmcoo_1D[:,t+1] = tempM[:]
CH_1D[:,t+1] = temph[:]
COH_1D[:,t+1] = tempH[:]
Col_1D[0,t+1] = Col_1D[1,t+1]
Cm_1D[0,t+1] = Cm_1D[1,t+1]
Drugcoo_1D[0,t+1] = Drugcoo_1D[1,t+1]
Colcoo_1D[0,t+1] = Colcoo_1D[1,t+1]
Cmcoo_1D[0,t+1] = Cmcoo_1D[1,t+1]
CH_1D[0,t+1] = CH_1D[1,t+1]
COH_1D[0,t+1] = COH_1D[1,t+1]
CH_1D[-1,t+1] = CH0
return(Drug_1D,Drugcoo_1D,Col_1D,Colcoo_1D,CH_1D,COH_1D,CH2O_1D,Cmcoo_1D,Cm_1D,
DDrug_1D,DCol_1D,DCH_1D,DCm_1D,t,A)
def Porosity( Maskfit_1D,invMaskfit_1D,VporeCol_1D,VporeDrug_1D,Vpore_1D,Rol_1D,Rm_1D,
Col_1D,Colcoo_1D,Cm_1D,Cmcoo_1D,Drug_1D,Drugcoo_1D,DCol_1D,DCm_1D,DDrug_1D,
DCH_1D,t,fDrug,Ce0,dt ,rb,dx ,CeSwell_1D,Mask_1D,Xc_1D,Xcend_1D,Xext_1D):
VporeCol_1D[:,t+1]= Maskfit_1D[:,t+1]*(
(( (Rol_1D[:,t+1]+Rm_1D[:,t+1])/Ce0)
- (Col_1D[:,t+1]+Colcoo_1D[:,t+1]-Col0
+ Cm_1D[:,t+1]+Cmcoo_1D[:,t+1]-Cm0) /Ce0 )
* (1-Xc0) + Xc_1D[:,t] )
VporeDrug_1D[:,t+1]= Maskfit_1D[:,t+1]*( 1 - (Drug_1D[:,t+1]+Drugcoo_1D[:,t+1])/((Drug_1D[:,0]+Drugcoo_1D[:,0])+tisMask_1D[:]))
Vpore_1D[:,t+1]=(VporeCol_1D[:,t+1])
Dpolymer_Col_Swell = Dpolymer_Col * np.exp( -Beta*(1- (CeSwell_1D[t+1])))
Dpolymer_Drug_Swell= Dpolymer_Drug* np.exp( -Beta*(1- (CeSwell_1D[t+1])))
Dpolymer_Cm_Swell = Dpolymer_Cm * np.exp( -Beta*(1- (CeSwell_1D[t+1])))
Dpolymer_CH_Swell = Dpolymer_CH * np.exp( -Beta*(1- (CeSwell_1D[t+1])))
if Dpolymer_Col_Swell > Dpore_Col: Dpolymer_Col_Swell = Dpore_Col
if Dpolymer_Drug_Swell > Dpore_Drug: Dpolymer_Drug_Swell = Dpore_Drug
if Dpolymer_Cm_Swell > Dpore_Cm: Dpolymer_Cm_Swell = Dpore_Cm
if Dpolymer_CH_Swell > Dpore_CH: Dpolymer_CH_Swell = Dpore_CH
DCol_1D[:,t+1] = invMaskfit_1D[:,t+1]*Dpore_Col + Maskfit_1D[:,t+1]*(Dpolymer_Col_Swell+(1.3*Vpore_1D[:,t+1]**2-0.3*Vpore_1D[:,t+1]**3)*(Dpore_Col-Dpolymer_Col_Swell))
DCm_1D[:,t+1] = invMaskfit_1D[:,t+1]*Dpore_Cm + Maskfit_1D[:,t+1]*(Dpolymer_Cm_Swell+(1.3*Vpore_1D[:,t+1]**2-0.3*Vpore_1D[:,t+1]**3)*(Dpore_Cm-Dpolymer_Cm_Swell))
DCH_1D[:,t+1] = invMaskfit_1D[:,t+1]*Dpore_CH + Maskfit_1D[:,t+1]*(Dpolymer_CH_Swell+(1.3*Vpore_1D[:,t+1]**2-0.3*Vpore_1D[:,t+1]**3)*(Dpore_CH-Dpolymer_CH_Swell))
DDrug_1D[:,t+1] = invMaskfit_1D[:,t+1]*Dpore_Drug + Maskfit_1D[:,t+1]*(Dpolymer_Drug_Swell+(1.3*Vpore_1D[:,t+1]**2-0.3*Vpore_1D[:,t+1]**3)*(Dpore_Drug-Dpolymer_Drug_Swell))
DDrug_1D[:,t+1] = np.copy(DDrug_1D[:,t+1])*(1-tisMask_1D) + Dpore_Drug_Tissue*(tisMask_1D)
return(Maskfit_1D,invMaskfit_1D,VporeCol_1D,VporeDrug_1D,Vpore_1D,Rol_1D,Rm_1D,
Col_1D,Colcoo_1D,Cm_1D,Cmcoo_1D,Drug_1D,Drugcoo_1D,DCol_1D,DCm_1D,DDrug_1D,
DCH_1D,t,fDrug)
def UpdateDiffusion(
Drug_1D,Drugcoo_1D,Col_1D,Colcoo_1D,CH_1D,COH_1D,CH2O_1D,Cmcoo_1D,Cm_1D,
DDrug_1D,DCol_1D,DCH_1D,DCm_1D,t,dt,rb,dx,invtisMask_1D,tisMask_1D,Rol_1D,Rm_1D ):
r = np.linspace(0,Lx,nx)
sig_xRC = np.zeros(nx)
sig_xLC = np.zeros(nx)
sig_xyC = np.zeros(nx)
rDCol_1D = np.zeros(nx)
rDDrug_1D = np.zeros(nx)
rDCH_1D = np.zeros(nx)
rDCm_1D = np.zeros(nx)
rDCol_1D[:] = DCol_1D[:,t+1]
rDDrug_1D[:] = DDrug_1D[:,t+1]
rDCH_1D[:] = DCH_1D[:,t+1]
rDCm_1D[:] = DCm_1D[:,t+1]
# =============================================================================
# Diffusion
# =============================================================================
w1 = 1# 2/3
w2 =1# 4/3
w3 = 0# 1/3
r = np.linspace(0,Lx,nx)#/Lx
dr = r[1]-r[0]
I = np.linspace(1,nx,nx)
sig_xRC[2:nx] = -(w1)*(dt/(dr)**2)*(rDCol_1D[1:nx-1]/2/I[1:nx-1]*I[2:nx]+rDCol_1D[2:nx]/2 )
sig_xLC[0:nx-2] = -(w1)*(dt/(dr)**2)*(rDCol_1D[1:nx-1]/2/I[1:nx-1]*I[:nx-2]+rDCol_1D[0:nx-2]/2 )
sig_xyC[1:nx-1] = 1 + (w1)*(dt/(dr)**2)*(rDCol_1D[1:nx-1] + rDCol_1D[2:nx]/2+ rDCol_1D[:nx-2]/2)
sig_xyC[0] = sig_xyC[-1] = 1
A1 = sparse.spdiags([sig_xyC[:],sig_xRC[:],sig_xLC[:]],[0,1,-1],nx,nx)
B1 = sparse.spdiags([(w2)*np.ones(nx)],[0],nx,nx)
C1 = sparse.spdiags([(w3)*np.ones(nx)],[0],nx,nx)
sig_xRC[2:nx] = -(w1)*(dt/(dr)**2)*(rDDrug_1D[1:nx-1]/2/I[1:nx-1]*I[2:nx]+rDDrug_1D[2:nx]/2 )
sig_xLC[0:nx-2] = -(w1)*(dt/(dr)**2)*(rDDrug_1D[1:nx-1]/2/I[1:nx-1]*I[:nx-2]+rDDrug_1D[0:nx-2]/2 )
sig_xyC[1:nx-1] = 1 + (w1)*(dt/(dr)**2)*(rDDrug_1D[1:nx-1] + rDDrug_1D[2:nx]/2+ rDDrug_1D[:nx-2]/2)
sig_xyC[0] = 1 + sig_xRC[2]#dt*(rDDrug_1D[0,t+1]+rDDrug_1D[1,t+1])/2*6*dt/dr**2
sig_xRC[1] = -sig_xRC[2] #dt*(rDDrug_1D[0,t+1]+rDDrug_1D[1,t+1])/2*6*dt/dr**2
sig_xyC[-1] = 1
# if t ==1: A = np.diag(sig_xLC[:-1],-1)+np.diag(sig_xyC,0)+np.diag(sig_xRC[1:],1)
A2 = sparse.spdiags([sig_xyC[:],sig_xRC[:],sig_xLC[:]],[0,1,-1],nx,nx)
B2 = sparse.spdiags([(w2)*np.ones(nx)],[0],nx,nx)
C2 = sparse.spdiags([(w3)*np.ones(nx)],[0],nx,nx)
sig_xRC[2:nx] = -(w1)*(dt/(dr)**2)*(rDCH_1D[1:nx-1]/2/I[1:nx-1]*I[2:nx]+rDCH_1D[2:nx]/2 )
sig_xLC[0:nx-2] = -(w1)*(dt/(dr)**2)*(rDCH_1D[1:nx-1]/2/I[1:nx-1]*I[:nx-2]+rDCH_1D[0:nx-2]/2 )
sig_xyC[1:nx-1] = 1 + (w1)*(dt/(dr)**2)*(rDCH_1D[1:nx-1] + rDCH_1D[2:nx]/2+ rDCH_1D[:nx-2]/2)
sig_xyC[0] = sig_xyC[-1] = 1
A3 = sparse.spdiags([sig_xyC[:],sig_xRC[:],sig_xLC[:]],[0,1,-1],nx,nx)
B3 = sparse.spdiags([(w2)*np.ones(nx)],[0],nx,nx)
C3 = sparse.spdiags([(w3)*np.ones(nx)],[0],nx,nx)
sig_xRC[2:nx] = -(w1)*(dt/(dr)**2)*(rDCm_1D[1:nx-1]/2/I[1:nx-1]*I[2:nx]+rDCm_1D[2:nx]/2 )
sig_xLC[0:nx-2] = -(w1)*(dt/(dr)**2)*(rDCm_1D[1:nx-1]/2/I[1:nx-1]*I[:nx-2]+rDCm_1D[0:nx-2]/2 )
sig_xyC[1:nx-1] = 1 + (w1)*(dt/(dr)**2)*(rDCm_1D[1:nx-1] + rDCm_1D[2:nx]/2+ rDCm_1D[:nx-2]/2)
sig_xyC[0] = sig_xyC[-1] = 1
A4 = sparse.spdiags([sig_xyC[:],sig_xRC[:],sig_xLC[:]],[0,1,-1],nx,nx)
B4 = sparse.spdiags([(w2)*np.ones(nx)],[0],nx,nx)
C4 = sparse.spdiags([(w3)*np.ones(nx)],[0],nx,nx)
# =============================================================================
# Solving The System of PDEs
# =============================================================================
bdrug = (B2.dot(Drug_1D[:,t]))
cdrug = (C2.dot(Drug_1D[:,t-1]))
Bdrug = bdrug-cdrug
bcol = (B1.dot(Col_1D[:,t]))
ccol = (C1.dot(Col_1D[:,t-1]))
Bcol = bcol-ccol
bcm = (B4.dot(Cm_1D[:,t]))
ccm = (C4.dot(Cm_1D[:,t-1]))
Bcm = bcm-ccm
bdrugcoo = (B2.dot(Drugcoo_1D[:,t]))
cdrugcoo = (C2.dot(Drugcoo_1D[:,t-1]))
Bdrugcoo = bdrugcoo-cdrugcoo
bcolcoo = (B1.dot(Colcoo_1D[:,t]))
ccolcoo = (C1.dot(Colcoo_1D[:,t-1]))
Bcolcoo = bcolcoo-ccolcoo
bcmcoo = (B4.dot(Cmcoo_1D[:,t]))
ccmcoo = (C4.dot(Cmcoo_1D[:,t-1]))
Bcmcoo = bcmcoo-ccmcoo
bch = (B3.dot(CH_1D[:,t]))
cch = (C3.dot(CH_1D[:,t-1]))
Bch = bch-cch
bch2o = (B3.dot(CH2O_1D[:,t]))
cch2o = (C3.dot(CH2O_1D[:,t-1]))
Bch2o = bch2o - cch2o
bcoh = (B3.dot(COH_1D[:,t]))
ccoh = (C3.dot(COH_1D[:,t-1]))
Bcoh = bcoh - ccoh
# =============================================================================
tempd= spsolve(A2, Bdrug)
tempc= spsolve(A1, Bcol)
tempm= spsolve(A4, Bcm)
tempD = spsolve(A2, Bdrugcoo)
tempC = spsolve(A1, Bcolcoo)
tempM = spsolve(A4, Bcmcoo)
temph= spsolve(A3, Bch)
tempH = spsolve(A3, Bcoh)
Drug_1D[:,t+1] = tempd[:]
Col_1D[:,t+1] = tempc[:]
Cm_1D[:,t+1] = tempm[:]
Drugcoo_1D[:,t+1] = tempD[:]
Colcoo_1D[:,t+1] = tempC[:]
Cmcoo_1D[:,t+1] = tempM[:]
COH_1D[:,t+1] = tempH[:]
Col_1D[0,t+1] = Col_1D[1,t+1]
Cm_1D[0,t+1] = Cm_1D[1,t+1]
Colcoo_1D[0,t+1] = Colcoo_1D[1,t+1]
Cmcoo_1D[0,t+1] = Cmcoo_1D[1,t+1]
COH_1D[0,t+1] = COH_1D[1,t+1]
CH_1D[:,t+1] = temph[:] # * invtisMask_1D[:] + CH0*tisMask_1D[:]
CH_1D[0,t+1] = CH_1D[1,t+1]
CH_1D[-1,t+1] = CH0
return(Drug_1D,Drugcoo_1D,Col_1D,Colcoo_1D,CH_1D,COH_1D,CH2O_1D,Cmcoo_1D,Cm_1D,
DDrug_1D,DCol_1D,DCH_1D,DCm_1D,t)
def Redistribution(Ce_1D,Col_1D,Colcoo_1D,Cm_1D,Cmcoo_1D,Cpol_1D,Cpolcoo_1D,Res_1D,Rrs_1D,Rol_1D,Rm_1D,Cend_1D
,Drug_1D,Drugcoo_1D,CSol_1D,CH2O_1D,CH_1D,COH_1D,t,rb,rb0,Cpol0,Ce0,Rb,CeSwell_1D):
OCopy_Ce_1D[:] = np.copy(Ce_1D[:,t-1]); OCopy_Col_1D[:] = np.copy(Col_1D[:,t-1]); OCopy_Colcoo_1D[:] = np.copy(Colcoo_1D[:,t-1])
OCopy_Cm_1D[:] = np.copy(Cm_1D[:,t-1]); OCopy_Cmcoo_1D[:] = np.copy(Cmcoo_1D[:,t-1]); OCopy_Cpol_1D[:] = np.copy(Cpol_1D[:,t-1]); OCopy_Cpolcoo_1D[:] = np.copy(Cpolcoo_1D[:,t-1])
OCopy_Res_1D[:] = np.copy(Res_1D[:,t-1]); OCopy_Rrs_1D[:] = np.copy(Rrs_1D[:,t-1]); OCopy_Rol_1D[:] = np.copy(Rol_1D[:,t-1]); OCopy_Rm_1D[:] = np.copy(Rm_1D[:,t-1])
OCopy_Cend_1D[:] = np.copy(Cend_1D[:,t-1]); OCopy_Drug_1D[:] = np.copy(Drug_1D[:,t-1]); OCopy_Drugcoo_1D[:] = np.copy(Drugcoo_1D[:,t-1]); OCopy_CSol_1D[:] = np.copy(CSol_1D[:,t-1])
OCopy_CH2O_1D[:] = np.copy(CH2O_1D[:,t-1]); OCopy_CH_1D[:] = np.copy(CH_1D[:,t-1]); OCopy_COH_1D[:] = np.copy(COH_1D[:,t-1]); OCopy_CSol_1D[:] = np.copy(CSol_1D[:,t-1])
Copy_Ce_1D[:] = np.copy(Ce_1D[:,t]); Copy_Col_1D[:] = np.copy(Col_1D[:,t]); Copy_Colcoo_1D[:] = np.copy(Colcoo_1D[:,t]); Copy_Cm_1D[:] = np.copy(Cm_1D[:,t])
Copy_Cmcoo_1D[:] = np.copy(Cmcoo_1D[:,t]); Copy_Cpol_1D[:] = np.copy(Cpol_1D[:,t]); Copy_Cpolcoo_1D[:] = np.copy(Cpolcoo_1D[:,t]); Copy_Res_1D[:] = np.copy(Res_1D[:,t])
Copy_Rrs_1D[:] = np.copy(Rrs_1D[:,t]); Copy_Rol_1D[:] = np.copy(Rol_1D[:,t]); Copy_Rm_1D[:] = np.copy(Rm_1D[:,t]); Copy_Cend_1D[:] = np.copy(Cend_1D[:,t])
Copy_Drug_1D[:] = np.copy(Drug_1D[:,t]); Copy_Drugcoo_1D[:] = np.copy(Drugcoo_1D[:,t]); Copy_CSol_1D[:] = np.copy(CSol_1D[:,t]); Copy_CH2O_1D[:] = np.copy(CH2O_1D[:,t])
Copy_CH_1D[:] = np.copy(CH_1D[:,t]); Copy_COH_1D[:] = np.copy(COH_1D[:,t]); Copy_CSol_1D[:] = np.copy(CSol_1D[:,t])
rbp = Rb[t-1]; rb = Rb[t]
Vb = (4/3*np.pi*rb**3); Vbp = (4/3*np.pi*rbp**3)
Vsi = (4/3*np.pi)*((rb0**3) - (rb0-rs)**3); Vsf = (4/3*np.pi)*((rb**3) - (rb - rs)**3)
CeSwell_1D[t] = Vsf/Vsi
Ce_1D[:,t] = 1/CeSwell_1D[t] * Copy_Ce_1D[:]; Col_1D[:,t] = Vbp/Vb * Copy_Col_1D[:]
Colcoo_1D[:,t] = Vbp/Vb * Copy_Colcoo_1D[:]; Cm_1D[:,t] = Vbp/Vb * Copy_Cm_1D[:]
Cmcoo_1D[:,t] = Vbp/Vb * Copy_Cmcoo_1D[:]; Cpol_1D[:,t] = Vbp/Vb * Copy_Cpol_1D[:]
Cpolcoo_1D[:,t] = Vbp/Vb * Copy_Cpolcoo_1D[:]; Drug_1D[:,t] = Vbp/Vb * Copy_Drug_1D[:]
Drugcoo_1D[:,t] = Vbp/Vb * Copy_Drugcoo_1D[:]; CSol_1D[:,t] = Vbp/Vb * Copy_CSol_1D[:]
CH2O_1D[:,t] = Vbp/Vb * Copy_CH2O_1D[:]; CH_1D[:,t] = Vbp/Vb * Copy_CH_1D[:]; CSol_1D[:,t] = Vbp/Vb * Copy_CSol_1D[:]
Ce_1D[:,t-1] = 1/CeSwell_1D[t] * OCopy_Ce_1D[:]; Col_1D[:,t-1] = Vbp/Vb * OCopy_Col_1D[:]
Colcoo_1D[:,t-1] = Vbp/Vb * OCopy_Colcoo_1D[:]; Cm_1D[:,t-1] = Vbp/Vb * OCopy_Cm_1D[:]
Cmcoo_1D[:,t-1] = Vbp/Vb * OCopy_Cmcoo_1D[:]; Cpol_1D[:,t-1] = Vbp/Vb * OCopy_Cpol_1D[:] ;Cpolcoo_1D[:,t-1] = Vbp/Vb * OCopy_Cpolcoo_1D[:]
Drug_1D[:,t-1] = Vbp/Vb * OCopy_Drug_1D[:]; Drugcoo_1D[:,t-1] = Vbp/Vb * OCopy_Drugcoo_1D[:]
CSol_1D[:,t-1] = Vbp/Vb * OCopy_CSol_1D[:]; CH_1D[:,t-1] = Vbp/Vb * OCopy_CH_1D[:]; CSol_1D[:,t-1] = Vbp/Vb * OCopy_CSol_1D[:]
return(Ce_1D,Col_1D,Colcoo_1D,Cm_1D,Cmcoo_1D,Cpol_1D,Cpolcoo_1D,Res_1D,Rrs_1D,Rol_1D,Rm_1D,Cend_1D
,Drug_1D,Drugcoo_1D,CSol_1D,CH2O_1D,CH_1D,COH_1D,DrugRelease_1D)
(DDrug_1D,Drug_1D,Drugcoo_1D,Maskfit_1D,
CSol_1D,CH2O_1D,DDrug_1D,A,omega,Ce0,invtisMask_1D,phi_SM,phi_PM,phi_WM,CH_1D,
nstart,nend,Ce_1D,Cend_1D,Xc_1D ,Rb,dX, TimeR, TimeT, DrugRelease_1D,CeSwell_1D,
Rrs_1D,Res_1D,Vpore_1D,Mn_1D,invMaskfit_1D,DCH_1D,
Cpol_1D, Cpolcoo_1D, Col_1D, Colcoo_1D, Cm_1D, Cmcoo_1D, Rol_1D,Rm_1D
)= ImplantModel(Variables)
N = len(Drug_1D[0,:])
DrugRel = np.zeros((nx,N))
for n in range(0,N):
DrugRel[:,n] = (Drug_1D[:,n]*Mask_1D) + (Drugcoo_1D[:,n]*Mask_1D)
SumDrug = np.sum(DrugRel,axis=0)
ProfileDrug = 1-SumDrug/SumDrug[0]
ProfileMn = np.zeros(N)
#
#np.savetxt('Drug_1D.txt',Drug_1D)
#np.savetxt('Drugcoo_1D.txt',Drugcoo_1D)
#np.savetxt('CH_1D.txt',CH_1D)
#np.savetxt('Maskfit_1D.txt',Maskfit_1D)
#np.savetxt('TimeR.txt',TimeR)
#np.savetxt('TimeT.txt',TimeT)
#np.savetxt('ProfileDrug.txt',ProfileDrug)
#
def PlotDrugF():
plt.figure(2,figsize=(8,6))
if ((Degr == 1)): plt.plot(TimeR[1:-1],(ProfileDrug[1:-1]),color='navy',linestyle='-',linewidth=2.5)
if ((Degr == 0)): plt.plot(TimeR[1:-1],(ProfileDrug[1:-1]),color='grey',linestyle='--',linewidth=2.5)
plt.errorbar(Data[0,:],Data[1,:]/100,yerr=Data[2,:]/100,color='k',marker='o',linestyle='',capsize=3,ecolor='gray')
plt.xlim((-1,(Data[0,-1]+1)))
plt.ylim((0,1))
plt.xlabel('Time (Days)',fontsize=20)
plt.ylabel('Drug Released (normalized)',fontsize=20)
def PlotDrug():
plt.figure(1,figsize=(8,6))
plt.plot(TimeR[1:],(ProfileDrug[:-1]),label='MW: %2.0f kDa' %(int(MW/1000))) #,color='royalblue',linestyle='--',linewidth=2.5)
# plt.plot(TimeR[1:],(ProfileDrug[:-1]))#,color='navy',linestyle='-',linewidth=2.5)
plt.errorbar(Data[0,:],Data[1,:]/100,yerr=Data[2,:]/100,color='k',marker='o',linestyle='',capsize=3,ecolor='gray')
#plt.legend(loc='best')
plt.xlim((-1,(Data[0,-1]+1)))
plt.ylim((0,1))
plt.xlabel('Time (Days)',fontsize=20)
plt.ylabel('Drug Released (normalized)',fontsize=20)
PlotDrug()
for ij in range(0,5):
print('\007')
time.sleep(1)
nt = len(TimeR)
MassTotal = np.zeros((nx,nt))
MassPoly = np.zeros((nx,nt))
MassDrug = np.zeros((nx,nt))
MassSol = np.zeros((nx,nt))
MassTotalSum = np.zeros(nt)
MassTotalPoly = np.zeros(nt)
MassTotalDrug = np.zeros(nt)
MassTotalSol = np.zeros(nt)
Vimplant = (4/3*np.pi*rb0**3)
PoreV = (4/3*np.pi*(rb-rs)**3)*0.74048 * (rp)**3/(rp+0.5*thick)**3
PolymerV = (4/3*np.pi*(rb-rs)**3) - PoreV + 4/3*np.pi*(rb**3 - (rb-rs)**3)
for n in range(0,nt):
MassPoly[:,n] = (
invtisMask_1D[:]*0*( (Col_1D[:,n]+Colcoo_1D[:,n])
+(Cm_1D[:,n]+Cmcoo_1D[:,n])
)/(np.sum(invtisMask_1D[:]))*Vimplant*MWpolyunit
+ (Ce_1D[:,n] )/(np.sum(Mask_1D[:]))*(Vimplant*(PolymerV/(PolymerV+PoreV)))*MWpolyunit)
MassDrug[:,n] = invtisMask_1D[:]*((Drug_1D[:,n]+Drugcoo_1D[:,n])
)/(np.sum(invtisMask_1D[:]))*Vimplant*Mw_d
MassSol[:,n] = invtisMask_1D[:]*(CSol_1D[:,n]
)/(np.sum(invtisMask_1D[:]))*Vimplant*Mw_s
MassTotal[:,n] = (MassPoly[:,n] + MassDrug[:,n] +MassSol[:,n]
)
MassTotalPoly = np.sum(MassPoly[:,:],axis=0)
MassTotalDrug = np.sum(MassDrug[:,:],axis=0)
MassTotalSol = np.sum(MassSol[:,:],axis=0)
MassTotalSum = (np.sum(MassTotal[:,:],axis=0))/np.sum(MassTotal[:,0],axis=0)
Times = [0,1,2.5,5,7.5,10,12.5,15,17.5,20,22.5]; Tim = 0
if Lt >=30: Times = [0,1,2.5,5,7.5,10,12.5,15,17.5,20,22.5,25,27.5,30]; Tim = 0
nd = len(Times)
Times_ind = np.zeros(nd)
for k in range(2,nt-1):
if (TimeR[k] >= Times[Tim]):
Times_ind[Tim] = k; Tim+=1
if Tim > nd-1: break
a = np.zeros((nx,nd)); b = np.zeros((nx,nd)); c = np.zeros((nx,nd))
for k in range(0, nd):
c[:,k] = phi_PM[:,int(Times_ind[k])]
a[:,k] = phi_SM[:,int(Times_ind[k])]
b[:,k] = phi_WM[:,int(Times_ind[k])]
Masks = np.zeros((nx,nd))
for i in range(0,nd-1):
Masks[:,i] = Mask_1D[:]
avg_a = np.sum(a*Masks,axis=0)/np.sum(Mask_1D)
avg_b = np.sum(b*Masks,axis=0)/np.sum(Mask_1D)
avg_c = np.sum(c*Masks,axis=0)/np.sum(Mask_1D)
Times = [0,0.5,1,2,3,4,5]; Tim = 0
nd = len(Times)
Times_ind = np.zeros(nd)
for k in range(2,nt-1):
if (TimeR[k] >= Times[Tim]):
Times_ind[Tim] = k; Tim+=1
if Tim > nd-1: break
avg_a_ear = np.zeros((nx,nd)); avg_b_ear = np.zeros((nx,nd)); avg_c_ear = np.zeros((nx,nd))
for k in range(0, nd):
avg_c_ear[:,k] = phi_PM[:,int(Times_ind[k])]
avg_a_ear[:,k] = phi_SM[:,int(Times_ind[k])]
avg_b_ear[:,k] = phi_WM[:,int(Times_ind[k])]
a = avg_a#avvg_a
b = avg_b#avg_b
c = avg_c#avg_c
# translate the data to cartesian corrds
with np.errstate(divide='ignore', invalid='ignore'):
x_deg = 0.5 * ( 2.*b+c ) / ( a+b+c )
y_deg = 0.5*np.sqrt(3) * c / (a+b+c)
#T = tri.Triangulation(x,y)
Masks = np.zeros((nx,nd))
for i in range(0,nd-1):
Masks[:,i] = Mask_1D[:]
a_ear = np.sum(avg_a_ear*Masks,axis=0)/np.sum(Mask_1D)
b_ear = np.sum(avg_b_ear*Masks,axis=0)/np.sum(Mask_1D)
c_ear = np.sum(avg_c_ear*Masks,axis=0)/np.sum(Mask_1D)
# translate the data to cartesian corrds
with np.errstate(divide='ignore', invalid='ignore'):
x_ear = 0.5 * ( 2.*b_ear + c_ear ) / ( a_ear+b_ear+c_ear )
y_ear = 0.5*np.sqrt(3) * c_ear / (a_ear+b_ear+ c_ear)
#T = tri.Triangulation(x_ear,y_ear)
TerData1 = np.array([[0.215,0.325],[0.182,0.301],[0.165,0.206],[0.130,0.182],[0.0936,0.0862],[0.0878,0.0567],[0.0822,0.0634],[0.0776,0.0487],[0.064,0.022],[0.0662,0.00726],[0.08202,0.01389],[0.0764,0.0166]])
TerData2 = np.array([[0.150,0.227],[0.131,0.20],[0.119,0.157],[0.0993,0.118],[0.0891,0.10],[0.078,0.077],[0.0641,0.0475],[0.0584,0.0354],[0.0482,0.0154],[0.0493,0.00734]])
TerData3 = np.array([[0.043,0.0144],[0.0483,0.0287],[0.0499,0.0449],[0.0599,0.06018],[0.066,0.0754],[0.0737,0.0754],[0.09374,0.124]])
ter4=np.array([[0.281720061, 0.017675118, 0.700604821],[0.316630745, 0.021686924, 0.661682331],[0.353163198, 0.027320502, 0.6195163],[0.386452111, 0.029710538, 0.583837351],[0.421362794, 0.033722345, 0.544914861],[0.456273478, 0.037734152, 0.505992371],[0.492588655, 0.040341464, 0.467069881],[0.520042314, 0.045323213, 0.434634472],[0.553331227, 0.04771325, 0.398955523],[0.585215646, 0.05150778, 0.363276574],[0.611047534, 0.054867759, 0.334084706],[0.648984483, 0.059096842, 0.291918675],[0.674816371, 0.062456821, 0.262726808],[0.705296295, 0.067655846, 0.227047859],[0.749503049, 0.075345746, 0.175151205],[0.783009238, 0.080762047, 0.136228715],[0.816949979, 0.092230878, 0.090819143],[0.839972879, 0.098399845, 0.061627276],[0.85251249, 0.105321479, 0.042166031],[0.859434124, 0.11786109, 0.022704786],[0.860303228, 0.129966149, 0.009730623]])
tl = len(ter4[:,0])
TerData4 = np.zeros((tl,2))
TerData4[:,0] = 0.5 * ( 2.*ter4[:,1]+ter4[:,2] ) / ( ter4[:,0]+ter4[:,1]+ter4[:,2] )
TerData4[:,1] = 0.5*np.sqrt(3) * ter4[:,2] / ( ter4[:,0]+ter4[:,1]+ter4[:,2] )
TitraMV = np.copy(TitraM)
TitraMV[0,:] = 0.6*TitraM[0,:]
WperV = np.linspace(0,0.6,1000)
TitraFV = np.zeros(len(WperV))
TitraFV[:] = 1 + (0-1)/(1+np.power(WperV[:]/Titr_fits[1,Ti],Titr_fits[2,Ti]))
for i in range(0,len(WperV)-1):
if TitraFV[i] >= 0.5:
break
Meta = WperV[i]
PolVal = np.zeros(1000)
WatVal = np.zeros(1000)
SolVal = np.zeros(1000)
PolVal = np.linspace(0,1,1000)
WatVal[:] = (1-PolVal)*Meta
SolVal = (1-PolVal) * (1-Meta)
test_data = np.zeros((1000,3))
test_data[:,:] = np.transpose(np.array([SolVal,WatVal,PolVal]) )
# barycentric coords: (a,b,c)
a=test_data[:,0]
b=test_data[:,1]
c=test_data[:,2]
# translate the data to cartesian corrds
xmeta = 0.5 * ( 2.*b+c ) / ( a+b+c )
ymeta = 0.5*np.sqrt(3) * c / (a+b+c)
def PlotTer():
# =============================================================================
#
# =============================================================================
plt.figure(6,figsize=(15,6))
plt.subplot(121)
# =============================================================================
# plt.tricontourf(x,y,T.triangles,v,cmap='Greys')
# =============================================================================
# create the grid
corners = np.array([[0, 0], [1, 0], [0.5, np.sqrt(3)*0.5]])
triangle = tri.Triangulation(corners[:, 0], corners[:, 1])
# creating the grid
refiner = tri.UniformTriRefiner(triangle)
trimesh = refiner.refine_triangulation(subdiv=3)
plt.triplot(trimesh,linestyle='-',color='lightgrey')
plt.title('Polymer',fontsize=14)
plt.xlabel('NMP H2O',fontsize=14)
plt.show()
if MW == 15000:
plt.plot(x_ear,y_ear,marker='o',linestyle='--',color='lightblue')
elif MW == 29000:
plt.plot(x_ear,y_ear,marker='o',linestyle='--',color='b')
elif MW == 53000:
plt.plot(x_ear,y_ear,marker='o',linestyle='--',color='darkblue')
##plt.plot(x,y,linestyle='--',color='b')
if MW == 53000:
plt.plot(xmeta,ymeta,color='k',linestyle='--',linewidth=2,label='53 kDa')
elif MW == 29000:
plt.plot(xmeta,ymeta,color='k',linestyle='--',linewidth=2,label='29 kDa')
elif MW == 15000:
plt.plot(xmeta,ymeta,color='k',linestyle='--',linewidth=2,label='15 kDa')
plt.subplot(122)
# =============================================================================
# plt.tricontourf(x,y,T.triangles,v,cmap='Greys')
# =============================================================================
# create the grid
corners = np.array([[0, 0], [1, 0], [0.5, np.sqrt(3)*0.5]])
triangle = tri.Triangulation(corners[:, 0], corners[:, 1])
# creating the grid
refiner = tri.UniformTriRefiner(triangle)
trimesh = refiner.refine_triangulation(subdiv=3)
plt.triplot(trimesh,linestyle='-',color='lightgrey')
plt.title('Polymer',fontsize=14)
plt.xlabel('NMP H2O',fontsize=14)
plt.show()
if MW == 15000:
plt.plot(x_deg,y_deg,marker='o',linestyle='--',color='lightblue')
elif MW == 29000:
plt.plot(x_deg,y_deg,marker='o',linestyle='--',color='b')
elif MW == 53000:
plt.plot(x_deg,y_deg,marker='o',linestyle='--',color='darkblue')
#plt.plot(x,y,linestyle='--',color='b')
if MW == 53000:
plt.plot(xmeta,ymeta,color='k',linestyle='--',linewidth=2,label='53 kDa')
elif MW == 29000:
plt.plot(xmeta,ymeta,color='k',linestyle='--',linewidth=2,label='29 kDa')
elif MW == 15000:
plt.plot(xmeta,ymeta,color='k',linestyle='--',linewidth=2,label='15 kDa')
corners = np.array([[0, 0], [1, 0], [0.5, np.sqrt(3)*0.5]])
triangle = tri.Triangulation(corners[:, 0], corners[:, 1])
refiner = tri.UniformTriRefiner(triangle)
trimesh = refiner.refine_triangulation(subdiv=3)
plt.triplot(trimesh,linestyle='-',color='lightgrey')
plt.title('Polymer')
plt.xlabel('NMP H2O',fontsize=14)
plt.show()