Permalink
Cannot retrieve contributors at this time
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_10um.py
Go to fileThis commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
2114 lines (1706 sloc)
94.9 KB
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| import numpy as np | |
| import matplotlib.pyplot as plt | |
| import scipy.sparse | |
| import scipy.sparse.linalg | |
| from scipy import sparse | |
| from scipy.ndimage import gaussian_filter | |
| from scipy.integrate import odeint | |
| import time | |
| # ============================================================================= | |
| # User Defined Parameters | |
| # ============================================================================= | |
| MW = 6500 # MW of polymer | |
| Opt = 1 # Set 0 for default, Opt 1 utilizes the fitted diffusivity of drug through polymer | |
| Degr = 1 # Set Degr = 0 for no degradation | |
| pH_bath = 7.4 | |
| wvP = 0.05 # % w/v | |
| wvD = 0.1 | |
| # Units: | |
| # days, kg, m | |
| Lt = 25#37 # Duration of Simulation | |
| Ltb = 3/24 # Time in Bath | |
| Ltl = 2 # Time Lyopho. (likely unecessary) | |
| C = 0.01 # Condition Number | |
| # ============================================================================= | |
| # Degradation Constants | |
| # ============================================================================= | |
| ke1 = 0.5*0.00394 | |
| ke2 = 0 | |
| kr1 = 6.645e-05 | |
| kr2 = 0.07171875 | |
| # ============================================================================= | |
| # Acid Dissociation Constants | |
| # ============================================================================= | |
| 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) | |
| # ============================================================================= | |
| # Material Properties | |
| # ============================================================================= | |
| # density of drug, polymer, water, and solvent | |
| rho_d = 1500 | |
| rho_p = 1340#250 | |
| rho_w = 997 | |
| rho_s = 1330 | |
| # Molecular weight of drug, polymer, water, solvent, oligomers and monomers | |
| Mw_w = 18 | |
| Mw_s = 84.93 | |
| Mw_d = 331.348 | |
| Mw_col = (90.08 + 58.04 - 1.008*2)/2 * 4 | |
| Mw_cm = (90.08 + 58.04 - 1.008*2)/2 | |
| # Misc. | |
| H2O_conc = 55.5*1000 | |
| CH0 = 1000*(10**(-pH_bath)) | |
| Col0 = 0 | |
| Cm0 = 0 | |
| M0 = 72 | |
| n = 0.5 | |
| m = 4 | |
| mm = MW/(Mw_cm) | |
| Mn0 = MW | |
| # ============================================================================= | |
| # Molecular weight dependent properties based on work by Pack et. al | |
| # ============================================================================= | |
| if MW == 6500: | |
| #print('6500 Da PLGA') | |
| rb = 10/2/1000/1000 | |
| encap = 0.436 | |
| Ti = 3 | |
| Data = np.array([[0.012239902,0.048959608,0.122399021,0.29375765,0.673194614,0.966952264,1.701346389,2.141982864,2.778457772,3.720930233,4.712362301,5.740514076,6.658506732,7.906976744,8.751529988,9.70624235] | |
| ,[0.251886792,0.314779874,0.385534591,0.596226415,0.76918239,0.789622642,0.85408805,0.915408805,0.979874214,0.970440252,0.973584906,0.978301887,0.989308176,0.994025157,0.997169811,1.000314465] | |
| ,[0,0,0.031446541,0.055031447,0.012578616,0.020440252,0.020440252,0.023584906,0.028301887,0.023584906,0.025157233,0.022012579,0.020440252,0.018867925,0.020440252,0.018867925]]) | |
| if MW == 21000: | |
| #print('21000 Da PLGA') | |
| rb = 9.8/2/1000/1000 | |
| encap = 0.484 | |
| Ti = 3 | |
| Data = np.array([[0.159118727,0.342717258,0.856793146,1.150550796,1.407588739,1.921664627,2.141982864,2.839657283,3.280293758,4.051407589,4.895960832,5.850673195,6.915544676,8.017135863,8.861689106,9.877600979,10.99143207,11.89718482,12.90085679,13.89228886,15.91187271,16.90330477,17.89473684,19.91432069,21.12607099,22.26438188,23.25581395,24.24724602,26.23011016,27.91921665,31.00367197] | |
| ,[0.077358491,0.140251572,0.247169811,0.262893082,0.289622642,0.355660377,0.369811321,0.421698113,0.446855346,0.495597484,0.528616352,0.572641509,0.624528302,0.660691824,0.682704403,0.71572327,0.74245283,0.756603774,0.773899371,0.791194969,0.806918239,0.836792453,0.877672956,0.916981132,0.931132075,0.940566038,0.953144654,0.961006289,0.979874214,0.989308176,1.001886792] | |
| ,[0.003144654,0.001572327,0.009433962,0.003144654,0.001572327,0.004716981,0.004716981,0.006289308,0.007861635,0.009433962,0.009433962,0.014150943,0.022012579,0.023584906,0.023584906,0.02672956,0.029874214,0.031446541,0.034591195,0.033018868,0.034591195,0.020440252,0.017295597,0.018867925,0.011006289,0.007861635,0.007861635,0.004716981,0.006289308,0.006289308,0.012578616]]) | |
| elif MW == 34000: | |
| #print('34000 Da PLGA') | |
| rb = 11.7/2/1000/1000 | |
| encap = 0.595 | |
| Ti = 3 | |
| Data = np.array([[1.701346389,2.141982864,2.76621787,3.720930233,4.712362301,5.740514076,6.658506732,7.906976744,8.751529988,9.596083231,10.69767442,11.83598531,12.64381885,13.70869033,14.81028152,15.76499388,17.71113831,19.73072215,20.72215422,21.750306,23.73317013,24.7246022,26.00979192,27.0379437,28.02937576,29.9755202,32.17870257,34.78580171] | |
| ,[0.149685535,0.18427673,0.225157233,0.278616352,0.316352201,0.35408805,0.387106918,0.426415094,0.446855346,0.470440252,0.509748428,0.547484277,0.567924528,0.611949686,0.644968553,0.659119497,0.722012579,0.761320755,0.786477987,0.857232704,0.891823899,0.910691824,0.920125786,0.929559748,0.948427673,0.968867925,0.99245283,1.003459119] | |
| ,[0.011006289,0.01572327,0.020440252,0.012578616,0.009433962,0.007861635,0.011006289,0.009433962,0.011006289,0.009433962,0.006289308,0.011006289,0.009433962,0.012578616,0.011006289,0.009433962,0.003144654,0.001572327,0.004716981,0.009433962,0.001572327,0.004716981,0.006289308,0.006289308,0.014150943,0.020440252,0.018867925,0.011006289]]) | |
| elif MW == 49000: | |
| #print('49000 Da PLGA') | |
| rb = 10.2/2/1000/1000 | |
| encap = 0.558 | |
| Ti = 5 | |
| Data = np.array([[0.856793146,1.138310894,1.407588739,1.921664627,2.141982864,2.802937576,3.280293758,4.051407589,4.895960832,5.887392901,6.915544676,8.017135863,8.8249694,9.853121175,10.95471236,11.90942472,12.90085679,13.89228886,15.8873929,16.90330477,17.89473684,19.91432069,21.08935129,22.26438188,23.25581395,24.24724602,26.23011016,27.91921665,31.00367197] | |
| ,[0.119811321,0.13081761,0.146540881,0.188993711,0.198427673,0.229874214,0.245597484,0.28490566,0.302201258,0.336792453,0.376100629,0.410691824,0.432704403,0.479874214,0.514465409,0.555345912,0.583647799,0.619811321,0.670125786,0.71572327,0.830503145,0.901257862,0.926415094,0.942138365,0.951572327,0.96572327,0.978301887,0.989308176,1.001886792] | |
| ,[0,0.011006289,0.017295597,0.014150943,0.007861635,0.011006289,0.007861635,0.009433962,0.007861635,0.009433962,0.012578616,0.014150943,0.018867925,0.012578616,0.009433962,0.01572327,0.012578616,0.014150943,0.012578616,0.025157233,0.033018868,0.020440252,0.017295597,0.01572327,0.018867925,0.01572327,0.014150943,0.004716981,0]]) | |
| elif MW == 67000: | |
| #print('67000 Da PLGA') | |
| rb = 10.1/2/1000/1000 | |
| encap = 0.66 | |
| Ti = 7 | |
| Data = np.array([[0,0.673194614,1.701346389,2.141982864,2.76621787,3.684210526,4.712362301,5.740514076,6.658506732,7.94369645,8.751529988,9.70624235,10.77111383,11.83598531,12.68053856,13.70869033,14.81028152,15.76499388,16.75642595,17.71113831,19.73072215,20.72215422,21.750306,23.73317013,24.76132191,26.04651163,29.9755202,32.21542228,34.78580171] | |
| ,[0,0.052201258,0.099371069,0.105660377,0.107232704,0.154402516,0.18427673,0.223584906,0.244025157,0.280188679,0.300628931,0.325786164,0.363522013,0.394968553,0.421698113,0.467295597,0.498742138,0.520754717,0.54591195,0.578930818,0.629245283,0.662264151,0.832075472,0.921698113,0.948427673,0.967295597,0.983018868,0.989308176,1.000314465] | |
| ,[0,0,0,0,0,0.012578616,0.012578616,0.012578616,0.011006289,0.011006289,0.009433962,0.011006289,0.009433962,0.009433962,0.014150943,0.009433962,0.006289308,0.007861635,0.009433962,0.011006289,0.012578616,0.011006289,0.017295597,0.020440252,0.020440252,0.020440252,0.01572327,0.017295597,0.018867925]]) | |
| if Degr == 0: ke1 = ke2 = kr1 = kr2 = 0 | |
| Lx = 2*rb | |
| nx = 101 | |
| dx = Lx/nx | |
| dt = dx**2 / (1e-6) | |
| nttt = 11 | |
| ntt = 51 | |
| kon = 100#0000 # doesn't matter switch to steady state | |
| if Degr==0:kon=0 | |
| 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 | |
| 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_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) | |
| 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) | |
| Dpore_H2O = (2e-9)*(24*3600) | |
| Dpolymer_H2O = (3e-10)*(24*3600) | |
| Dpore_CH = Dpore_H2O | |
| Dpolymer_CH = Dpolymer_H2O | |
| if Opt == 1: | |
| if MW == 6500: | |
| Dpolymer_Drug =(70e-18)*60*60*24 | |
| elif MW == 21000: | |
| Dpolymer_Drug = (3.75e-18)*60*60*24 #(6e-18)*60*60*24 | |
| elif MW == 34000: | |
| Dpolymer_Drug = (3.0e-18)*60*60*24 #(6.5e-18)*60*60*24 | |
| elif MW == 49000: | |
| Dpolymer_Drug = (2.35e-18)*60*60*24 | |
| elif MW == 67000: | |
| Dpolymer_Drug = (1.5e-18)*60*60*24 | |
| fDrug = 0 | |
| rp = 0 | |
| rs = 0 | |
| # ============================================================================= | |
| # Crystallinity | |
| # ============================================================================= | |
| Vc = 4.19e-24 | |
| p = 0.01 | |
| etaA = 6.022e23 | |
| lam = 1 | |
| Xmax = 0.665 | |
| Xc0 =0#.1*Xmax | |
| 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]]) | |
| 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 TFunc60_116(phi_W,MW): | |
| return (abs(116-MW/1000)/(116-60) * (1+(0-1)/(1+np.power(phi_W[:]/Titr_fits[1,7],Titr_fits[2,7]))) | |
| + abs(60-MW/1000)/(116-60) * (1+(0-1)/(1+np.power(phi_W[:]/Titr_fits[1,9],Titr_fits[2,9])))) | |
| 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 PlotW(): | |
| plt.figure(14) | |
| plt.plot(TitraM[0,:],TitraM[Ti,:],'ok',label='Experimental Data') | |
| plt.plot(WperV/0.6,TFunc(WperV)) | |
| Nt = 2 | |
| LT = Ltb + Ltl + Lt | |
| Xext_1D = np.zeros((nx,Nt)) | |
| Xc_1D = np.zeros((nx,Nt)) | |
| Xcend_1D = np.zeros((nx,Nt)) | |
| # ============================================================================= | |
| # Initialization (poor code replace with = []) | |
| # ============================================================================= | |
| 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)) | |
| pinchColTotal_1D = np.zeros((nx,ntt*Nt)) | |
| pinchCHTotal_1D = np.zeros((nx,ntt*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)) | |
| DegWeight = np.zeros((nx,Nt)) | |
| phi_SM = np.zeros((nx,Nt)) | |
| phi_PM = np.zeros((nx,Nt)) | |
| phi_WM = np.zeros((nx,Nt)) | |
| # ============================================================================= | |
| # Time/Space | |
| # ============================================================================= | |
| x = np.linspace(0,Lx,nx) | |
| 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,DCH2O_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,DegWeight,DrugRelease_1D,TimeR,TimeT | |
| ,phi_SM,phi_WM,phi_PM) | |
| # ============================================================================= | |
| # Model | |
| # ============================================================================= | |
| 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,DCH2O_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,DegWeight,DrugRelease_1D,TimeR,TimeT | |
| ,phi_SM,phi_WM,phi_PM) = 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) | |
| # ============================================================================= | |
| # 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 | |
| )= Initialize( | |
| Mask_1D,invMask_1D,Maskfit_1D,invMaskfit_1D,invtisMask_1D,tisMask_1D,CSol_1D,CH2O_1D) | |
| # ============================================================================= | |
| # Solidification of Implant through Phase Inversion in the Bath | |
| # ============================================================================= | |
| t = 1 | |
| c = 1 | |
| while TimeT[t] <= Ltb: | |
| if TimeT[t] > c: | |
| print(f'Days = {TimeT[t]: 0.3f}') | |
| c += 1 | |
| dt = 1* dx**2 / (np.sum(Mask_1D*DDrug_1D[:,t])/np.sum(Mask_1D)) | |
| # if dt <= 1e-7: dt = 1e-7 | |
| TimeBold = Ltb * (np.linspace(0,1.1,1001)**2) | |
| # TimeBold = Ltb * (np.linspace(0,1.1,10001)**2) | |
| dt = TimeBold[t]-TimeBold[t-1] | |
| TimeT = np.append(TimeT,0) | |
| TimeT[t+1] = TimeT[t] + dt | |
| phi_SM = np.append(phi_SM,np.zeros((nx,1)),axis=1) | |
| phi_WM = np.append(phi_WM,np.zeros((nx,1)),axis=1) | |
| phi_PM = np.append(phi_PM,np.zeros((nx,1)),axis=1) | |
| DrugRelease_1D = np.append(DrugRelease_1D,np.zeros((nx,1)),axis=1) | |
| Ce_1D = np.append(Ce_1D,np.zeros((nx,1)),axis=1) | |
| Cend_1D = np.append(Cend_1D,np.zeros((nx,1)),axis=1) | |
| Res_1D = np.append(Res_1D,np.zeros((nx,1)),axis=1) | |
| Rrs_1D = np.append(Rrs_1D,np.zeros((nx,1)),axis=1) | |
| Rol_1D = np.append(Rol_1D,np.zeros((nx,1)),axis=1) | |
| Rm_1D = np.append(Rm_1D,np.zeros((nx,1)),axis=1) | |
| Cm_1D = np.append(Cm_1D,np.zeros((nx,1)),axis=1) | |
| Cmcoo_1D = np.append(Cmcoo_1D,np.zeros((nx,1)),axis=1) | |
| Col_1D = np.append(Col_1D,np.zeros((nx,1)),axis=1) | |
| Colcoo_1D = np.append(Colcoo_1D,np.zeros((nx,1)),axis=1) | |
| Cpol_1D = np.append(Cpol_1D,np.zeros((nx,1)),axis=1) | |
| Cpolcoo_1D = np.append(Cpolcoo_1D,np.zeros((nx,1)),axis=1) | |
| Drug_1D = np.append(Drug_1D,np.zeros((nx,1)),axis=1) | |
| Drugcoo_1D = np.append(Drugcoo_1D,np.zeros((nx,1)),axis=1) | |
| CH_1D = np.append(CH_1D,np.zeros((nx,1)),axis=1) | |
| DDrug_1D = np.append(DDrug_1D,np.zeros((nx,1)),axis=1) | |
| DCSol_1D = np.append(DCSol_1D,np.zeros((nx,1)),axis=1) | |
| DCm_1D = np.append(DCm_1D,np.zeros((nx,1)),axis=1) | |
| DCol_1D = np.append(DCol_1D,np.zeros((nx,1)),axis=1) | |
| DCH2O_1D = np.append(DCH2O_1D,np.zeros((nx,1)),axis=1) | |
| DCH_1D = np.append(DCH_1D,np.zeros((nx,1)),axis=1) | |
| VporeCol_1D = np.append(VporeCol_1D,np.zeros((nx,1)),axis=1) | |
| VporeDrug_1D = np.append(VporeDrug_1D,np.zeros((nx,1)),axis=1) | |
| Vpore_1D = np.append(Vpore_1D,np.zeros((nx,1)),axis=1) | |
| Xc_1D = np.append(Xc_1D,np.zeros((nx,1)),axis=1) | |
| Xcend_1D = np.append(Xcend_1D,np.zeros((nx,1)),axis=1) | |
| Xext_1D = np.append(Xext_1D,np.zeros((nx,1)),axis=1) | |
| Maskfit_1D = np.append(Maskfit_1D,np.zeros((nx,1)),axis=1) | |
| invMaskfit_1D = np.append(invMaskfit_1D,np.zeros((nx,1)),axis=1) | |
| CH2O_1D = np.append(CH2O_1D,np.zeros((nx,1)),axis=1) | |
| COH_1D = np.append(COH_1D,np.zeros((nx,1)),axis=1) | |
| CSol_1D = np.append(CSol_1D,np.zeros((nx,1)),axis=1) | |
| DegWeight = np.append(DegWeight,np.zeros((nx,1)),axis=1) | |
| Mn_1D = np.append(Mn_1D,np.zeros((nx,1)),axis=1) | |
| (Maskfit_1D,invMaskfit_1D,CSol_1D,CH2O_1D,DCSol_1D,DCH2O_1D,phi_SM,phi_PM,phi_WM | |
| ) = Solidification( | |
| Maskfit_1D,invMaskfit_1D,CSol_1D,CH2O_1D,DCSol_1D,DCH2O_1D,Sol_conc,t,dt ,rb,dx ,Vpore_1D,phi_SM,phi_PM,phi_WM,Ce0,invtisMask_1D ,Mn_1D,Cend_1D ) | |
| # ============================================================================= | |
| # Formation in Bath, Deg-Diff | |
| # ============================================================================= | |
| # 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,DegWeight | |
| ) = 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,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 | |
| ) = 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 ) | |
| # 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 ,Xc_1D,Mask_1D,invMask_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 ) | |
| #Outputs | |
| DrugRelease_1D[:,t] =invtisMask_1D[:] * (Drug_1D[:,t] +Drugcoo_1D[:,t]) | |
| t += 1 | |
| # ============================================================================= | |
| # Lyopholize, deg-diff | |
| # ============================================================================= | |
| CH2O_1D[:,t-1] = 0*tisMask_1D[:] + CH2O_1D[:,t-1] * invtisMask_1D[:] | |
| CH2O_1D[:,t] = 0*tisMask_1D[:] + CH2O_1D[:,t] * invtisMask_1D[:] | |
| while TimeT[t] <= (Ltb+Ltl): | |
| if TimeT[t] > c: | |
| print(f'Days = {TimeT[t]: 0.3f}') | |
| c += 1 | |
| dt = 1* dx**2 / Dpolymer_Drug#(np.sum(Mask_1D*DDrug_1D[:,t])/np.sum(Mask_1D)) | |
| TimeT = np.append(TimeT,0) | |
| TimeT[t+1] = TimeT[t] + dt | |
| phi_SM = np.append(phi_SM,np.zeros((nx,1)),axis=1) | |
| phi_WM = np.append(phi_WM,np.zeros((nx,1)),axis=1) | |
| phi_PM = np.append(phi_PM,np.zeros((nx,1)),axis=1) | |
| DrugRelease_1D = np.append(DrugRelease_1D,np.zeros((nx,1)),axis=1) | |
| Ce_1D = np.append(Ce_1D,np.zeros((nx,1)),axis=1) | |
| Cend_1D = np.append(Cend_1D,np.zeros((nx,1)),axis=1) | |
| Res_1D = np.append(Res_1D,np.zeros((nx,1)),axis=1) | |
| Rrs_1D = np.append(Rrs_1D,np.zeros((nx,1)),axis=1) | |
| Rol_1D = np.append(Rol_1D,np.zeros((nx,1)),axis=1) | |
| Rm_1D = np.append(Rm_1D,np.zeros((nx,1)),axis=1) | |
| Col_1D = np.append(Col_1D,np.zeros((nx,1)),axis=1) | |
| Colcoo_1D = np.append(Colcoo_1D,np.zeros((nx,1)),axis=1) | |
| Cm_1D = np.append(Cm_1D,np.zeros((nx,1)),axis=1) | |
| Cmcoo_1D = np.append(Cmcoo_1D,np.zeros((nx,1)),axis=1) | |
| Cpol_1D = np.append(Cpol_1D,np.zeros((nx,1)),axis=1) | |
| Cpolcoo_1D = np.append(Cpolcoo_1D,np.zeros((nx,1)),axis=1) | |
| Drug_1D = np.append(Drug_1D,np.zeros((nx,1)),axis=1) | |
| Drugcoo_1D = np.append(Drugcoo_1D,np.zeros((nx,1)),axis=1) | |
| CH_1D = np.append(CH_1D,np.zeros((nx,1)),axis=1) | |
| DDrug_1D = np.append(DDrug_1D,np.zeros((nx,1)),axis=1) | |
| DCSol_1D = np.append(DCSol_1D,np.zeros((nx,1)),axis=1) | |
| DCm_1D = np.append(DCm_1D,np.zeros((nx,1)),axis=1) | |
| DCol_1D = np.append(DCol_1D,np.zeros((nx,1)),axis=1) | |
| DCH2O_1D = np.append(DCH2O_1D,np.zeros((nx,1)),axis=1) | |
| DCH_1D = np.append(DCH_1D,np.zeros((nx,1)),axis=1) | |
| VporeCol_1D = np.append(VporeCol_1D,np.zeros((nx,1)),axis=1) | |
| VporeDrug_1D = np.append(VporeDrug_1D,np.zeros((nx,1)),axis=1) | |
| Vpore_1D = np.append(Vpore_1D,np.zeros((nx,1)),axis=1) | |
| Xc_1D = np.append(Xc_1D,np.zeros((nx,1)),axis=1) | |
| Xcend_1D = np.append(Xcend_1D,np.zeros((nx,1)),axis=1) | |
| Xext_1D = np.append(Xext_1D,np.zeros((nx,1)),axis=1) | |
| Maskfit_1D = np.append(Maskfit_1D,np.zeros((nx,1)),axis=1) | |
| invMaskfit_1D = np.append(invMaskfit_1D,np.zeros((nx,1)),axis=1) | |
| CH2O_1D = np.append(CH2O_1D,np.zeros((nx,1)),axis=1) | |
| COH_1D = np.append(COH_1D,np.zeros((nx,1)),axis=1) | |
| CSol_1D = np.append(CSol_1D,np.zeros((nx,1)),axis=1) | |
| DegWeight = np.append(DegWeight,np.zeros((nx,1)),axis=1) | |
| Mn_1D = np.append(Mn_1D,np.zeros((nx,1)),axis=1) | |
| (Maskfit_1D,invMaskfit_1D,CSol_1D,CH2O_1D,DCSol_1D,DCH2O_1D,phi_SM,phi_PM,phi_WM | |
| ) = Solidification( | |
| Maskfit_1D,invMaskfit_1D,CSol_1D,CH2O_1D,DCSol_1D,DCH2O_1D,Sol_conc,t,dt ,rb,dx ,Vpore_1D,phi_SM,phi_PM,phi_WM,Ce0,invtisMask_1D ,Mn_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,DegWeight | |
| ) = 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,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 | |
| ) = LYODiffusion( | |
| 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,DCH2O_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 ,Xc_1D,Mask_1D,invMask_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 | |
| ) = LYOUpdateDiffusion( | |
| 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,DCH2O_1D ) | |
| t+=1 | |
| # ============================================================================= | |
| # Begin Time Iteration | |
| # ============================================================================= | |
| # ============================================================================= | |
| nstart = t | |
| CH2O_1D[:,t] = H2O_conc*tisMask_1D[:] + 0*invtisMask_1D[:] | |
| CH2O_1D[:,t-1] = H2O_conc*tisMask_1D[:] + 0*invtisMask_1D[:] | |
| while TimeT[t] <= (LT): | |
| if TimeT[t] > c: | |
| print(f'Days = {TimeT[t]: 0.3f}') | |
| c += 1 | |
| drugleft = np.sum(Mask_1D*(Drug_1D[:,t-1]+Drugcoo_1D[:,t-1]))/np.sum(Mask_1D*(Drug_1D[:,nstart]+Drugcoo_1D[:,nstart])) | |
| dt = C * dx**2 / (np.sum(Mask_1D*DDrug_1D[:,t])/np.sum(Mask_1D)) | |
| # dt = 0.1* x**2 / (np.sum(Mask_1D*DDrug_1D[:,t])/np.sum(Mask_1D)) | |
| if drugleft <= 0.01: dt = 0.01 | |
| if dt <= 1e-2: dt = 1e-2 #Converged | |
| # dt = 1e-3 | |
| TimeT = np.append(TimeT,0) | |
| TimeT[t+1] = TimeT[t] + dt | |
| TimeR = np.append(TimeR,0) | |
| TimeR[t-nstart+1] = TimeR[t-nstart] + dt | |
| phi_SM = np.append(phi_SM,np.zeros((nx,1)),axis=1) | |
| phi_WM = np.append(phi_WM,np.zeros((nx,1)),axis=1) | |
| phi_PM = np.append(phi_PM,np.zeros((nx,1)),axis=1) | |
| DrugRelease_1D = np.append(DrugRelease_1D,np.zeros((nx,1)),axis=1) | |
| Ce_1D = np.append(Ce_1D,np.zeros((nx,1)),axis=1) | |
| Cend_1D = np.append(Cend_1D,np.zeros((nx,1)),axis=1) | |
| Res_1D = np.append(Res_1D,np.zeros((nx,1)),axis=1) | |
| Rrs_1D = np.append(Rrs_1D,np.zeros((nx,1)),axis=1) | |
| Rol_1D = np.append(Rol_1D,np.zeros((nx,1)),axis=1) | |
| Rm_1D = np.append(Rm_1D,np.zeros((nx,1)),axis=1) | |
| Col_1D = np.append(Col_1D,np.zeros((nx,1)),axis=1) | |
| Colcoo_1D = np.append(Colcoo_1D,np.zeros((nx,1)),axis=1) | |
| Cm_1D = np.append(Cm_1D,np.zeros((nx,1)),axis=1) | |
| Cmcoo_1D = np.append(Cmcoo_1D,np.zeros((nx,1)),axis=1) | |
| Cpol_1D = np.append(Cpol_1D,np.zeros((nx,1)),axis=1) | |
| Cpolcoo_1D = np.append(Cpolcoo_1D,np.zeros((nx,1)),axis=1) | |
| Drug_1D = np.append(Drug_1D,np.zeros((nx,1)),axis=1) | |
| Drugcoo_1D = np.append(Drugcoo_1D,np.zeros((nx,1)),axis=1) | |
| CH_1D = np.append(CH_1D,np.zeros((nx,1)),axis=1) | |
| DDrug_1D = np.append(DDrug_1D,np.zeros((nx,1)),axis=1) | |
| DCSol_1D = np.append(DCSol_1D,np.zeros((nx,1)),axis=1) | |
| DCm_1D = np.append(DCm_1D,np.zeros((nx,1)),axis=1) | |
| DCol_1D = np.append(DCol_1D,np.zeros((nx,1)),axis=1) | |
| DCH2O_1D = np.append(DCH2O_1D,np.zeros((nx,1)),axis=1) | |
| DCH_1D = np.append(DCH_1D,np.zeros((nx,1)),axis=1) | |
| VporeCol_1D = np.append(VporeCol_1D,np.zeros((nx,1)),axis=1) | |
| VporeDrug_1D = np.append(VporeDrug_1D,np.zeros((nx,1)),axis=1) | |
| Vpore_1D = np.append(Vpore_1D,np.zeros((nx,1)),axis=1) | |
| Xc_1D = np.append(Xc_1D,np.zeros((nx,1)),axis=1) | |
| Xcend_1D = np.append(Xcend_1D,np.zeros((nx,1)),axis=1) | |
| Xext_1D = np.append(Xext_1D,np.zeros((nx,1)),axis=1) | |
| Maskfit_1D = np.append(Maskfit_1D,np.zeros((nx,1)),axis=1) | |
| invMaskfit_1D = np.append(invMaskfit_1D,np.zeros((nx,1)),axis=1) | |
| CH2O_1D = np.append(CH2O_1D,np.zeros((nx,1)),axis=1) | |
| COH_1D = np.append(COH_1D,np.zeros((nx,1)),axis=1) | |
| CSol_1D = np.append(CSol_1D,np.zeros((nx,1)),axis=1) | |
| DegWeight = np.append(DegWeight,np.zeros((nx,1)),axis=1) | |
| Mn_1D = np.append(Mn_1D,np.zeros((nx,1)),axis=1) | |
| (Maskfit_1D,invMaskfit_1D,CSol_1D,CH2O_1D,DCSol_1D,DCH2O_1D,phi_SM,phi_PM,phi_WM | |
| ) = Solidification( | |
| Maskfit_1D,invMaskfit_1D,CSol_1D,CH2O_1D,DCSol_1D,DCH2O_1D,Sol_conc,t,dt ,rb,dx ,Vpore_1D,phi_SM,phi_PM,phi_WM,Ce0,invtisMask_1D ,Mn_1D,Cend_1D ) | |
| Maskfit_1D[:,t+1] = Mask_1D[:] | |
| invMaskfit_1D[:,t+1] = invMask_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,DegWeight | |
| ) = 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,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 | |
| ) = 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 ) | |
| # 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,Xc_1D,Mask_1D,invMask_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 ) | |
| #Outputs | |
| DrugRelease_1D[:,t] =invtisMask_1D[:] * (Drug_1D[:,t] +Drugcoo_1D[:,t]) | |
| t+=1 | |
| nend = t | |
| # ============================================================================= | |
| # ============================================================================= | |
| # Output Variables | |
| # ============================================================================= | |
| SumDrug = np.sum(DrugRelease_1D*Mw_d*dx ,axis = 0) | |
| SumMn = np.sum(Mn_1D,axis = 0) | |
| SumXc = np.sum(Xc_1D,axis = 0) | |
| ProfileMn = np.zeros(t+1) | |
| ProfileDrug = np.zeros(t+1) | |
| ProfileXc = np.zeros(t+1) | |
| ProfileMn[:-nstart] = (SumMn[(nstart):])/SumMn[nstart] | |
| ProfileDrug[:-nstart] =100* (SumDrug[nstart:])/SumDrug[nstart] | |
| ProfileXc[:-nstart] = (SumXc[nstart:])/SumXc[nstart] | |
| return(ProfileDrug,ProfileMn,ProfileXc,DDrug_1D,Drug_1D,Drugcoo_1D,Maskfit_1D,CSol_1D,CH2O_1D,TimeR,TimeT,nstart, | |
| Vpore_1D,Col_1D,Colcoo_1D,DCol_1D,Ce_1D,Xc_1D,Rol_1D,Rm_1D,Mn_1D,nend,CH_1D) | |
| # ============================================================================= | |
| # ============================================================================= | |
| # ============================================================================= | |
| # # # User Defined Functions | |
| # ============================================================================= | |
| # ============================================================================= | |
| # ============================================================================= | |
| def Maskfunc( | |
| rb,rp,rs,x,Lx,Mask_1D,invMask_1D,tisMask_1D,invtisMask_1D): | |
| #Removes everything except core region, where interior pores exist | |
| tempc =((x[0:]-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 | |
| # Mask_1D[0:int(2.5/2/1000/1000/rb*nx)] = 0 | |
| # Mask_1D[int(15/2/1000/1000/rb*nx):int(18/2/1000/1000/rb*nx)] = 0 | |
| M = (Mask_1D[:] < 1) | |
| invMask_1D = M*1 | |
| tis = ((x[:]-0*0.5*Lx)**2) > (rb)**2 | |
| tisMask_1D[:] = tis*1 | |
| Mm = (tisMask_1D[:] < 1) | |
| invtisMask_1D = Mm*1 | |
| 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): | |
| Vimplant = (4/3*3.141592)*rb**3 | |
| Vsolution = 2e-5 # arbitray | |
| mp = wvP *(Vsolution*1e3) | |
| ms = Vsolution*rho_s | |
| vp = mp/rho_p | |
| md = wvD * mp | |
| vd = md/rho_d | |
| Vtot = Vsolution+vp+vd | |
| cp = (mp*1000)/MW/Vtot; cs = (ms*1000)/Mw_s/Vtot; cd = (md*1000)/Mw_d/Vtot | |
| phi_poly = MW*cp/1000/rho_p | |
| phi_sol = Mw_s*cs/1000/rho_s | |
| phi_drug = Mw_d*cd/1000/rho_d | |
| ############################################################################## | |
| Ce0 = (1-Xc0)* mm*cp | |
| 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 = cs | |
| phi_P =invtisMask_1D* phi_poly | |
| phi_S = invtisMask_1D* phi_sol | |
| 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 | |
| ############################################################################## | |
| Drug_1D[:,0] = cd * invtisMask_1D[:] | |
| Drug_1D[:,1] = cd * invtisMask_1D[:] | |
| Ce_1D[:,0] = Ce0*Mask_1D[:] | |
| Ce_1D[:,1] = Ce0 *Mask_1D[:] | |
| Cend_1D[:,0] = 2*Cpol0*Mask_1D[:] | |
| Cend_1D[:,1] = 2*Cpol0*Mask_1D[:] | |
| 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 | |
| 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] = 0 | |
| COH_1D[-1,0] = 0 | |
| CH2O_1D[-1,0] = H2O_conc | |
| 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] = CH0 | |
| COH_1D[-1,1] = 0 | |
| CH2O_1D[-1,1] = H2O_conc | |
| CSol_1D[-1,1] = 0 | |
| Cpol_1D[-1,1] = 0 | |
| Cpolcoo_1D[-1,1] = 0 | |
| Drug_1D[0,0] = 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 | |
| CH2O_1D[0,0] = H2O_conc | |
| CSol_1D[0,0] = 0 | |
| Cpol_1D[0,0] = 0 | |
| Cpolcoo_1D[0,0] = 0 | |
| Drug_1D[0,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 | |
| CH2O_1D[0,1] = H2O_conc | |
| 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) | |
| def Solidification(Maskfit_1D,invMaskfit_1D,CSol_1D,CH2O_1D,DCSol_1D,DCH2O_1D,Sol_conc,t,dt ,rb,dx,Vpore_1D,phi_SM,phi_PM,phi_WM ,Ce0,invtisMask_1D,Mn_1D,Cend_1D): | |
| r = np.linspace(0,Lx,nx) | |
| CH2O_1D[0,t] = CH2O_1D[1,t] | |
| CSol_1D[0,t] = CSol_1D[1,t] | |
| phi_SM[:,t] = invtisMask_1D*(Mw_s*CSol_1D[:,t] /1000/rho_s) | |
| phi_PM[:,t] =Mask_1D* (MW *0.5*Cend_1D[:,t-1] /1000/rho_p) | |
| phi_WM[:,t] = (1 - phi_SM[:,t] - phi_PM[:,t]) | |
| phi_WM[phi_WM<=0] = 0 | |
| # 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))/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 | |
| 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 = invMaskfit_1D[:,t]*1/D12Bar_theo + Maskfit_1D[:,t]*(Dpolymer_Sol+(1.3*Vpore_1D[:,t-1]**2-0.3*Vpore_1D[:,t-1]**3)*(1/D12Bar_theo-Dpolymer_Sol)) | |
| D11bar = invMaskfit_1D[:,t]*1/D11Bar_theo + Maskfit_1D[:,t]*(Dpolymer_Sol+(1.3*Vpore_1D[:,t-1]**2-0.3*Vpore_1D[:,t-1]**3)*(1/D11Bar_theo-Dpolymer_Sol)) | |
| D22bar = invMaskfit_1D[:,t]*1/D22Bar_theo + Maskfit_1D[:,t]*(Dpolymer_Sol+(1.3*Vpore_1D[:,t-1]**2-0.3*Vpore_1D[:,t-1]**3)*(1/D22Bar_theo-Dpolymer_Sol)) | |
| D12bar = invMaskfit_1D[:,t]*1/D12Bar_theo + Maskfit_1D[:,t+1]*(1/D12Bar_poly_theo + (1.3*Vpore_1D[:,t]**2-0.3*Vpore_1D[:,t]**3)*(1/D12Bar_theo - 1/D12Bar_poly_theo))#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]/I[0:-2] | |
| +(D12bar[2:nx]+D12bar[1:nx-1])/2 ) | |
| sig_xLC[0:nx-2] = -(w1)*(dt/(dr)**2)*( | |
| -D12bar[1:nx-1]/I[0:-2] | |
| +(D12bar[1:nx-1]+D12bar[0:nx-2])/2 ) | |
| sig_xyC[1:nx-1] = 1 + (w1)*(dt/(dr)**2)*( | |
| +(D12bar[2:nx]+D12bar[1:nx-1])/2 | |
| +(D12bar[1:nx-1]+D12bar[0:nx-2])/2 | |
| ) | |
| sig_xyC[0] = 1 +dt*D12bar[0]*6*dt/dr**2 | |
| sig_xRC[1] = -dt*D12bar[0]*6*dt/dr**2 | |
| sig_xyC[-1] = 1#-dr*(1-1/r[-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]/I[0:-2] | |
| +(D12bar[2:nx]+D12bar[1:nx-1])/2 ) | |
| sig_xLC[0:nx-2] = -(w1)*(dt/(dr)**2)*( | |
| -D12bar[1:nx-1]/I[0:-2] | |
| +(D12bar[1:nx-1]+D12bar[0:nx-2])/2 ) | |
| sig_xyC[1:nx-1] = 1 + (w1)*(dt/(dr)**2)*( | |
| +(D12bar[2:nx]+D12bar[1:nx-1])/2 | |
| +(D12bar[1:nx-1]+D12bar[0:nx-2])/2 | |
| ) | |
| sig_xyC[0] = 1 +dt*D12bar[0]*6*dt/dr**2 | |
| sig_xRC[1] = -dt*D12bar[0]*6*dt/dr**2 | |
| sig_xyC[-1] = 1#-dr*(1-1/r[-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) | |
| # Solve System of Equations Governing Diffusion | |
| bh2o = (B5.dot(CH2O_1D[:,t]) ) | |
| ch2o = (C5.dot(CH2O_1D[:,t-1]) ) | |
| Bh2o = bh2o-ch2o | |
| bsol = (B6.dot(CSol_1D[:,t] ) ) | |
| csol = (C6.dot(CSol_1D[:,t-1] ) ) | |
| Bsol = bsol-csol | |
| temph = scipy.sparse.linalg.spsolve(A5, Bh2o) | |
| temps = scipy.sparse.linalg.spsolve(A6, Bsol) | |
| CH2O_1D[:,t+1] = temph[:] | |
| CSol_1D[:,t+1] = temps[:] | |
| # ============================================================================= | |
| # Solidification | |
| # ============================================================================= | |
| Maskfit_1D[:,t+1] = TFunc(phi_WM[:,t]) * Mask_1D[:] | |
| if MW <= 33000: Maskfit_1D[:,t+1] = TFunc8_33(phi_WM[:,t],MW) * Mask_1D[:] | |
| elif ((MW>33000)&(MW<=50000)): Maskfit_1D[:,t+1] = TFunc33_50(phi_WM[:,t],MW) * Mask_1D[:] | |
| elif ((MW>50000)&(MW<=60000)): Maskfit_1D[:,t+1] = TFunc50_60(phi_WM[:,t],MW) * Mask_1D[:] | |
| elif MW > 60000: Maskfit_1D[:,t+1] = TFunc60_116(phi_WM[:,t],MW) * Mask_1D[:] | |
| Maskfit_1D[(Maskfit_1D[:,t+1])>=0.999,t+1] = 1 | |
| Maskfit_1D[(Maskfit_1D[:,t+1])>=1,t+1] = 1 | |
| invMaskfit_1D[:,t+1] = 1 - Maskfit_1D[:,t+1] | |
| # ============================================================================= | |
| # Update Diffusion | |
| # ============================================================================= | |
| Dpolymer_Sol_Swell = Dpolymer_Sol | |
| D12bar = invMaskfit_1D[:,t+1]*1/D12Bar_theo + Maskfit_1D[:,t+1]*(Dpolymer_Sol_Swell+(1.3*Vpore_1D[:,t]**2-0.3*Vpore_1D[:,t]**3)*(1/D12Bar_theo-Dpolymer_Sol_Swell)) | |
| D11bar = invMaskfit_1D[:,t+1]*1/D11Bar_theo + Maskfit_1D[:,t+1]*(Dpolymer_Sol+(1.3*Vpore_1D[:,t]**2-0.3*Vpore_1D[:,t]**3)*(1/D11Bar_theo-Dpolymer_Sol)) | |
| D22bar = invMaskfit_1D[:,t+1]*1/D22Bar_theo + Maskfit_1D[:,t+1]*(Dpolymer_Sol+(1.3*Vpore_1D[:,t]**2-0.3*Vpore_1D[:,t]**3)*(1/D22Bar_theo-Dpolymer_Sol)) | |
| D12bar = invMaskfit_1D[:,t+1]*1/D12Bar_theo + Maskfit_1D[:,t+1]*(1/D12Bar_poly_theo + (1.3*Vpore_1D[:,t]**2-0.3*Vpore_1D[:,t]**3)*(1/D12Bar_theo - 1/D12Bar_poly_theo))#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]/I[0:-2] | |
| +(D12bar[2:nx]+D12bar[1:nx-1])/2 ) | |
| sig_xLC[0:nx-2] = -(w1)*(dt/(dr)**2)*( | |
| -D12bar[1:nx-1]/I[0:-2] | |
| +(D12bar[1:nx-1]+D12bar[0:nx-2])/2 ) | |
| sig_xyC[1:nx-1] = 1 + (w1)*(dt/(dr)**2)*( | |
| +(D12bar[2:nx]+D12bar[1:nx-1])/2 | |
| +(D12bar[1:nx-1]+D12bar[0:nx-2])/2 | |
| ) | |
| sig_xyC[0] = 1 +dt*D12bar[0]*6*dt/dr**2 | |
| sig_xRC[1] = -dt*D12bar[0]*6*dt/dr**2 | |
| sig_xyC[-1] = 1#-dr*(1-1/r[-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]/I[0:-2] | |
| +(D12bar[2:nx]+D12bar[1:nx-1])/2 ) | |
| sig_xLC[0:nx-2] = -(w1)*(dt/(dr)**2)*( | |
| -D12bar[1:nx-1]/I[0:-2] | |
| +(D12bar[1:nx-1]+D12bar[0:nx-2])/2 ) | |
| sig_xyC[1:nx-1] = 1 + (w1)*(dt/(dr)**2)*( | |
| +(D12bar[2:nx]+D12bar[1:nx-1])/2 | |
| +(D12bar[1:nx-1]+D12bar[0:nx-2])/2 | |
| ) | |
| sig_xyC[0] = 1 +dt*D12bar[0]*6*dt/dr**2 | |
| sig_xRC[1] = -dt*D12bar[0]*6*dt/dr**2 | |
| sig_xyC[-1] = 1#-dr*(1-1/r[-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* (scipy.sparse.linalg.spsolve(A5, Bh2o)- phi_PM[:,t]) | |
| temps = rho_s/Mw_s*1000* scipy.sparse.linalg.spsolve(A6, Bsol) | |
| CH2O_1D[:,t+1] = temph[:] | |
| CSol_1D[:,t+1] = temps[:] | |
| CH2O_1D[-1,t+1] = CH2O_1D[-1,0] | |
| return Maskfit_1D,invMaskfit_1D,CSol_1D,CH2O_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,cd,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 = (tf - t0) / float(ntt) | |
| h = (dt) / float(ntt) | |
| col[:,0] = Col0 | |
| colcoo[:,0] = Colcoo0 | |
| cm[:,0] = Cm0 | |
| cmcoo[:,0] = Cmcoo0 | |
| cpol[:,0] = Cpol0 | |
| cpolcoo[:,0] = Cpolcoo0 | |
| drug[:,0] = cd | |
| drugcoo[:,0] = Drugcoo0 | |
| ch[:,0] = CH0 | |
| # ut[0] = t0 | |
| 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] = t0 + (i+1)*h | |
| 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] | |
| # ttimeR,pinchCol_1D,pinchColcoo_1D,pinchCm_1D,pinchCmcoo_1D,pinchCpol_1D,pinchCpolcoo_1D,pinchDrug_1D,pinchDrugcoo_1D,pinchCH_1D = RK(fCol,fColcoo,fCm,fCmcoo,fCpol,fCpolcoo,fDrug,fDrugcoo,fCH,pinchCol0,pinchColcoo0,pinchCm0,pinchCmcoo0,pinchCpol0,pinchCpolcoo0,pinchDrug0,pinchDrugcoo0,pinchCH0,dt,ntt,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,CH_1D ): | |
| # ============================================================================= | |
| # Recalculate the rate of scission for new time step based on C_H+ | |
| # ============================================================================= | |
| Xext_1D[:,t] = Mask_1D[:]*(p*etaA*(Rrs_1D[:,t]+Res_1D[:,t])*Vc) | |
| if lam == 1: | |
| Xc_1D[:,t] = Mask_1D[:]*(Xmax - (Xmax - Xc0)*np.exp(-Xext_1D[:,t])) | |
| else: | |
| Xc_1D[:,t] = Mask_1D[:]*(Xmax - ((lam-1)*Xext_1D[:,t] + (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])) | |
| Cend_1D[:,t] = Mask_1D[:]* 2*( Cpol0 + (Rrs_1D[:,t] - Rol_1D[:,t]/m)) | |
| 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] = 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] = Mask_1D[index] * ((Ce_1D[index,t]+omega*Xc_1D[index,t]) * M0)/(0.5*Cend_1D[index,t] + invMask_1D[index]) | |
| 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] | |
| DegWeight[:,t+1] = (CH2O_1D[:,t+1]/(np.max(CH2O_1D[:,t+1]))) | |
| for ttt in range(0,nttt-1): | |
| dttt = dt/nttt | |
| Rrs_temp[:,ttt+1] = Rrs_temp[:,ttt] + DegWeight[:,t+1]*(dttt*(Ce_1D[:,t] * (kr1 + kr2*CH_1D[:,t]))) | |
| Res_temp[:,ttt+1] = Res_temp[:,ttt] + DegWeight[:,t+1]*(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] | |
| 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) | |
| 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 ): | |
| sig_xRC = np.zeros(nx) | |
| sig_xLC = np.zeros(nx) | |
| sig_xyC = np.zeros(nx) | |
| r = np.linspace(0,Lx,nx)#/Lx | |
| dr = r[1]-r[0] | |
| I = np.linspace(1,nx,nx) | |
| # ============================================================================= | |
| # Diffusion | |
| # ============================================================================= | |
| w1 = 1#2/3 | |
| w2 = 1#4/3 | |
| w3 = 0#1/3 | |
| sig_xRC[2:nx] = -(w1)*(dt/(dr)**2)* ( | |
| DCol_1D[1:nx-1,t]/I[0:-2] | |
| +(DCol_1D[2:nx,t]+DCol_1D[1:nx-1,t])/2 ) | |
| sig_xLC[0:nx-2] = -(w1)*(dt/(dr)**2)*( | |
| -DCol_1D[1:nx-1,t]/I[0:-2] | |
| +(DCol_1D[1:nx-1,t]+DCol_1D[0:nx-2,t])/2 ) | |
| sig_xyC[1:nx-1] = 1 + (w1)*(dt/(dr)**2)*( | |
| +(DCol_1D[2:nx,t]+DCol_1D[1:nx-1,t])/2 | |
| +(DCol_1D[1:nx-1,t]+DCol_1D[0:nx-2,t])/2 | |
| ) | |
| sig_xyC[0] = 1 +dt*DCol_1D[0,t]*6*dt/dr**2 | |
| sig_xRC[1] = -dt*DCol_1D[0,t]*6*dt/dr**2 | |
| sig_xyC[-1] = 1#-dr*(1-1/r[-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)* ( | |
| DDrug_1D[1:nx-1,t]/I[0:-2] | |
| +(DDrug_1D[2:nx,t]+DDrug_1D[1:nx-1,t])/2 ) | |
| sig_xLC[0:nx-2] = -(w1)*(dt/(dr)**2)*( | |
| -DDrug_1D[1:nx-1,t]/I[0:-2] | |
| +(DDrug_1D[1:nx-1,t]+DDrug_1D[0:nx-2,t])/2 ) | |
| sig_xyC[1:nx-1] = 1 + (w1)*(dt/(dr)**2)*( | |
| +(DDrug_1D[2:nx,t]+DDrug_1D[1:nx-1,t])/2 | |
| +(DDrug_1D[1:nx-1,t]+DDrug_1D[0:nx-2,t])/2 | |
| ) | |
| sig_xyC[0] = 1 +dt*DDrug_1D[0,t]*6*dt/dr**2 | |
| sig_xRC[1] = -dt*DDrug_1D[0,t]*6*dt/dr**2 | |
| sig_xyC[-1] = 1#-dr*(1-1/r[-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)* ( | |
| DCH_1D[1:nx-1,t]/I[0:-2] | |
| +(DCH_1D[2:nx,t]+DCH_1D[1:nx-1,t])/2 ) | |
| sig_xLC[0:nx-2] = -(w1)*(dt/(dr)**2)*( | |
| -DCH_1D[1:nx-1,t]/I[0:-2] | |
| +(DCH_1D[1:nx-1,t]+DCH_1D[0:nx-2,t])/2 ) | |
| sig_xyC[1:nx-1] = 1 + (w1)*(dt/(dr)**2)*( | |
| +(DCH_1D[2:nx,t]+DCH_1D[1:nx-1,t])/2 | |
| +(DCH_1D[1:nx-1,t]+DCH_1D[0:nx-2,t])/2 | |
| ) | |
| sig_xyC[0] = 1 +dt*DCH_1D[0,t]*6*dt/dr**2 | |
| sig_xRC[1] = -dt*DCH_1D[0,t]*6*dt/dr**2 | |
| sig_xyC[-1] = 1#-dr*(1-1/r[-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 = np.zeros(nx) | |
| sig_xLC = np.zeros(nx) | |
| sig_xyC = np.zeros(nx) | |
| sig_xRC[2:nx] = -(w1)*(dt/(dr)**2)* ( | |
| DCm_1D[1:nx-1,t]/I[0:-2] | |
| +(DCm_1D[2:nx,t]+DCm_1D[1:nx-1,t])/2 ) | |
| sig_xLC[0:nx-2] = -(w1)*(dt/(dr)**2)*( | |
| -DCm_1D[1:nx-1,t]/I[0:-2] | |
| +(DCm_1D[1:nx-1,t]+DCm_1D[0:nx-2,t])/2 ) | |
| sig_xyC[1:nx-1] = 1 + (w1)*(dt/(dr)**2)*( | |
| +(DCm_1D[2:nx,t]+DCm_1D[1:nx-1,t])/2 | |
| +(DCm_1D[1:nx-1,t]+DCm_1D[0:nx-2,t])/2 | |
| ) | |
| sig_xyC[0] = 1 +dt*DCm_1D[0,t]*6*dt/dr**2 | |
| sig_xRC[1] = -dt*DCm_1D[0,t]*6*dt/dr**2 | |
| sig_xyC[-1] = 1#-dr*(1-1/r[-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 | |
| # ============================================================================= | |
| Drug_1D[:,t+1] = scipy.sparse.linalg.spsolve(A2, Bdrug) | |
| Col_1D[:,t+1] = scipy.sparse.linalg.spsolve(A1, Bcol) | |
| Cm_1D[:,t+1] = scipy.sparse.linalg.spsolve(A4, Bcm) | |
| Drugcoo_1D[:,t+1] = scipy.sparse.linalg.spsolve(A2, Bdrugcoo) | |
| Colcoo_1D[:,t+1] = scipy.sparse.linalg.spsolve(A1, Bcolcoo) | |
| Cmcoo_1D[:,t+1] = scipy.sparse.linalg.spsolve(A4, Bcmcoo) | |
| CH_1D[:,t+1] = scipy.sparse.linalg.spsolve(A3, Bch) | |
| CH2O_1D[:,t+1] = scipy.sparse.linalg.spsolve(A3, Bch2o) | |
| COH_1D[:,t+1] = scipy.sparse.linalg.spsolve(A3, Bcoh) | |
| Drug_1D[0,t+1] = Drug_1D[1,t+1] | |
| 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] | |
| # CH2O_1D[:,t+1] = scipy.sparse.linalg.spsolve(A3, Bch2o) | |
| COH_1D[0,t+1] = COH_1D[1,t+1] | |
| 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 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,Xc_1D,Mask_1D,invMask_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+1] | |
| ) | |
| with np.errstate(divide='ignore', invalid='ignore'): | |
| 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]) | |
| DCol_1D[:,t+1] = invMaskfit_1D[:,t+1]*Dpore_Col + Maskfit_1D[:,t+1]*(Dpolymer_Col+(1.3*Vpore_1D[:,t+1]**2-0.3*Vpore_1D[:,t+1]**3)*(Dpore_Col-Dpolymer_Col)) | |
| DCm_1D[:,t+1] = invMaskfit_1D[:,t+1]*Dpore_Cm + Maskfit_1D[:,t+1]*(Dpolymer_Cm+(1.3*Vpore_1D[:,t+1]**2-0.3*Vpore_1D[:,t+1]**3)*(Dpore_Cm-Dpolymer_Cm)) | |
| DDrug_1D[:,t+1] = invMaskfit_1D[:,t+1]*Dpore_Drug + Maskfit_1D[:,t+1]*(Dpolymer_Drug+(1.3*Vpore_1D[:,t+1]**2-0.3*Vpore_1D[:,t+1]**3)*(Dpore_Drug-Dpolymer_Drug)) | |
| DCH_1D[:,t+1] = invMaskfit_1D[:,t+1]*Dpore_CH + Maskfit_1D[:,t+1]*(Dpolymer_CH+(1.3*Vpore_1D[:,t+1]**2-0.3*Vpore_1D[:,t+1]**3)*(Dpore_CH-Dpolymer_CH)) | |
| 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 ): | |
| r = np.linspace(dx,Lx+dx,nx)#/Lx | |
| dr = r[1]-r[0] | |
| I = np.linspace(1,nx,nx) | |
| sig_xRC = np.zeros(nx) | |
| sig_xLC = np.zeros(nx) | |
| sig_xyC = np.zeros(nx) | |
| # ============================================================================= | |
| # Diffusion | |
| # ============================================================================= | |
| w1 = 1#2/3 | |
| w2 = 1#4/3 | |
| w3 = 0#1/3 | |
| sig_xRC[2:nx] = -(w1)*(dt/(dr)**2)* ( | |
| DCol_1D[1:nx-1,t+1]/I[0:-2] | |
| +(DCol_1D[2:nx,t+1]+DCol_1D[1:nx-1,t+1])/2 ) | |
| sig_xLC[0:nx-2] = -(w1)*(dt/(dr)**2)*( | |
| -DCol_1D[1:nx-1,t+1]/I[0:-2] | |
| +(DCol_1D[1:nx-1,t+1]+DCol_1D[0:nx-2,t+1])/2 ) | |
| sig_xyC[1:nx-1] = 1 + (w1)*(dt/(dr)**2)*( | |
| +(DCol_1D[2:nx,t+1]+DCol_1D[1:nx-1,t+1])/2 | |
| +(DCol_1D[1:nx-1,t+1]+DCol_1D[0:nx-2,t+1])/2 | |
| ) | |
| sig_xyC[0] = 1 +dt*DCol_1D[0,t+1]*6*dt/dr**2 | |
| sig_xRC[1] = -dt*DCol_1D[0,t+1]*6*dt/dr**2 | |
| sig_xyC[-1] = 1#-dr*(1-1/r[-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)* ( | |
| DDrug_1D[1:nx-1,t+1]/I[0:-2] | |
| +(DDrug_1D[2:nx,t+1]+DDrug_1D[1:nx-1,t+1])/2 ) | |
| sig_xLC[0:nx-2] = -(w1)*(dt/(dr)**2)*( | |
| -DDrug_1D[1:nx-1,t+1]/I[0:-2] | |
| +(DDrug_1D[1:nx-1,t+1]+DDrug_1D[0:nx-2,t+1])/2 ) | |
| sig_xyC[1:nx-1] = 1 + (w1)*(dt/(dr)**2)*( | |
| +(DDrug_1D[2:nx,t+1]+DDrug_1D[1:nx-1,t+1])/2 | |
| +(DDrug_1D[1:nx-1,t+1]+DDrug_1D[0:nx-2,t+1])/2 | |
| ) | |
| sig_xyC[0] = 1 +dt*DDrug_1D[0,t+1]*6*dt/dr**2 | |
| sig_xRC[1] = -dt*DDrug_1D[0,t+1]*6*dt/dr**2 | |
| sig_xyC[-1] = 1#-dr*(1-1/r[-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)* ( | |
| DCH_1D[1:nx-1,t+1]/I[0:-2] | |
| +(DCH_1D[2:nx,t+1]+DCH_1D[1:nx-1,t+1])/2 ) | |
| sig_xLC[0:nx-2] = -(w1)*(dt/(dr)**2)*( | |
| -DCH_1D[1:nx-1,t+1]/I[0:-2] | |
| +(DCH_1D[1:nx-1,t+1]+DCH_1D[0:nx-2,t+1])/2 ) | |
| sig_xyC[1:nx-1] = 1 + (w1)*(dt/(dr)**2)*( | |
| +(DCH_1D[2:nx,t+1]+DCH_1D[1:nx-1,t+1])/2 | |
| +(DCH_1D[1:nx-1,t+1]+DCH_1D[0:nx-2,t+1])/2 | |
| ) | |
| sig_xyC[0] = 1 +dt*DCH_1D[0,t+1]*6*dt/dr**2 | |
| sig_xRC[1] = -dt*DCH_1D[0,t+1]*6*dt/dr**2 | |
| sig_xyC[-1] = 1#-dr*(1-1/r[-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 = np.zeros(nx) | |
| sig_xLC = np.zeros(nx) | |
| sig_xyC = np.zeros(nx) | |
| sig_xRC[2:nx] = -(w1)*(dt/(dr)**2)* ( | |
| DCm_1D[1:nx-1,t+1]/I[0:-2] | |
| +(DCm_1D[2:nx,t+1]+DCm_1D[1:nx-1,t+1])/2 ) | |
| sig_xLC[0:nx-2] = -(w1)*(dt/(dr)**2)*( | |
| -DCm_1D[1:nx-1,t+1]/I[0:-2] | |
| +(DCm_1D[1:nx-1,t+1]+DCm_1D[0:nx-2,t+1])/2 ) | |
| sig_xyC[1:nx-1] = 1 + (w1)*(dt/(dr)**2)*( | |
| +(DCm_1D[2:nx,t+1]+DCm_1D[1:nx-1,t+1])/2 | |
| +(DCm_1D[1:nx-1,t+1]+DCm_1D[0:nx-2,t+1])/2 | |
| ) | |
| sig_xyC[0] = 1 +dt*DCm_1D[0,t+1]*6*dt/dr**2 | |
| sig_xRC[1] = -dt*DCm_1D[0,t+1]*6*dt/dr**2 | |
| sig_xyC[-1] = 1#-dr*(1-1/r[-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 | |
| # ============================================================================= | |
| Drug_1D[:,t+1] = scipy.sparse.linalg.spsolve(A2, Bdrug) | |
| Col_1D[:,t+1] = scipy.sparse.linalg.spsolve(A1, Bcol) | |
| Cm_1D[:,t+1] = scipy.sparse.linalg.spsolve(A4, Bcm) | |
| Drugcoo_1D[:,t+1] = scipy.sparse.linalg.spsolve(A2, Bdrugcoo) | |
| Colcoo_1D[:,t+1] = scipy.sparse.linalg.spsolve(A1, Bcolcoo) | |
| Cmcoo_1D[:,t+1] = scipy.sparse.linalg.spsolve(A4, Bcmcoo) | |
| CH_1D[:,t+1] = scipy.sparse.linalg.spsolve(A3, Bch) | |
| CH2O_1D[:,t+1] = scipy.sparse.linalg.spsolve(A3, Bch2o) | |
| COH_1D[:,t+1] = scipy.sparse.linalg.spsolve(A3, Bcoh) | |
| Drug_1D[0,t+1] = Drug_1D[1,t+1] | |
| 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] | |
| # CH2O_1D[:,t+1] = scipy.sparse.linalg.spsolve(A3, Bch2o) | |
| COH_1D[0,t+1] = COH_1D[1,t+1] | |
| 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 LYODiffusion( | |
| 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,DCH2O_1D ): | |
| sig_xRC = np.zeros(nx) | |
| sig_xLC = np.zeros(nx) | |
| sig_xyC = np.zeros(nx) | |
| # ============================================================================= | |
| # Diffusion | |
| # ============================================================================= | |
| w1 = 1#2/3 | |
| w2 = 1#4/3 | |
| w3 = 0#1/3 | |
| sig_xRC = np.zeros(nx) | |
| sig_xLC = np.zeros(nx) | |
| sig_xyC = np.zeros(nx) | |
| sig_xRC[2:int(rb/Lx*nx+1)] = -(w1)*(dt/(dx)**2)*(DCol_1D[2:int(rb/Lx*nx+1),t] + DCol_1D[1:int(rb/Lx*nx+1)-1,t])/2 #* ((x[2:int(rb/Lx*nx+1)]+x[1:int(rb/Lx*nx+1)-1])/2)**2/(x[1:int(rb/Lx*nx+1)-1]**2) | |
| sig_xLC[0:int(rb/Lx*nx-1)] = -(w1)*(dt/(dx)**2)*(DCol_1D[1:int(rb/Lx*nx-1)+1,t] + DCol_1D[0:int(rb/Lx*nx-1),t])/2 #* ((x[0:int(rb/Lx*nx-1)]+x[1:int(rb/Lx*nx-1)+1])/2)**2/(x[1:int(rb/Lx*nx-1)+1]**2) | |
| sig_xLC[0] = 0 | |
| sig_xLC[int(rb/Lx*nx-1)] = -(w1)*(dt/(dx)**2)*(DCol_1D[int(rb/Lx*nx-1)+1,t] + DCol_1D[int(rb/Lx*nx-1),t])/2 #* ((x[int(rb/Lx*nx-1)]+x[int(rb/Lx*nx-1)+1])/2)**2/(x[int(rb/Lx*nx-1)+1]**2) | |
| sig_xyC[1:int(rb/Lx*nx)] = 1 - sig_xRC[2:int(rb/Lx*nx+1)] - sig_xLC[0:int(rb/Lx*nx-1)] | |
| sig_xyC[0] = 1 | |
| sig_xyC[int(rb/Lx*nx)] = 1 - sig_xLC[int(rb/Lx*nx)-1] | |
| sig_xyC[int(rb/Lx*nx+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 = np.zeros(nx) | |
| sig_xLC = np.zeros(nx) | |
| sig_xyC = np.zeros(nx) | |
| sig_xRC[2:int(rb/Lx*nx+1)] = -(w1)*(dt/(dx)**2)*(DDrug_1D[2:int(rb/Lx*nx+1),t] + DDrug_1D[1:int(rb/Lx*nx+1)-1,t])/2 #* ((x[2:int(rb/Lx*nx+1)]+x[1:int(rb/Lx*nx+1)-1])/2)**2/(x[1:int(rb/Lx*nx+1)-1]**2) | |
| sig_xLC[0:int(rb/Lx*nx-1)] = -(w1)*(dt/(dx)**2)*(DDrug_1D[1:int(rb/Lx*nx-1)+1,t] + DDrug_1D[0:int(rb/Lx*nx-1),t])/2 #* ((x[0:int(rb/Lx*nx-1)]+x[1:int(rb/Lx*nx-1)+1])/2)**2/(x[1:int(rb/Lx*nx-1)+1]**2) | |
| sig_xLC[0] = 0 | |
| sig_xLC[int(rb/Lx*nx-1)] = -(w1)*(dt/(dx)**2)*(DDrug_1D[int(rb/Lx*nx-1)+1,t] + DDrug_1D[int(rb/Lx*nx-1),t])/2 #* ((x[int(rb/Lx*nx-1)]+x[int(rb/Lx*nx-1)+1])/2)**2/(x[int(rb/Lx*nx-1)+1]**2) | |
| sig_xyC[1:int(rb/Lx*nx)] = 1 - sig_xRC[2:int(rb/Lx*nx+1)] - sig_xLC[0:int(rb/Lx*nx-1)] | |
| sig_xyC[0] = 1 | |
| sig_xyC[int(rb/Lx*nx)] = 1 - sig_xLC[int(rb/Lx*nx)-1] | |
| sig_xyC[int(rb/Lx*nx+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 = np.zeros(nx) | |
| sig_xLC = np.zeros(nx) | |
| sig_xyC = np.zeros(nx) | |
| sig_xRC[2:int(rb/Lx*nx+1)] = -(w1)*(dt/(dx)**2)*(DCH_1D[2:int(rb/Lx*nx+1),t] + DCH_1D[1:int(rb/Lx*nx+1)-1,t])/2 #* ((x[2:int(rb/Lx*nx+1)]+x[1:int(rb/Lx*nx+1)-1])/2)**2/(x[1:int(rb/Lx*nx+1)-1]**2) | |
| sig_xLC[0:int(rb/Lx*nx-1)] = -(w1)*(dt/(dx)**2)*(DCH_1D[1:int(rb/Lx*nx-1)+1,t] + DCH_1D[0:int(rb/Lx*nx-1),t])/2 #* ((x[0:int(rb/Lx*nx-1)]+x[1:int(rb/Lx*nx-1)+1])/2)**2/(x[1:int(rb/Lx*nx-1)+1]**2) | |
| sig_xLC[0] = 0 | |
| sig_xLC[int(rb/Lx*nx-1)] = -(w1)*(dt/(dx)**2)*(DCH_1D[int(rb/Lx*nx-1)+1,t] + DCH_1D[int(rb/Lx*nx-1),t])/2 #* ((x[int(rb/Lx*nx-1)]+x[int(rb/Lx*nx-1)+1])/2)**2/(x[int(rb/Lx*nx-1)+1]**2) | |
| sig_xyC[1:int(rb/Lx*nx)] = 1 - sig_xRC[2:int(rb/Lx*nx+1)] - sig_xLC[0:int(rb/Lx*nx-1)] | |
| sig_xyC[0] = 1 | |
| sig_xyC[int(rb/Lx*nx)] = 1 - sig_xLC[int(rb/Lx*nx)-1] | |
| sig_xyC[int(rb/Lx*nx+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 = np.zeros(nx) | |
| sig_xLC = np.zeros(nx) | |
| sig_xyC = np.zeros(nx) | |
| sig_xRC[2:int(rb/Lx*nx+1)] = -(w1)*(dt/(dx)**2)*(DCm_1D[2:int(rb/Lx*nx+1),t] + DCm_1D[1:int(rb/Lx*nx+1)-1,t])/2 #* ((x[2:int(rb/Lx*nx+1)]+x[1:int(rb/Lx*nx+1)-1])/2)**2/(x[1:int(rb/Lx*nx+1)-1]**2) | |
| sig_xLC[0:int(rb/Lx*nx-1)] = -(w1)*(dt/(dx)**2)*(DCm_1D[1:int(rb/Lx*nx-1)+1,t] + DCm_1D[0:int(rb/Lx*nx-1),t])/2 #* ((x[0:int(rb/Lx*nx-1)]+x[1:int(rb/Lx*nx-1)+1])/2)**2/(x[1:int(rb/Lx*nx-1)+1]**2) | |
| sig_xLC[0] = 0 | |
| sig_xLC[int(rb/Lx*nx-1)] = -(w1)*(dt/(dx)**2)*(DCm_1D[int(rb/Lx*nx-1)+1,t] + DCm_1D[int(rb/Lx*nx-1),t])/2 #* ((x[int(rb/Lx*nx-1)]+x[int(rb/Lx*nx-1)+1])/2)**2/(x[int(rb/Lx*nx-1)+1]**2) | |
| sig_xyC[1:int(rb/Lx*nx)] = 1 - sig_xRC[2:int(rb/Lx*nx+1)] - sig_xLC[0:int(rb/Lx*nx-1)] | |
| sig_xyC[0] = 1 | |
| sig_xyC[int(rb/Lx*nx)] = 1 - sig_xLC[int(rb/Lx*nx)-1] | |
| sig_xyC[int(rb/Lx*nx+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) | |
| sig_xRC = np.zeros(nx) | |
| sig_xLC = np.zeros(nx) | |
| sig_xyC = np.zeros(nx) | |
| sig_xRC[2:nx] = -(w1)*(dt/(dx)**2)*(DCH2O_1D[2:nx,t] + DCH2O_1D[1:nx-1,t])/2 #* ((x[2:nx]+x[1:nx-1])/2)**2/(x[1:nx-1]**2) | |
| sig_xLC[0:nx-2] = -(w1)*(dt/(dx)**2)*(DCH2O_1D[1:nx-1,t] + DCH2O_1D[0:nx-2,t])/2 #* ((x[0:nx-2]+x[1:nx-1])/2)**2/(x[1:nx-1]**2) | |
| sig_xLC[0] = 0 | |
| sig_xyC[1:nx-1] = 1 - sig_xRC[2:nx] - sig_xLC[0:nx-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) | |
| # ============================================================================= | |
| # 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 = (B5.dot(CH2O_1D[:,t])) | |
| cch2o = (C5.dot(CH2O_1D[:,t-1])) | |
| Bch2o = bch2o - cch2o | |
| bcoh = (B3.dot(COH_1D[:,t])) | |
| ccoh = (C3.dot(COH_1D[:,t-1])) | |
| Bcoh = bcoh - ccoh | |
| # ============================================================================= | |
| Drug_1D[:,t+1] = scipy.sparse.linalg.spsolve(A2, Bdrug) | |
| Col_1D[:,t+1] = scipy.sparse.linalg.spsolve(A1, Bcol) #+ dt*(Rol_1D[:,t+1]-Rol_1D[:,t])/dt#0.35* ((a*b*(Rrs_1D[:,t+1]/Ce0)**(b-1))*(Rrs_1D[:,t+1]-Rrs_1D[:,t])/dt) | |
| Cm_1D[:,t+1] = scipy.sparse.linalg.spsolve(A4, Bcm) #+ dt* ((Rm_1D[:,t+1]-Rm_1D[:,t])/dt) | |
| Drugcoo_1D[:,t+1] = scipy.sparse.linalg.spsolve(A2, Bdrugcoo) | |
| Colcoo_1D[:,t+1] = scipy.sparse.linalg.spsolve(A1, Bcolcoo) | |
| Cmcoo_1D[:,t+1] = scipy.sparse.linalg.spsolve(A4, Bcmcoo) | |
| CH_1D[:,t+1] = scipy.sparse.linalg.spsolve(A3, Bch) | |
| CH2O_1D[:,t+1] = scipy.sparse.linalg.spsolve(A5, Bch2o) | |
| COH_1D[:,t+1] = scipy.sparse.linalg.spsolve(A3, Bcoh) | |
| 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 LYOUpdateDiffusion( | |
| 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 ,DCH2O_1D ): | |
| sig_xRC = np.zeros(nx) | |
| sig_xLC = np.zeros(nx) | |
| sig_xyC = np.zeros(nx) | |
| # ============================================================================= | |
| # Diffusion | |
| # ============================================================================= | |
| w1 = 1#2/3 | |
| w2 = 1#4/3 | |
| w3 = 0#1/3 | |
| sig_xRC = np.zeros(nx) | |
| sig_xLC = np.zeros(nx) | |
| sig_xyC = np.zeros(nx) | |
| sig_xRC[2:int(rb/Lx*nx+1)] = -(w1)*(dt/(dx)**2)*(DCol_1D[2:int(rb/Lx*nx+1),t+1] + DCol_1D[1:int(rb/Lx*nx+1)-1,t+1])/2 #* ((x[2:int(rb/Lx*nx+1)]+x[1:int(rb/Lx*nx+1)-1])/2)**2/(x[1:int(rb/Lx*nx+1)-1]**2) | |
| sig_xLC[0:int(rb/Lx*nx-1)] = -(w1)*(dt/(dx)**2)*(DCol_1D[1:int(rb/Lx*nx-1)+1,t+1] + DCol_1D[0:int(rb/Lx*nx-1),t+1])/2 #* ((x[0:int(rb/Lx*nx-1)]+x[1:int(rb/Lx*nx-1)+1])/2)**2/(x[1:int(rb/Lx*nx-1)+1]**2) | |
| sig_xLC[0] = 0 | |
| sig_xLC[int(rb/Lx*nx-1)] = -(w1)*(dt/(dx)**2)*(DCol_1D[int(rb/Lx*nx-1)+1,t+1] + DCol_1D[int(rb/Lx*nx-1),t+1])/2 #* ((x[int(rb/Lx*nx-1)]+x[int(rb/Lx*nx-1)+1])/2)**2/(x[int(rb/Lx*nx-1)+1]**2) | |
| sig_xyC[1:int(rb/Lx*nx)] = 1 - sig_xRC[2:int(rb/Lx*nx+1)] - sig_xLC[0:int(rb/Lx*nx-1)] | |
| sig_xyC[0] = 1 | |
| sig_xyC[int(rb/Lx*nx)] = 1 - sig_xLC[int(rb/Lx*nx)-1] | |
| sig_xyC[int(rb/Lx*nx+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 = np.zeros(nx) | |
| sig_xLC = np.zeros(nx) | |
| sig_xyC = np.zeros(nx) | |
| sig_xRC[2:int(rb/Lx*nx+1)] = -(w1)*(dt/(dx)**2)*(DDrug_1D[2:int(rb/Lx*nx+1),t+1] + DDrug_1D[1:int(rb/Lx*nx+1)-1,t+1])/2 #* ((x[2:int(rb/Lx*nx+1)]+x[1:int(rb/Lx*nx+1)-1])/2)**2/(x[1:int(rb/Lx*nx+1)-1]**2) | |
| sig_xLC[0:int(rb/Lx*nx-1)] = -(w1)*(dt/(dx)**2)*(DDrug_1D[1:int(rb/Lx*nx-1)+1,t+1] + DDrug_1D[0:int(rb/Lx*nx-1),t+1])/2 #* ((x[0:int(rb/Lx*nx-1)]+x[1:int(rb/Lx*nx-1)+1])/2)**2/(x[1:int(rb/Lx*nx-1)+1]**2) | |
| sig_xLC[0] = 0 | |
| sig_xLC[int(rb/Lx*nx-1)] = -(w1)*(dt/(dx)**2)*(DDrug_1D[int(rb/Lx*nx-1)+1,t+1] + DDrug_1D[int(rb/Lx*nx-1),t+1])/2 #* ((x[int(rb/Lx*nx-1)]+x[int(rb/Lx*nx-1)+1])/2)**2/(x[int(rb/Lx*nx-1)+1]**2) | |
| sig_xyC[1:int(rb/Lx*nx)] = 1 - sig_xRC[2:int(rb/Lx*nx+1)] - sig_xLC[0:int(rb/Lx*nx-1)] | |
| sig_xyC[0] = 1 | |
| sig_xyC[int(rb/Lx*nx)] = 1 - sig_xLC[int(rb/Lx*nx)-1] | |
| sig_xyC[int(rb/Lx*nx+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 = np.zeros(nx) | |
| sig_xLC = np.zeros(nx) | |
| sig_xyC = np.zeros(nx) | |
| sig_xRC[2:int(rb/Lx*nx+1)] = -(w1)*(dt/(dx)**2)*(DCH_1D[2:int(rb/Lx*nx+1),t+1] + DCH_1D[1:int(rb/Lx*nx+1)-1,t+1])/2 #* ((x[2:int(rb/Lx*nx+1)]+x[1:int(rb/Lx*nx+1)-1])/2)**2/(x[1:int(rb/Lx*nx+1)-1]**2) | |
| sig_xLC[0:int(rb/Lx*nx-1)] = -(w1)*(dt/(dx)**2)*(DCH_1D[1:int(rb/Lx*nx-1)+1,t+1] + DCH_1D[0:int(rb/Lx*nx-1),t+1])/2 #* ((x[0:int(rb/Lx*nx-1)]+x[1:int(rb/Lx*nx-1)+1])/2)**2/(x[1:int(rb/Lx*nx-1)+1]**2) | |
| sig_xLC[0] = 0 | |
| sig_xLC[int(rb/Lx*nx-1)] = -(w1)*(dt/(dx)**2)*(DCH_1D[int(rb/Lx*nx-1)+1,t+1] + DCH_1D[int(rb/Lx*nx-1),t+1])/2 #* ((x[int(rb/Lx*nx-1)]+x[int(rb/Lx*nx-1)+1])/2)**2/(x[int(rb/Lx*nx-1)+1]**2) | |
| sig_xyC[1:int(rb/Lx*nx)] = 1 - sig_xRC[2:int(rb/Lx*nx+1)] - sig_xLC[0:int(rb/Lx*nx-1)] | |
| sig_xyC[0] = 1 | |
| sig_xyC[int(rb/Lx*nx)] = 1 - sig_xLC[int(rb/Lx*nx)-1] | |
| sig_xyC[int(rb/Lx*nx+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 = np.zeros(nx) | |
| sig_xLC = np.zeros(nx) | |
| sig_xyC = np.zeros(nx) | |
| sig_xRC[2:int(rb/Lx*nx+1)] = -(w1)*(dt/(dx)**2)*(DCm_1D[2:int(rb/Lx*nx+1),t+1] + DCm_1D[1:int(rb/Lx*nx+1)-1,t+1])/2 #* ((x[2:int(rb/Lx*nx+1)]+x[1:int(rb/Lx*nx+1)-1])/2)**2/(x[1:int(rb/Lx*nx+1)-1]**2) | |
| sig_xLC[0:int(rb/Lx*nx-1)] = -(w1)*(dt/(dx)**2)*(DCm_1D[1:int(rb/Lx*nx-1)+1,t+1] + DCm_1D[0:int(rb/Lx*nx-1),t+1])/2 #* ((x[0:int(rb/Lx*nx-1)]+x[1:int(rb/Lx*nx-1)+1])/2)**2/(x[1:int(rb/Lx*nx-1)+1]**2) | |
| sig_xLC[0] = 0 | |
| sig_xLC[int(rb/Lx*nx-1)] = -(w1)*(dt/(dx)**2)*(DCm_1D[int(rb/Lx*nx-1)+1,t+1] + DCm_1D[int(rb/Lx*nx-1),t+1])/2 #* ((x[int(rb/Lx*nx-1)]+x[int(rb/Lx*nx-1)+1])/2)**2/(x[int(rb/Lx*nx-1)+1]**2) | |
| sig_xyC[1:int(rb/Lx*nx)] = 1 - sig_xRC[2:int(rb/Lx*nx+1)] - sig_xLC[0:int(rb/Lx*nx-1)] | |
| sig_xyC[0] = 1 | |
| sig_xyC[int(rb/Lx*nx)] = 1 - sig_xLC[int(rb/Lx*nx)-1] | |
| sig_xyC[int(rb/Lx*nx+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) | |
| sig_xRC = np.zeros(nx) | |
| sig_xLC = np.zeros(nx) | |
| sig_xyC = np.zeros(nx) | |
| sig_xRC[2:nx] = -(w1)*(dt/(dx)**2)*(DCH2O_1D[2:nx,t+1] + DCH2O_1D[1:nx-1,t+1])/2 #* ((x[2:nx]+x[1:nx-1])/2)**2/(x[1:nx-1]**2) | |
| sig_xLC[0:nx-2] = -(w1)*(dt/(dx)**2)*(DCH2O_1D[1:nx-1,t+1] + DCH2O_1D[0:nx-2,t+1])/2 #* ((x[0:nx-2]+x[1:nx-1])/2)**2/(x[1:nx-1]**2) | |
| sig_xLC[0] = 0 | |
| sig_xyC[1:nx-1] = 1 - sig_xRC[2:nx] - sig_xLC[0:nx-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) | |
| # ============================================================================= | |
| # 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 = (B5.dot(CH2O_1D[:,t])) | |
| cch2o = (C5.dot(CH2O_1D[:,t-1])) | |
| Bch2o = bch2o - cch2o | |
| bcoh = (B3.dot(COH_1D[:,t])) | |
| ccoh = (C3.dot(COH_1D[:,t-1])) | |
| Bcoh = bcoh - ccoh | |
| # ============================================================================= | |
| Drug_1D[:,t+1] = scipy.sparse.linalg.spsolve(A2, Bdrug) | |
| Col_1D[:,t+1] = scipy.sparse.linalg.spsolve(A1, Bcol) # + dt*(Rol_1D[:,t+1]-Rol_1D[:,t])/dt#0.35* ((a*b*(Rrs_1D[:,t+1]/Ce0)**(b-1))*(Rrs_1D[:,t+1]-Rrs_1D[:,t])/dt) | |
| Cm_1D[:,t+1] = scipy.sparse.linalg.spsolve(A4, Bcm) #+ dt* ((Rm_1D[:,t+1]-Rm_1D[:,t])/dt) | |
| Drugcoo_1D[:,t+1] = scipy.sparse.linalg.spsolve(A2, Bdrugcoo) | |
| Colcoo_1D[:,t+1] = scipy.sparse.linalg.spsolve(A1, Bcolcoo) | |
| Cmcoo_1D[:,t+1] = scipy.sparse.linalg.spsolve(A4, Bcmcoo) | |
| CH_1D[:,t+1] = scipy.sparse.linalg.spsolve(A3, Bch) | |
| CH2O_1D[:,t+1] = scipy.sparse.linalg.spsolve(A5, Bch2o) | |
| COH_1D[:,t+1] = scipy.sparse.linalg.spsolve(A3, Bcoh) | |
| 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) | |
| start = time.process_time() | |
| ProfileDrug,ProfileMn,ProfileXc,DDrug_1D,Drug_1D,Drugcoo_1D,Maskfit_1D,CSol_1D,CH2O_1D,TimeR,TimeT,nstart,Vpore_1D,Col_1D,Colcoo_1D,DCol_1D,Ce_1D,Xc_1D,Rol_1D,Rm_1D,Mn_1D,nend,CH_1D = ImplantModel(Variables) | |
| end = time.process_time() | |
| WallTime = end-start | |
| print('WallTime:', (end-start)/60, '(min)') | |
| FinalDrug = ProfileDrug[:nend] | |
| FinalMN = ProfileMn[:nend] | |
| #fig = plt.figure(23, figsize = (15,4)) | |
| #plt.subplot(121) | |
| #plt.plot(TimeR[:-2],ProfileXc[:-(ntb+ntl+2)]) | |
| ##plt.ylim((0,100)) | |
| #plt.xlim((0,60)) | |
| #plt.title('Crystallinity') | |
| #plt.subplot(122) | |
| #plt.plot(TimeR[:-2],Vpore_1D[20,:-(ntb+ntl+2)]) | |
| ##plt.ylim((0,100)) | |
| #plt.xlim((0,60)) | |
| #plt.title('Porosity') | |
| # | |
| def PlotDrugF(): | |
| plt.figure(22,figsize=(14,6)) | |
| if Degr == 1: | |
| if Opt == 1: | |
| ax1 = plt.subplot(121) | |
| # if Opt == 1: plt.plot(TimeR[:-1],(100 - ProfileDrug[:-(ntb+ntl+1)])/100,color='royalblue',linestyle='--',linewidth=2.5) | |
| if Opt == 1: plt.plot(TimeR[:-1],(100 - ProfileDrug[:-(nstart+1)])/100,color='navy',linestyle='--',linewidth=2.5) | |
| # plt.plot(TimeR[:-1],1 - ProfileDrug[:-(ntb+ntl+1)]/100,c='lightgray',linestyle='--',linewidth=2) | |
| #plt.errorbar(Data[0,:],Data[1,:],yerr=Data[2,:],color='k',marker='o') | |
| #plt.plot(Data[0,:],Data[1,:],'ok') | |
| # plt.errorbar(Data[0,:],Data[1,:],yerr=Data[2,:],color='k',marker='o',linestyle='',capsize=3,ecolor='gray') | |
| plt.xlim((-1,(Data[0,-1]+1))) | |
| plt.xlabel('Time (Days)', fontsize=20) | |
| plt.ylabel('Drug Release (normalized)',fontsize=20) | |
| ax1.xaxis.set_tick_params(labelsize=14) | |
| ax1.yaxis.set_tick_params(labelsize=14) | |
| elif Opt == 0: | |
| ax1 = plt.subplot(121) | |
| plt.plot(TimeR[:-1],1 - ProfileDrug[:-(nstart)]/100,c='navy',linewidth=2) | |
| # if Opt == 0: plt.plot(TimeR[:-1],(100 - ProfileDrug[:-(nstart)])/100)#,color='navy',linestyle='-',linewidth=2.5) | |
| #plt.errorbar(Data[0,:],Data[1,:],yerr=Data[2,:],color='k',marker='o') | |
| #plt.plot(Data[0,:],Data[1,:],'ok') | |
| # plt.errorbar(Data[0,:],Data[1,:],yerr=Data[2,:],color='k',marker='o',linestyle='',capsize=3,ecolor='gray') | |
| plt.xlim((-1,(Data[0,-1]+1))) | |
| plt.xlabel('Time (Days)', fontsize=20) | |
| plt.ylabel('Drug Release (normalized)',fontsize=20) | |
| ax1.xaxis.set_tick_params(labelsize=14) | |
| ax1.yaxis.set_tick_params(labelsize=14) | |
| elif Degr == 0: | |
| ax1 = plt.subplot(121) | |
| plt.plot(TimeR[:-1],(100 - ProfileDrug[:-(nstart)])/100,color='grey',linestyle='--',linewidth=2.5) | |
| # plt.plot(TimeR[:-1],1 - ProfileDrug[:-(ntb+ntl+1)]/100,c='lightgray',linestyle='--',linewidth=2) | |
| #plt.errorbar(Data[0,:],Data[1,:],yerr=Data[2,:],color='k',marker='o') | |
| #plt.plot(Data[0,:],Data[1,:],'ok') | |
| plt.errorbar(Data[0,:],Data[1,:],yerr=Data[2,:],color='k',marker='o',linestyle='',capsize=3,ecolor='gray') | |
| plt.xlim((-1,(Data[0,-1]+1))) | |
| plt.xlabel('Time (Days)', fontsize=20) | |
| plt.ylabel('Drug Release (normalized)',fontsize=20) | |
| ax1.xaxis.set_tick_params(labelsize=14) | |
| ax1.yaxis.set_tick_params(labelsize=14) | |
| #def PlotDrug(): | |
| # plt.figure(1,figsize=(8,6)) | |
| # | |
| # KMs = [ke1m, ke2m, kr1m, kr2m] | |
| # plt.plot(TimeR[:-2],(100 - ProfileDrug[:-(nstart+1)])/100,label='MW: %2.0f, KMs: %1.2f %1.2f %1.2f %1.2f' %(int(MW/1000),KMs[0],KMs[1],KMs[2],KMs[3])) #,color='royalblue',linestyle='--',linewidth=2.5) | |
| # plt.errorbar(Data[0,:],Data[1,:],yerr=Data[2,:],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 Release (normalized)',fontsize=20) | |
| # plt.legend(loc='best') | |
| # plt.show() | |
| #PlotDrug() | |
| def PlotDrug(): | |
| fig, ax1 = plt.subplots(1, 1, figsize=(7,6)) | |
| if Opt == 1: ax1.plot(TimeR[:-2],1 - FinalDrug[:nend-nstart]/100)#,label='MW: %2.0f, KMs: %1.2f %1.2f %1.2f %1.2f' %(int(MW/1000),KMs[0],KMs[1],KMs[2],KMs[3])) #,color='royalblue',linestyle='--',linewidth=2.5) | |
| else: ax1.plot(TimeR[:-2],1 - FinalDrug[:nend-nstart]/100)#,label='MW: %2.0f, KMs: %1.2f %1.2f %1.2f %1.2f' %(int(MW/1000),KMs[0],KMs[1],KMs[2],KMs[3])) #,color='royalblue',linestyle='--',linewidth=2.5)#)#,color='navy',linestyle='-',linewidth=2.5) | |
| ax1.errorbar(Data[0,:],Data[1,:],yerr=Data[2,:],color='k',marker='o',linestyle='',capsize=3,ecolor='gray') | |
| #plt.plot(Data[0,:],Data[1,:],'ok') | |
| ax1.set_xlim((-1,(Data[0,-1]+1))) | |
| ax1.set_ylim((0,1)) | |
| ax1.set_xlabel('Time (Days)', fontsize=20) | |
| ax1.set_ylabel('Drug Release (normalized)',fontsize=20) | |
| ax1.xaxis.set_tick_params(labelsize=14) | |
| ax1.yaxis.set_tick_params(labelsize=14) | |
| #plt.legend(loc='best') | |
| plt.show() | |
| PlotDrug() | |
| encapT = (np.sum(Drug_1D[:,nstart]+Drugcoo_1D[:,nstart],axis=0))/(np.sum(Drug_1D[:,0]+Drugcoo_1D[:,0],axis=0)) | |
| print(ke1,ke2,kr1,kr2) | |
| print('Exp:',encap) | |
| print('Theo:',encapT) | |
| # fig = plt.figure(2) | |
| # plt.plot(CH2O_1D[:,0:5]) | |
| print('\007') |