Comments and print statements updated in terzaghi1D

Author: JananiPSrinivasan

examples/matlab/terzaghi1D.m

@@ -12,14 +12,18 @@
 m = 50;            % Number of cells
 
 % Spatial discretization
-a = 0; 
-b = l;
-dx = (b - a)/m;
-xgrid = [a a+dx/2 : dx : b-dx/2 b];  % Staggered grid
+a = 0;                  % Left boundary of the domain [m]
+b = l;                  % Right boundary of the domain [m]
+dx = (b - a)/m;         % Cell width (uniform grid spacing) [m]
+xgrid = [a a+dx/2 : dx : b-dx/2 b];  % Staggered grid with ghost nodes at boundaries
+% xgrid has m+2 points: 
+% - ghost cell at a (Dirichlet BC)
+% - m internal nodes
+% - ghost cell at b (Neumann BC)
 
 % Times to evaluate (in hours)
-times_hr = [1, 10, 40, 70];
-times_sec = times_hr * 3600;
+times_hr = [1, 10, 40, 70];         % Time snapshots in hours for comparison
+times_sec = times_hr * 3600;        % Convert time points to seconds for simulation
 
 % Fluid properties
 K = 1e-12;        % Permeability [m^2]
@@ -32,41 +36,78 @@
 G = grad(k, m, dx);
 L(1,:) = 0; L(end,:) = 0;
 
-p_numerical = zeros(length(xgrid), length(times_sec));
+p_numerical = zeros(length(xgrid), length(times_sec)); % Pressure field [Pa]
+q_numerical = zeros(m+1, length(times_sec));  % Darcy flux [m/s] (size = m+1)
+
+% Loop over each specified final time
 for i = 1:length(times_sec)
     t_final = times_sec(i);
     dt = 0.03;
     nsteps = round(t_final / dt);
 
-    % Uniform Initial Condition
+    % Uniform Initial Condition : p(x,0) = P0
     p = P0 * ones(size(xgrid))';
 
+    % Time integration using Forward Euler
     for step = 1:nsteps
         p = p + dt * cf * (L * p);
-        p(1) = 0;             % Dirichlet BC
-        p(end) = p(end-1);    % Neumann BC
+        p(1) = 0;             % Dirichlet BC at x = 0
+        p(end) = p(end-1);    % Neumann BC at x = L
     end
 
     p_numerical(:,i) = p;
+
+    % Compute Darcy flux (numerical)
+    dpdx = G * p;
+    q = - (K / mu) * (dpdx - rho * g);
+    q_numerical(:,i) = q(1:m+1);  
+
+    % Print numerical results
+    fprintf('\nNumerical results at t = %.2f hr:\n', times_hr(i));
+    fprintf('|     x (m)     | Numerical p [Pa] | Darcy Flux [m/s] |\n');
+    fprintf('|---------------|------------------|------------------|\n');
+    for j = 1:m+1
+        fprintf('| %13.6f | %16.6e | %16.6e |\n', xgrid(j), p_numerical(j,i), q_numerical(j,i));
+    end
 end
 
 %% Analytical Solution
-N_max = 100;
-p_analytical = zeros(length(xgrid), length(times_sec));
+N_max = 100; % Number of Fourier series terms
+p_analytical = zeros(length(xgrid), length(times_sec)); % Pressure field [Pa]
+q_analytical = zeros(m+1, length(times_sec));  % Darcy flux [m/s] (size = m+1)
+
+% Loop over all time snapshots
 for i = 1:length(times_sec)
     t = times_sec(i);
     p = zeros(size(xgrid));
+
+     % Compute analytical solution using truncated Fourier sine series
     for N = 0:N_max
-        n = 2*N + 1;
+        n = 2*N + 1; % Odd terms only (satisfies boundary conditions)
         term = (1/n) * sin(n*pi*xgrid/(2*l)) .* exp(-(n^2)*(pi^2)*cf*t/(4*l^2));
         p = p + term;
     end
-    p_analytical(:,i) = (4/pi) * P0 * p;
+    p = (4/pi) * P0 * p;
+    p_analytical(:,i) = p;
+
+    % Compute Darcy flux (analytical)
+    dpdx = gradient(p, dx);  % basic gradient
+    q = - (K / mu) * (dpdx - rho * g);
+    q_analytical(:,i) = q(1:m+1);
+
+    % Print analytical results
+    fprintf('\nAnalytical results at t = %.2f hr:\n', times_hr(i));
+    fprintf('|     x (m)     | Analytical p [Pa] | Darcy Flux [m/s] |\n');
+    fprintf('|---------------|-------------------|------------------|\n');
+    for j = 1:m+1
+        fprintf('| %13.6f | %18.6e | %16.6e |\n', xgrid(j), p_analytical(j,i), q_analytical(j,i));
+    end
 end
 
 %% Combined Plot
 figure('Name','Terzaghi one-dimensional consolidation');
 set(gcf,'Color','white');
+
 % Top subplot: MOLE numerical
 subplot(2,1,1);
 hold on;
@@ -94,7 +135,15 @@
 
 %% Print relative L2 errors
 fprintf('\nRelative L2 Errors (Numerical vs Analytical):\n');
+fprintf('|     Time [hr]    |   Pressure Error   |   Darcy Flux Error |\n');
+fprintf('|------------------|--------------------|---------------------|\n');
 for i = 1:length(times_hr)
-    rel_err = norm(p_numerical(:,i) - p_analytical(:,i)) / norm(p_analytical(:,i));
-    fprintf('t = %.2f hr: %.6e\n', times_hr(i), rel_err);
+    % L2 error for pressure
+    rel_err_p = norm(p_numerical(:,i) - p_analytical(:,i)) / norm(p_analytical(:,i));
+    
+    % L2 error for flux
+    rel_err_q = norm(q_numerical(:,i) - q_analytical(:,i)) / norm(q_analytical(:,i));
+    
+    % Print in table format
+    fprintf('| %16.2f | %18.6e | %19.6e |\n', times_hr(i), rel_err_p, rel_err_q);
 end