lots of bugs fixed for the 1d spherical

simulkade · Mar 2, 2018 · 7681b2b · 7681b2b
1 parent b08d3d0
commit 7681b2b
Show file tree

Hide file tree

Showing 26 changed files with 41 additions and 41 deletions.
diff --git a/Boundary/boundaryCondition.m b/Boundary/boundaryCondition.m
@@ -30,7 +30,7 @@
 % See the license file
 
 d = BC.domain.dimension;
-if (d ==1) || (d==1.5)
+if (d ==1) || (d==1.5) || (d==1.8)
 	[BCMatrix, BCRHS] = boundaryCondition1D(BC);
 elseif (d == 2) || (d == 2.5)
 	[BCMatrix, BCRHS] = boundaryCondition2D(BC);

diff --git a/Boundary/cellBoundary.m b/Boundary/cellBoundary.m
@@ -22,7 +22,7 @@
 % extract data from the mesh structure
 d = BC.domain.dimension;
 
-if (d ==1) || (d==1.5)
+if (d ==1) || (d==1.5) || (d==1.8)
 	phiBC = cellBoundary1D(phi, BC);
 elseif (d == 2) || (d == 2.5)
 	phiBC = cellBoundary2D(phi, BC);

diff --git a/Boundary/createBC.m b/Boundary/createBC.m
@@ -25,7 +25,7 @@
 % See the license file
 
 d = meshvar.dimension;
-if (d == 1) || (d==1.5)
+if (d == 1) || (d==1.5) || (d==1.8)
 	BC = createBC1D(meshvar);
 elseif (d == 2) || (d == 2.5) || (d == 2.8)
 	BC = createBC2D(meshvar);

diff --git a/Calculus/ddtTerm.m b/Calculus/ddtTerm.m
@@ -23,7 +23,7 @@
 % See the license file
 
 d = MeshStructure.dimension;
-if (d ==1) || (d==1.5)
+if (d ==1) || (d==1.5) || (d==1.8)
 	ddt = ddtTerm1D(MeshStructure, dt, phi);
 elseif (d == 2) || (d == 2.5) || (d==2.8)
 	ddt = ddtTerm2D(MeshStructure, dt, phi);

diff --git a/Calculus/gradientCellTerm.m b/Calculus/gradientCellTerm.m
@@ -21,7 +21,7 @@
 
 phi_face=linearMean(phi);
 d = phi.domain.dimension;
-if (d ==1) || (d==1.5)
+if (d ==1) || (d==1.5) || (d==1.8)
 	Nx = phi.domain.dims(1);
 	DX = phi.domain.cellsize.x(2:end-1);
 	xvalue = (phi_face.xvalue(2:Nx+1)-phi_face.xvalue(1:Nx))./DX;

diff --git a/Calculus/gradientTerm.m b/Calculus/gradientTerm.m
@@ -21,7 +21,7 @@
 % See the license file
 
 d = phi.domain.dimension;
-if (d ==1) || (d==1.5)
+if (d ==1) || (d==1.5) || (d==1.8)
 	faceGrad = gradientTerm1D(phi);
 elseif (d == 2) || (d == 2.5)
 	faceGrad = gradientTerm2D(phi);

diff --git a/Discretization/constantSourceTerm.m b/Discretization/constantSourceTerm.m
@@ -20,7 +20,7 @@
 % See the license file
 
 d = phi.domain.dimension;
-if (d ==1) || (d==1.5)
+if (d ==1) || (d==1.5) || (d==1.8)
 	RHS = constantSourceTerm1D(phi);
 elseif (d == 2) || (d == 2.5) || (d==2.8)
 	RHS = constantSourceTerm2D(phi);

diff --git a/Discretization/diffusionTermSpherical1D.m b/Discretization/diffusionTermSpherical1D.m
@@ -21,38 +21,38 @@
 % See the license file
 
 % extract data from the mesh structure
-Nx = D.domain.dims(1);
-G = 1:Nx+2;
+Nr = D.domain.dims(1);
+G = 1:Nr+2;
 DX = D.domain.cellsize.x;
 dx = 0.5*(DX(1:end-1)+DX(2:end));
 % rp = D.domain.cellcenters.x;
 rf = D.domain.facecenters.x;
 
 % define the vectors to store the sparse matrix data
-iix = zeros(3*(Nx+2),1);
-jjx = zeros(3*(Nx+2),1);
-sx = zeros(3*(Nx+2),1);
+iix = zeros(3*(Nr+2),1);
+jjx = zeros(3*(Nr+2),1);
+sx = zeros(3*(Nr+2),1);
 
 % extract the velocity data
 % note: size(Dx) = [1:m+1, 1:n] and size(Dy) = [1:m, 1:n+1]
 Dx = D.xvalue;
 
 % reassign the east, west, north, and south velocity vectors for the
 % code readability
-De = Dx(2:Nx+1).*rf(2:Nx+1).^2./(1/3*(rf(2:Nr+1).^3-rf(1:Nr).^3).*dx(2:Nx+1));
-Dw = Dx(1:Nx).*rf(1:Nx).^2./(1/3*(rf(2:Nr+1).^3-rf(1:Nr).^3).*dx(1:Nx));
+De = Dx(2:Nr+1).*rf(2:Nr+1).^2./(1/3*(rf(2:Nr+1).^3-rf(1:Nr).^3).*dx(2:Nr+1));
+Dw = Dx(1:Nr).*rf(1:Nr).^2./(1/3*(rf(2:Nr+1).^3-rf(1:Nr).^3).*dx(1:Nr));
 
 % calculate the coefficients for the internal cells
-AE = reshape(De,Nx,1);
-AW = reshape(Dw,Nx,1);
-APx = reshape(-(De+Dw),Nx,1);
+AE = reshape(De,Nr,1);
+AW = reshape(Dw,Nr,1);
+APx = reshape(-(De+Dw),Nr,1);
 
 % build the sparse matrix based on the numbering system
-rowx_index = reshape(G(2:Nx+1),Nx,1); % main diagonal x
-iix(1:3*Nx) = repmat(rowx_index,3,1);
-jjx(1:3*Nx) = [reshape(G(1:Nx),Nx,1); reshape(G(2:Nx+1),Nx,1); reshape(G(3:Nx+2),Nx,1)];
-sx(1:3*Nx) = [AW; APx; AE];
+rowx_index = reshape(G(2:Nr+1),Nr,1); % main diagonal x
+iix(1:3*Nr) = repmat(rowx_index,3,1);
+jjx(1:3*Nr) = [reshape(G(1:Nr),Nr,1); reshape(G(2:Nr+1),Nr,1); reshape(G(3:Nr+2),Nr,1)];
+sx(1:3*Nr) = [AW; APx; AE];
 
 % build the sparse matrix
-kx = 3*Nx;
-M = sparse(iix(1:kx), jjx(1:kx), sx(1:kx), Nx+2, Nx+2);
+kx = 3*Nr;
+M = sparse(iix(1:kx), jjx(1:kx), sx(1:kx), Nr+2, Nr+2);
diff --git a/Discretization/linearSourceTerm.m b/Discretization/linearSourceTerm.m
@@ -21,7 +21,7 @@
 % See the license file
 
 d = k.domain.dimension;
-if (d ==1) || (d==1.5)
+if (d ==1) || (d==1.5) || (d==1.5)
 	M = linearSourceTerm1D(k);
 elseif (d == 2) || (d == 2.5) || (d==2.8)
 	M = linearSourceTerm2D(k);

diff --git a/Discretization/transientTerm.m b/Discretization/transientTerm.m
@@ -24,7 +24,7 @@
 % See the license file
 
 d = phi_old.domain.dimension;
-if (d ==1) || (d==1.5)
+if (d ==1) || (d==1.5) || (d==1.8)
 	% extract data from the mesh structure
     Nx = phi_old.domain.dims(1);
     G = 1:Nx+2;

diff --git a/MeshGeneration/createFaceVariable.m b/MeshGeneration/createFaceVariable.m
@@ -23,7 +23,7 @@
 d = meshvar.dimension;
 mn = meshvar.dims;
 
-if (d ==1) || (d==1.5)
+if (d ==1) || (d==1.5) || (d==1.8)
 	xvalue = faceval(1).*ones(mn(1)+1,1);
     yvalue=[];
     zvalue=[];

diff --git a/MeshGeneration/createMeshSpherical1D.m b/MeshGeneration/createMeshSpherical1D.m
@@ -74,7 +74,7 @@
   FaceLocation.z= [0.0];
 end
 
-MS=MeshStructure(1.5, ...
+MS=MeshStructure(1.8, ...
   [Nr,1], ...
   CellSize, ...
   CellLocation, ...

diff --git a/Utilities/arithmeticMean.m b/Utilities/arithmeticMean.m
@@ -23,7 +23,7 @@
 % extract data from the mesh structure
 
 d = phi.domain.dimension;
-if (d ==1) || (d==1.5)
+if (d ==1) || (d==1.5) || (d==1.8)
     dx = phi.domain.cellsize.x;
     xvalue=(dx(1:end-1).*phi.value(1:end-1)+dx(2:end).*phi.value(2:end))./(dx(2:end)+dx(1:end-1));
     yvalue=[];

diff --git a/Utilities/cellLocations.m b/Utilities/cellLocations.m
@@ -27,7 +27,7 @@
 Z= createCellVariable(m, 0);
 d = m.dimension;
 switch d
-        case {1, 1.5}
+case {1, 1.5, 1.8}
     	X=createCellVariable(m, m.cellcenters.x);
     case {2, 2.5, 2.8}
         X= createCellVariable(m, repmat(m.cellcenters.x, 1, N(2)));

diff --git a/Utilities/excludeGhostRHS.m b/Utilities/excludeGhostRHS.m
@@ -22,7 +22,7 @@
 % check the size of the variable and the mesh dimension
 d = MS.dimension;
 
-if (d ==1) || (d==1.5)
+if (d ==1) || (d==1.5) || (d==1.8)
 	Nx=MS.dims(1);
 	G = 1:Nx+2;
 	RHSout = RHS(reshape(G(2:end-1),Nx,1));

diff --git a/Utilities/faceLocations.m b/Utilities/faceLocations.m
@@ -26,7 +26,7 @@
 N=m.dims;
 d = m.dimension;
 switch d
-    case {1, 1.5}
+case {1, 1.5, 1.8}
     	X.xvalue= m.facecenters.x;
     case {2, 2.5, 2.8}
         X.xvalue= repmat(m.facecenters.x, 1, N(2));

diff --git a/Utilities/geometricMean.m b/Utilities/geometricMean.m
@@ -23,7 +23,7 @@
 % extract data from the mesh structure
 
 d = phi.domain.dimension;
-if (d ==1) || (d==1.5)
+if (d ==1) || (d==1.5) || (d==1.8)
     dx = phi.domain.cellsize.x;
     xvalue=exp((dx(1:end-1).*log(phi.value(1:end-1))+dx(2:end).*log(phi.value(2:end)))./(dx(2:end)+dx(1:end-1)));
     yvalue=[];

diff --git a/Utilities/harmonicMean.m b/Utilities/harmonicMean.m
@@ -23,7 +23,7 @@
 % extract data from the mesh structure
 
 d = phi.domain.dimension;
-if (d ==1) || (d==1.5)
+if (d ==1) || (d==1.5) || (d==1.8)
     dx = phi.domain.cellsize.x;
     xvalue=phi.value(1:end-1).*phi.value(2:end).*(dx(2:end)+dx(1:end-1))...
         ./(dx(2:end).*phi.value(1:end-1)+dx(1:end-1).*phi.value(2:end));

diff --git a/Utilities/internalCells.m b/Utilities/internalCells.m
@@ -23,7 +23,7 @@
 d = phi.domain.dimension;
 N = phi.domain.dims;
 
-if (d ==1) || (d==1.5)
+if (d ==1) || (d==1.5) || (d==1.8)
 	cellvar = phi.value(2:N(1)+1);
 elseif (d == 2) || (d == 2.5) || (d==2.8)
 	cellvar = phi.value(2:N(1)+1, 2:N(2)+1);

diff --git a/Utilities/linearMean.m b/Utilities/linearMean.m
@@ -23,7 +23,7 @@
 % extract data from the mesh structure
 
 d = phi.domain.dimension;
-if (d ==1) || (d==1.5)
+if (d ==1) || (d==1.5) || (d==1.8)
     dx = phi.domain.cellsize.x;
     xvalue=(dx(2:end).*phi.value(1:end-1)+dx(1:end-1).*phi.value(2:end))./(dx(2:end)+dx(1:end-1));
     yvalue=[];

diff --git a/Utilities/maskCells.m b/Utilities/maskCells.m
@@ -23,7 +23,7 @@
 d = meshstruct.dimension;
 domain_size = meshstruct.numberofcells+2; % 2 is added for the ghost cells
 M_size = size(M);
-if d ==1 || d==1.5
+if d ==1 || d==1.5 || (d==1.8)
     i = sub2ind([domain_size 1], cellIndex(:,1));
 elseif d==2 || d==2.5 || d==2.8
     i = sub2ind(domain_size, cellIndex(:,1), cellIndex(:,2));

diff --git a/Utilities/normCellVector.m b/Utilities/normCellVector.m
@@ -21,7 +21,7 @@
 % check the size of the variable and the mesh dimension
 d = cellvec.domain.dimension;
 
-if (d ==1) || (d==1.5)
+if (d ==1) || (d==1.5) || (d==1.8)
 	cellnormval = abs(cellvec.xvalue);
 elseif (d == 2) || (d == 2.5) || (d==2.8)
 	cellnormval = realsqrt(cellvec.xvalue.*cellvec.xvalue+ ...

diff --git a/Utilities/reshapeCell.m b/Utilities/reshapeCell.m
@@ -24,7 +24,7 @@
 d = MS.dimension;
 N = MS.dims;
 
-if (d ==1) || (d==1.5)
+if (d ==1) || (d==1.5) || (d==1.8)
 	cellvar = reshape(phi, N(1)+2, 1);
 elseif (d == 2) || (d == 2.5) || (d==2.8)
 	cellvar = reshape(phi, N(1)+2, N(2)+2);

diff --git a/Utilities/reshapeInternalCell.m b/Utilities/reshapeInternalCell.m
@@ -25,7 +25,7 @@
 d = MS.dimension;
 N = MS.dims;
 
-if (d ==1) || (d==1.5)
+if (d ==1) || (d==1.5) || (d==1.8)
 	cellvar = reshape(phi, N(1), 1);
 elseif (d == 2) || (d == 2.5) || (d==2.8)
 	cellvar = reshape(phi, N(1), N(2));

diff --git a/Utilities/tvdMean.m b/Utilities/tvdMean.m
@@ -24,7 +24,7 @@
 % extract data from the mesh structure
 
 d = phi.domain.dimension;
-if (d ==1) || (d==1.5)
+if (d ==1) || (d==1.5) || (d==1.8)
 	phiFaceAverage = tvdMean1D(phi, u, FL);
 elseif (d == 2) || (d == 2.5) || (d==2.8)
 	phiFaceAverage = tvdMean2D(phi, u, FL);

diff --git a/Utilities/upwindMean.m b/Utilities/upwindMean.m
@@ -24,7 +24,7 @@
 % extract data from the mesh structure
 
 d = phi.domain.dimension;
-if (d ==1) || (d==1.5)
+if (d ==1) || (d==1.5) || (d==1.8)
 	phiFaceAverage = upwindMean1D(phi, u);
 elseif (d == 2) || (d == 2.5) || (d==2.8)
 	phiFaceAverage = upwindMean2D(phi, u);