updated terzaghi1D

Author: JananiPSrinivasan

examples/matlab/compare_terzaghi_solutions.m

@@ -1,96 +1,100 @@
-% Numerical solution of Terzaghi's 1D consolidation equation using MOLE
+% Terzaghi 1D: using MOLE operators
 clc;
 close all;
 
-addpath('../../src/matlab');  % Update this path to your MOLE operators
+addpath('../../src/matlab');  % MOLE operator path
 
-% Parameters
+%% Parameters
 P0 = 10e6;         % Face load in Pascals
-cf = 1e-4;         % Consolidation (diffusion) constant
-L = 25;            % Domain length in meters
-k = 2;             % Order of accuracy
+cf = 1e-4;         % Diffusion constant
+l = 25;            % Domain length in meters
+k = 2;             % MOLE operator order
 m = 50;            % Number of cells
-a = 0; b = L;
+
+% Spatial discretization
+a = 0; 
+b = l;
 dx = (b - a)/m;
 xgrid = [a a+dx/2 : dx : b-dx/2 b];  % Staggered grid
 
 % Times to evaluate (in hours)
 times_hr = [1, 10, 40, 70];
 times_sec = times_hr * 3600;
 
-% Fluid and medium properties
-K = 1e-12;       % Permeability [m^2]
-mu = 1e-3;       % Dynamic viscosity [Pa·s]
-rho = 1000;      % Fluid density [kg/m^3]
-g = 9.81;        % Gravity [m/s^2]
-
-% MOLE operators
-L = lap(k, m, dx);  % Laplacian
-G = grad(k, m, dx);   % Gradient
+% Fluid properties
+K = 1e-12;        % Permeability [m^2]
+mu = 1e-3;        % Dynamic viscosity [Pa·s]
+rho = 1000;       % Fluid density [kg/m^3]
+g = 9.81;         % Gravity [m/s^2]
 
-% Enforce zero rows for boundary (to avoid Laplacian affecting boundary values)
-L(1,:) = 0;
-L(end,:) = 0;
+%% Numerical (MOLE) Solution
+L = lap(k, m, dx);
+G = grad(k, m, dx);
+L(1,:) = 0; L(end,:) = 0;
 
-% Store results
-p_all = zeros(length(xgrid), length(times_sec));
-q_all = zeros(length(xgrid), length(times_sec));
-
-% Plot pressure
-figure(1); hold on;
+p_numerical = zeros(length(xgrid), length(times_sec));
 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 = P0 * ones(size(xgrid))';
 
-    % Time stepping loop
     for step = 1:nsteps
         p = p + dt * cf * (L * p);
-        % Boundary Conditions
-        p(1) = 0;              % Dirichlet at x = 0
-        p(end) = p(end-1);     % Neumann at x = L (zero gradient)
+        p(1) = 0;             % Dirichlet BC
+        p(end) = p(end-1);    % Neumann BC
     end
-    
-    % Store final pressure
-    p_all(:,i) = p;
 
-    % Compute gradient and Darcy flux
-    dpdx = G * p;
-    q = - (K / mu) * (dpdx - rho * g);
-    q_all(1:m+1,i) = q(1:m+1);  % Exclude ghost point
+    p_numerical(:,i) = p;
+end
 
-    % Plot pressure profile
-    semilogy(xgrid, p/1e6, '-o', 'DisplayName', ['t = ' num2str(times_hr(i)) ' hr']);
+%% Analytical Solution
+N_max = 100;
+p_analytical = zeros(length(xgrid), length(times_sec));
+for i = 1:length(times_sec)
+    t = times_sec(i);
+    p = zeros(size(xgrid));
+    for N = 0:N_max
+        n = 2*N + 1;
+        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;
 end
-title('Pressure p(x,t) at various times');
-xlabel('x (m)');
+
+%% Combined Plot
+figure('Name','Terzaghi one-dimensional consolidation');
+set(gcf,'Color','white');
+% Top subplot: MOLE numerical
+subplot(2,1,1);
+hold on;
+for i = 1:length(times_sec)
+    semilogy(xgrid, p_numerical(:,i)/1e6, '-o', 'DisplayName', [num2str(times_hr(i)) ' hr']);
+end
+title('MOLE Numerical Solution');
 ylabel('Excess Pore Pressure p(x,t) [MPa]');
+axis([0 l 1e-3 10]);
 legend('Location', 'southeast');
 grid on;
 
-% Plot Darcy flux
-figure(2); hold on;
-x_flux = linspace(a + dx/2, b - dx/2, m+1);
+% Bottom subplot: Analytical
+subplot(2,1,2);
+hold on;
 for i = 1:length(times_sec)
-    q = q_all(1:m+1,i);
-    semilogy(x_flux, q, '-s', 'DisplayName', ['t = ' num2str(times_hr(i)) ' hr']);
+    semilogy(xgrid, p_analytical(:,i)/1e6, '--s', 'DisplayName', [num2str(times_hr(i)) ' hr']);
 end
-title('Darcy Flux q(x,t) at various times');
+title('Analytical Benchmark Solution');
 xlabel('x (m)');
-ylabel('q(x,t) [m/s]');
+ylabel('Excess Pore Pressure p(x,t) [MPa]');
+axis([0 l 1e-3 10]);
 legend('Location', 'southeast');
 grid on;
 
-% Print final values
+%% Print relative L2 errors
+fprintf('\nRelative L2 Errors (Numerical vs Analytical):\n');
 for i = 1:length(times_hr)
-    fprintf('\nFinal pressure p(x,t = %.2f hr):\n', times_hr(i));
-    disp(p_all(:,i)');
-    fprintf('Final Darcy flux q(x,t = %.2f hr):\n', times_hr(i));
-    disp(q_all(1:m+1,i)');
+    rel_err = norm(p_numerical(:,i) - p_analytical(:,i)) / norm(p_analytical(:,i));
+    fprintf('t = %.2f hr: %.6e\n', times_hr(i), rel_err);
 end
-
-% Save pressure for comparison
-p_numerical = p_all;