Uploaded R03

New kernel ('cell-centered'); fixed memory leak; added CG restart
lanl · Nov 5, 2024 · 50b6788 · 50b6788
1 parent 0f9024c
commit 50b6788
Show file tree

Hide file tree

Showing 3 changed files with 35 additions and 103 deletions.
diff --git a/OmniMatlab/example/Test_OSMODI.m b/OmniMatlab/example/Test_OSMODI.m
@@ -1,31 +1,3 @@
-% © 2024. Triad National Security, LLC. All rights reserved.
-% This program was produced under U.S. Government contract
-% 89233218CNA000001 for Los Alamos National Laboratory (LANL), which is
-% operated by Triad National Security, LLC for the U.S. Department of
-% Energy/National Nuclear Security Administration. All rights in the
-% program are reserved by Triad National Security, LLC, and the U.S.
-% Department of Energy/National Nuclear Security Administration. The
-% Government is granted for itself and others acting on its behalf a
-% nonexclusive, paid-up, irrevocable worldwide license in this material
-% to reproduce, prepare. derivative works, distribute copies to the
-% public, perform publicly and display publicly, and to permit
-% others to do so.
-%
-% This program is free software: you can redistribute it and/or modify it
-% under the terms of the GNU General Public License as published by the
-% Free Software Foundation, either version 3 of the License, or (at your
-% option) any later version.
-%
-% This program is distributed in the hope that it will be useful, but
-% WITHOUT ANY WARRANTY; without even the implied warranty of
-% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-% General Public License for more details.
-%
-% You should have received a copy of the GNU General Public License along
-% with this program. If not, see <https://www.gnu.org/licenses/>.
-
-%Minimal Working Example code for usage example
-
 clear; clc; close all;
 
 Nx=10; Ny=20; Nz=30;
@@ -38,7 +10,11 @@
 opts.SolverToleranceAbs=1e-4;
 opts.SolverDevice='GPU';
 opts.Verbose=1;
+opts.Kernel='cell-centered';
+%opts.Kernel='face-crossing';
 
 [P, CGS]=OSMODI(single(Sx),single(Sy),single(Sz),single(delta),opts);
 
+sum(P(:)-single(Sz(:)))
+
 
diff --git a/OmniMatlab/example/Test_Taylor2D.m b/OmniMatlab/example/Test_Taylor2D.m
@@ -1,35 +1,8 @@
-% © 2024. Triad National Security, LLC. All rights reserved.
-% This program was produced under U.S. Government contract
-% 89233218CNA000001 for Los Alamos National Laboratory (LANL), which is
-% operated by Triad National Security, LLC for the U.S. Department of
-% Energy/National Nuclear Security Administration. All rights in the
-% program are reserved by Triad National Security, LLC, and the U.S.
-% Department of Energy/National Nuclear Security Administration. The
-% Government is granted for itself and others acting on its behalf a
-% nonexclusive, paid-up, irrevocable worldwide license in this material
-% to reproduce, prepare. derivative works, distribute copies to the
-% public, perform publicly and display publicly, and to permit
-% others to do so.
-%
-% This program is free software: you can redistribute it and/or modify it
-% under the terms of the GNU General Public License as published by the
-% Free Software Foundation, either version 3 of the License, or (at your
-% option) any later version.
-%
-% This program is distributed in the hope that it will be useful, but
-% WITHOUT ANY WARRANTY; without even the implied warranty of
-% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-% General Public License for more details.
-%
-% You should have received a copy of the GNU General Public License along
-% with this program. If not, see <https://www.gnu.org/licenses/>.
-
-
 clear; clc; close all;
 
 %Generates a Taylor vortex according to Charonko 2010 DOI 10.1088/0957-0233/21/10/105401
 
-t=.3; %Time, s
+t=.6; %Time, s
 dt=0.001; %dt for time derivative
 H=1e-6; %Vortex strength, m^2
 nu=1e-6; %Kinematic Viscosity, m^2/s
@@ -40,9 +13,9 @@
 T0=L0/U0; %Characteristic Time, s
 P0=rho*U0^2; %Characteristic Pressure, Pa
 
-NptsX=1001; NptsY=NptsX; r=NptsY/NptsX;
+NptsX=2101; NptsY=NptsX; r=NptsY/NptsX;
 
-DomainSize=3;
+DomainSize=2;
 x=linspace(-DomainSize*L0,DomainSize*L0,NptsX); y=linspace(-DomainSize*L0,DomainSize*L0,NptsY); dx=x(2)-x(1); dy=y(2)-y(1);
 [X,Y]=ndgrid(x,y);
 R=sqrt(X.^2+Y.^2);
@@ -67,34 +40,37 @@
 Sy=(-rho*(dvdt+Ux.*DVDX+Uy.*DVDY))/(P0/L0);
 Sz=zeros(size(Sx));
 
-%% ----THIS IS WHERE THE OSMODI CODE IS USED----
-% Make sure the OSMODI.mexw64 file is inside this folder or inside of a Matlab search folder
 opts.SolverToleranceRel=1e-4;
 opts.SolverToleranceAbs=1e-6;
 opts.SolverDevice='GPU';
-opts.Verbose=0;
+opts.Verbose=1;
+%opts.Kernel='face-crossing';
+opts.Kernel='cell-centered';
 
-[P_OSMODI, CGS]=OSMODI(single(Sx),single(Sy),single(Sz),single([dx dy]/L0),opts);
+%Sx(10:20,20:40)=nan;
 
-disp(['Time Taken for just GPU computation: ' num2str(CGS(end,3)-CGS(1,3)) 's'])
+%[P_OSMODI, CGS]=OSMODI_Kernels(single(Sx),single(Sy),single(Sz),single([dx dy]/L0),opts);
+%[P_OSMODI, CGS]=OSMODI(double(Sx),double(Sy),double(Sz),double([dx dy]/L0),opts);
+[P_OSMODI, CGS]=OSMODI(Sx,Sy,Sz,[dx dy]/L0,opts);
 
+disp(['Time Taken for just ' opts.SolverDevice ' computation: ' num2str(CGS(end,3)-CGS(1,3)) 's'])
 %%
 
 P_OSMODI=P_OSMODI-P_OSMODI(2,2); %Removes corner value 
 
-figure('color','w','position',[-1.2182e+03 421.8000 915.2000 361.6000]);
+figure('color','w');
 subplot(1,2,1)
 imagesc(P_OSMODI'); colorbar; title('Solver'); daspect([dy dx 1])
 subplot(1,2,2)
 imagesc(P'/(rho*U0^2)); colorbar; title('Truth'); daspect([dy dx 1])
 
 Ptruth=P/(rho*U0^2);
 Err=(P_OSMODI' - Ptruth')./NormPmax;
-figure('color','w'); imagesc(Err); colorbar; title('Error'); daspect([dy dx 1])
+figure('color','w'); imagesc(Err); colorbar; title(['Error (' opts.SolverDevice ')']); daspect([dy dx 1])
 
 ErrorPercent=100*sqrt(sum(Err(:).^2)./numel(Err));
 disp(['Error = ' num2str(ErrorPercent, '%0.3f') '%'])
 
-figure('color','w','Position', [1.3538e+03 530.6000 560 420]); semilogy(CGS(:,1), CGS(:,2))
-title('CG solver convergence')
+% figure('color','w'); semilogy(CGS(:,1), CGS(:,2))
+% title('CG solver convergence')
 
diff --git a/OmniMatlab/example/Test_Taylor3D.m b/OmniMatlab/example/Test_Taylor3D.m
@@ -1,30 +1,3 @@
-% © 2024. Triad National Security, LLC. All rights reserved.
-% This program was produced under U.S. Government contract
-% 89233218CNA000001 for Los Alamos National Laboratory (LANL), which is
-% operated by Triad National Security, LLC for the U.S. Department of
-% Energy/National Nuclear Security Administration. All rights in the
-% program are reserved by Triad National Security, LLC, and the U.S.
-% Department of Energy/National Nuclear Security Administration. The
-% Government is granted for itself and others acting on its behalf a
-% nonexclusive, paid-up, irrevocable worldwide license in this material
-% to reproduce, prepare. derivative works, distribute copies to the
-% public, perform publicly and display publicly, and to permit
-% others to do so.
-%
-% This program is free software: you can redistribute it and/or modify it
-% under the terms of the GNU General Public License as published by the
-% Free Software Foundation, either version 3 of the License, or (at your
-% option) any later version.
-%
-% This program is distributed in the hope that it will be useful, but
-% WITHOUT ANY WARRANTY; without even the implied warranty of
-% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-% General Public License for more details.
-%
-% You should have received a copy of the GNU General Public License along
-% with this program. If not, see <https://www.gnu.org/licenses/>.
-
-
 clear; clc; close all;
 
 %Generates a Taylor vortex according to Charonko 2010 DOI 10.1088/0957-0233/21/10/105401 3D version
@@ -39,7 +12,7 @@
 T0=L0/U0; %Characteristic Time, s
 P0=rho*U0^2; %Characteristic Pressure, Pa
 
-NptsX=101; NptsY=NptsX; NptsZ=NptsX;
+NptsX=201; NptsY=NptsX; NptsZ=NptsX;
 
 x=linspace(-2*L0,2*L0,NptsX); y=linspace(-2*L0,2*L0,NptsY); z=linspace(-2*L0,2*L0,NptsZ);
 dx=x(2)-x(1); dy=y(2)-y(1); dz=z(2)-z(1);
@@ -66,20 +39,26 @@
 Sy=(-rho*(dvdt+Ux.*DVDX+Uy.*DVDY))/(P0/L0);
 Sz=zeros(size(Sx));
 
-
-%% ----THIS IS WHERE THE OSMODI CODE IS USED----
-% Make sure the OSMODI.mexw64 file is inside this folder or inside of a Matlab search folder
-
 opts.SolverToleranceRel=1e-4;
 opts.SolverToleranceAbs=1e-6;
 opts.SolverDevice='GPU';
+%opts.Kernel='face-crossing';
+opts.Kernel='cell-centered';
 opts.Verbose=1;
 
+M=zeros(size(Sx));
+M(26:40,26:40,26:40)=1; %Makes sure the degeneracies are handled correctly
+M=M & rand(size(Sx))>0.5;
+Sx(M)=nan;
+
 disp('Starting OSMODI...')
 tic
-[P_OSMODI, CGS]=OSMODI(single(Sx),single(Sy),single(Sz),single([dx dy dz]/L0),opts);
+%[P_OSMODI, CGS]=OSMODI(single(Sx),single(Sy),single(Sz),single([dx dy dz]/L0),opts);
+[P_OSMODI, CGS]=OSMODI(Sx,Sy,Sz,[dx dy dz]/L0,opts);
 toc
 
+%figure; imagesc(squeeze(P_OSMODI(:,ceil(NptsX/2),:)));colorbar
+
 disp(['Time Taken for just GPU computation: ' num2str(CGS(end,3)-CGS(1,3)) 's'])
 
 
@@ -105,8 +84,9 @@
 ErrorPercent=100*sqrt(sum(Err(:).^2)./numel(Err));
 disp(['Error = ' num2str(ErrorPercent, '%0.3f') '%'])
 
-figure('color','w','Position', [1.3538e+03 530.6000 560 420]); semilogy(CGS(:,1), CGS(:,2))
-title('CG solver convergence')
+% figure('color','w','Position', [1.3538e+03 530.6000 560 420]); semilogy(CGS(:,1), CGS(:,2))
+% title('CG solver convergence')
 
-pause(2)
+% pause(2)
 
+clear;