Description of BTPDE

0001 function [TOUT,YOUT,MF_cmpts,MF_allcmpts,difftime,elapsed_time] ...
0002     = BTPDE(experiment,mymesh,DIFF_cmpts,kappa_bdys,IC_cmpts)
0003 
0004 % solve Bloch-Torrey equation
0005 %
0006 % Input:
0007 %     1. experiment is a structure with 10 elements:
0008 %         ngdir_total
0009 %         gdir
0010 %         sdeltavec
0011 %         bdeltavec
0012 %         seqvec
0013 %         npervec
0014 %         rtol
0015 %         atol
0016 %         qvalues
0017 %         bvalues
0018 %     2. mymesh is a structure with 10 elements:
0019 %         Nnode
0020 %         Nele
0021 %         Nface
0022 %         Pts_cmpt_reorder
0023 %         Ele_cmpt_reorder
0024 %         Pts_ind
0025 %         Pts_boundary_reorder
0026 %         Fac_boundary_reorder
0027 %         Nboundary
0028 %         Ncmpt
0029 %     3. DIFF_cmpts
0030 %     4. kappa_bdys
0031 %     5. IC_cmpts
0032 %
0033 % Output:
0034 %     1. TOUT
0035 %     2. YOUT
0036 %     3. MF_cmpts
0037 %     4. MF_allcmpts
0038 %     5. difftime
0039 %     6. elapsed_time
0040 
0041 global FEM_M FEM_K FEM_A FEM_Q FEM_G
0042 global QVAL UG 
0043 global BDELTA SDELTA SEQ OGSEPER 
0044   
0045 bvalues = experiment.bvalues;
0046 qvalues = experiment.qvalues;
0047 gdir = experiment.gdir;
0048 sdeltavec = experiment.sdeltavec;
0049 bdeltavec = experiment.bdeltavec;
0050 seqvec = experiment.seqvec;
0051 npervec = experiment.npervec;
0052 ODEsolve_rtol = experiment.rtol;
0053 ODEsolve_atol = experiment.atol;
0054 yes = 1;  no = 0;
0055 
0056 disp(['Magnetization: Solving Bloch-Torrey PDE']);
0057 
0058 UG = gdir';
0059 UG = UG/norm(UG);
0060 
0061 if (mymesh.Ncmpt == 1 | abs(max(kappa_bdys)) <= 1e-16)
0062     DO_COUPLING = no;
0063 else
0064     DO_COUPLING = yes;
0065 end
0066 
0067 nexperi = length(sdeltavec);
0068 %disp(['Simulating ',num2str(nexperi),' experiments']);
0069 nb = size(qvalues,2);
0070 elapsed_time = zeros(nb, nexperi);
0071 
0072 for icmpt = 1:mymesh.Ncmpt
0073     %disp(['Working on cmpt ', num2str(icmpt)]);
0074     coordinates = mymesh.Pts_cmpt_reorder{icmpt};
0075     elements = mymesh.Ele_cmpt_reorder{icmpt};
0076     facets = [];
0077     GX = -sqrt(-1)*UG*coordinates;
0078     FEM_MAT{icmpt}.Q = sparse(length(coordinates),length(coordinates));
0079     for iboundary = 1:mymesh.Nboundary
0080         if (kappa_bdys(iboundary) ~= 0)
0081             neumann = mymesh.Fac_boundary_reorder{icmpt}{iboundary}';
0082             neumann_nodes = unique(neumann);
0083             %zvalues = mymesh.Pts_cmpt_reorder{1}(3,neumann_nodes);
0084             %mvalues = sqrt((20-abs(zvalues)));
0085             coeffs_flux_matrix = zeros(max(neumann_nodes),1);
0086             coeffs_flux_matrix(neumann_nodes) = kappa_bdys(iboundary); %*mvalues
0087             if ~isempty(neumann)
0088                 FEM_MAT{icmpt}.Q = FEM_MAT{icmpt}.Q + flux_matrixP1_3D(neumann,coordinates',coeffs_flux_matrix);
0089             end
0090         end
0091     end
0092     [FEM_MAT{icmpt}.K,volumes] = stiffness_matrixP1_3D(elements',coordinates',DIFF_cmpts(icmpt));
0093     FEM_MAT{icmpt}.M = mass_matrixP1_3D(elements',volumes);
0094     FEM_MAT{icmpt}.A = mass_matrixP1_3D(elements',volumes,GX');
0095     
0096     % calculate the IC
0097     IC_Pts{icmpt} = IC_cmpts(icmpt)*ones(size(mymesh.Pts_cmpt_reorder{icmpt},2),1);
0098 end
0099 
0100 if (DO_COUPLING == yes)
0101     tic
0102     [FEMcouple_MAT,FEMcouple_ind0,FEMcouple_indf] ...
0103         = generate_FEM_coupling(FEM_MAT,mymesh.Ncmpt,mymesh.Nboundary,...
0104         mymesh.Pts_cmpt_reorder,mymesh.Ele_cmpt_reorder,mymesh.Pts_ind,mymesh.Pts_boundary_reorder,mymesh.Fac_boundary_reorder);
0105     %stop
0106     toc
0107     % IC for coupled case
0108     IC_couple = zeros(size(FEMcouple_MAT.M,1),1);
0109     for icmpt = 1:mymesh.Ncmpt
0110         IC_couple(FEMcouple_ind0(icmpt):FEMcouple_indf(icmpt),1) = IC_Pts{icmpt};
0111     end
0112 else
0113     FEMcouple_MAT = [];
0114     FEMcouple_ind0 = [];
0115     FEMcouple_indf = [];    
0116 end
0117 
0118 %% solve ODE
0119 for iexperi = 1:nexperi
0120     SDELTA = sdeltavec(iexperi);
0121     BDELTA = bdeltavec(iexperi);
0122     TE = SDELTA+BDELTA;
0123     SEQ = seqvec(iexperi);% for choosing case PGSE, OGSEcos or OGSEsin
0124     omega = 2*pi*npervec(iexperi)/SDELTA;
0125     OGSEPER = 1./omega*2*pi;%% set up number for OGSE
0126     
0127     disp(['    Experiment ',num2str(iexperi)]);
0128     %disp([SDELTA,BDELTA,npervec(iexperi)]);
0129     
0130     TLIST = [0,TE];
0131     for ib = 1:nb
0132         b_start_time = clock;
0133         % global variable setting QVAL for ODE time stepping
0134         QVAL = qvalues(iexperi,ib);               
0135         disp(['      qvalue ',num2str(QVAL,'%.1e')]);        
0136         difftime(iexperi) = seqdifftime;
0137         
0138         %% Solving for case of coupling between compartments.
0139         if (DO_COUPLING == yes)
0140             FEM_M = FEMcouple_MAT.M;
0141             FEM_K = FEMcouple_MAT.K;
0142             FEM_A = FEMcouple_MAT.A; %*seqprofile(state.time)*QVAL;
0143             FEM_Q = FEMcouple_MAT.Q;
0144             FEM_G = sparse(zeros(size(FEM_M,1),1));
0145             
0146             options = odeset('Mass',FEM_M,'AbsTol',ODEsolve_atol,'RelTol',ODEsolve_rtol,'Vectorized','on','Stats','off',...
0147                 'Jacobian',@odejac_bt_includeb);            
0148             disp('***Coupled: start ode solve ode23t'); tic
0149             sol = ode23t(@odefun_bt_includeb,TLIST,IC_couple,options);
0150             disp('***Coupled: end ode solve ode23t'); toc
0151             for icmpt = 1:mymesh.Ncmpt
0152                 YOUT{iexperi}{ib}{icmpt} = sol.y(FEMcouple_ind0(icmpt):FEMcouple_indf(icmpt),:);
0153                 TOUT{iexperi}{ib}{icmpt} = sol.x;
0154                 MT{iexperi}{ib}{icmpt} = sum(FEM_MAT{icmpt}.M*YOUT{iexperi}{ib}{icmpt},1);
0155             end
0156         else
0157             %% Solving for case of no coupling between compartments.
0158             for icmpt = 1:mymesh.Ncmpt
0159                 FEM_M = FEM_MAT{icmpt}.M;
0160                 FEM_K = FEM_MAT{icmpt}.K;
0161                 FEM_A = FEM_MAT{icmpt}.A;
0162                 FEM_Q = FEM_MAT{icmpt}.Q;
0163                 FEM_G = sparse(zeros(size(FEM_M,1),1));
0164                 
0165                 options = odeset('Mass',FEM_M,'AbsTol',ODEsolve_atol,'RelTol',ODEsolve_rtol,'Vectorized','on','Stats','off',...
0166                     'Jacobian',@odejac_bt_includeb);
0167                 %disp('***Uncoupled: start ode solver ode23t'); tic
0168                 ICC = IC_Pts{icmpt};
0169                 if (max(abs(ICC))<=1e-16)
0170                     sol.y = zeros(size(ICC,1),size(TLIST,2));
0171                     sol.x = TLIST;
0172                 else
0173                     sol = ode23t(@odefun_bt_includeb,TLIST,ICC,options);
0174                 end
0175                 %disp('***Uncoupled: end ode solver ode23t'); toc
0176                 YOUT{iexperi}{ib}{icmpt} = sol.y;
0177                 TOUT{iexperi}{ib}{icmpt} = sol.x;
0178                 MT{iexperi}{ib}{icmpt} = sum(FEM_MAT{icmpt}.M*YOUT{iexperi}{ib}{icmpt},1);
0179             end            
0180         end
0181         elapsed_time(ib, iexperi)=etime(clock, b_start_time);
0182     end
0183 end
0184 
0185 for iexperi = 1:nexperi
0186     bvec = bvalues(iexperi,:);  
0187     nb = length(bvec);
0188     for ib = 1:nb
0189         for icmpt = 1:mymesh.Ncmpt
0190             MF_cmpts(icmpt,iexperi,ib) = MT{iexperi}{ib}{icmpt}(end);
0191             M0(icmpt,iexperi,ib) = MT{iexperi}{ib}{icmpt}(1);
0192         end
0193         MF_allcmpts(iexperi,ib) = 0;
0194         for icmpt = 1:mymesh.Ncmpt
0195             MF_allcmpts(iexperi,ib) = MF_allcmpts(iexperi,ib) + MF_cmpts(icmpt,iexperi,ib);
0196         end
0197         M0_allcmpts(iexperi,ib) = 0;
0198         for icmpt = 1:mymesh.Ncmpt
0199             M0_allcmpts(iexperi,ib) = M0_allcmpts(iexperi,ib) + M0(icmpt,iexperi,ib);
0200         end
0201     end   
0202     ib0 = find(abs(bvec)<=1e-16);
0203     ibn0 = find(abs(bvec)>1e-16);
0204     if (length(ib0) >= 1)
0205         for icmpt = 1:mymesh.Ncmpt
0206             S0(icmpt,iexperi) = MF_cmpts(icmpt,iexperi,ib0(1));
0207         end
0208         S0_allcmpts(iexperi) = MF_allcmpts(iexperi,ib0(1));
0209     else
0210         S0(1:mymesh.Ncmpt,iexperi) = nan;
0211         S0_allcmpts(iexperi) = nan;
0212     end
0213 end
BTPDE

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
