quicwind/mat_flux.m

   1 function Mmat=mat_flux(X)
   2 % M=mat_flux(X)
   3 % Compute wind flux matrix for a hexa grid assuming all sides in the y
   4 % direction are straight vertical while horizontal sides may be slanted
   5 % in:
   6 %     X{l}(i,j,k) coordinate l of node (i,j,k), if l=1,2 does not depend on k
   7 % out:
   8 %     Mmat        matrix computing fluxes in normal directions on mesh
   9 %
  10
  11 % Mesh computations
  12 x = X{1}; y = X{2}; z = X{3};
  13 % Q: this assumes that we are working on a cube based on the X direction?
  14 [nx1, ny1, nz1] = size(x);
  15 nx = nx1 - 1; ny = ny1 - 1; nz = nz1 - 1;
  16 % TODO: error check to make sure we are only working in 3 dimensions?
  17 n = nx * ny * nz1 + nx * ny1 * nz + nx1 * ny * nz;
  18 s = cell_sizes(X);
  19 % continue one layer up
  20 s.dz_at_u(:,:,nz+1)=s.dz_at_u(:,:,nz);
  21 s.dz_at_v(:,:,nz+1)=s.dz_at_v(:,:,nz);
  22
  23
  24 % slope in x and y directions
  25 dzdx = zeros(nx,ny+1,nz+1);
  26 dzdy = zeros(nx+1,ny,nz+1);
  27 for k=1:nz+1
  28     for j=1:ny+1
  29         for i=1:nx
  30             dzdx(i,j,k) = ((z(i+1,j,k)-z(i,j,k))/(x(i+1,j,k)-x(i,j,k)));
  31         end
  32     end
  33     for j=1:ny
  34         for i=1:nx+1
  35             dzdy(i,j,k) = ((z(i,j+1,k)-z(i,j,k))/(y(i,j+1,k)-y(i,j,k)));
  36         end
  37     end
  38 end
  39
  40 nnz_est = n;
  41 I = zeros(nnz_est,1);
  42 J = zeros(nnz_est,1);
  43 M = zeros(nnz_est,1);
  44
  45 M_col = 1;
  46 M_row = 1;
  47 M_entry = 1;
  48 row_offset = nx1 * ny * nz + ny1 * nx * nz;
  49 for ux=1:numel(X)
  50
  51     % Size of the grid for the current direction
  52     grid_size = [nx,ny,nz];
  53     grid_size(ux) = grid_size(ux) + 1;
  54     X_len = grid_size(1); Y_len = grid_size(2); Z_len = grid_size(3);
  55
  56     for kx=1:Z_len
  57         for jx=1:Y_len
  58             for ix=1:X_len
  59                 if ux == 1
  60                     % U elements
  61                     M(M_entry) = s.area_u(ix,jx,kx);
  62                     I(M_entry) = M_row;
  63                     J(M_entry) = M_col;
  64                     M_entry = M_entry + 1;
  65
  66                     % Edge cases
  67                     is_top_level = (kx == Z_len);
  68                     is_low_level = (kx == 1);
  69                     is_first_row = (ix == 1);
  70                     is_last_row = (ix == X_len);
  71
  72                     if not(is_top_level)
  73                         if not(is_first_row)
  74                             M(M_entry) = -s.area_w(ix - 1, jx, kx + 1) * ...
  75                                 0.5 * (dzdx(ix - 1, jx, kx + 1) + dzdx(ix - 1, jx + 1, kx + 1)) * ...
  76                                 0.5 * s.dz_at_u(ix - 1, jx, kx) / (s.dz_at_u(ix - 1, jx, kx + 1) + s.dz_at_u(ix - 1, jx, kx));
  77                             I(M_entry) = row_offset + (ix - 1) + (jx - 1) * nx + kx * nx * ny;
  78                             J(M_entry) = M_col;
  79                             M_entry = M_entry + 1;
  80                         end
  81                         if not(is_last_row)
  82                             M(M_entry) = -s.area_w(ix, jx, kx + 1) * ...
  83                                 0.5 * (dzdx(ix, jx, kx + 1) + dzdx(ix, jx + 1, kx + 1)) * ...
  84                                 0.5 * s.dz_at_u(ix, jx, kx) / (s.dz_at_u(ix, jx, kx + 1) + s.dz_at_u(ix, jx, kx));
  85                             I(M_entry) = row_offset + ix + (jx - 1) * nx + kx * nx * ny;
  86                             J(M_entry) = M_col;
  87                             M_entry = M_entry + 1;
  88                         end
  89
  90                         if not(is_low_level)
  91                             if not(is_first_row)
  92                                 M(M_entry) = -s.area_w(ix - 1, jx, kx) * ...
  93                                     0.5 * (dzdx(ix - 1, jx, kx) + dzdx(ix - 1, jx + 1, kx)) * ...
  94                                     0.5 * s.dz_at_u(ix - 1, jx, kx) / (s.dz_at_u(ix - 1, jx, kx - 1) + s.dz_at_u(ix - 1, jx, kx));
  95                                 I(M_entry) = row_offset + (ix - 1) + (jx - 1) * nx + (kx - 1) * nx * ny;
  96                                 J(M_entry) = M_col;
  97                                 M_entry = M_entry + 1;
  98                             end
  99                             if not(is_last_row)
 100                                 M(M_entry) = -s.area_w(ix, jx, kx) * ...
 101                                     0.5 * (dzdx(ix, jx, kx) + dzdx(ix, jx + 1, kx)) * ...
 102                                     0.5 * s.dz_at_u(ix, jx, kx) / (s.dz_at_u(ix, jx, kx - 1) + s.dz_at_u(ix, jx, kx));
 103                                 I(M_entry) = row_offset + ix + (jx - 1) * nx + (kx - 1) * nx * ny;
 104                                 J(M_entry) = M_col;
 105                                 M_entry = M_entry + 1;
 106                             end
 107                         end
 108                     else
 109                         if not(is_first_row)
 110                             M(M_entry) = -s.area_w(ix - 1, jx, kx) * ...
 111                                 0.5 * (dzdx(ix - 1, jx, kx) + dzdx(ix - 1, jx + 1, kx)) * ...
 112                                 0.5 * s.dz_at_u(ix - 1, jx, kx) / (s.dz_at_u(ix - 1, jx, kx) + s.dz_at_u(ix - 1, jx, kx - 1));
 113                             I(M_entry) = row_offset + (ix - 1) + (jx - 1) * nx + (kx - 1) * nx * ny;
 114                             J(M_entry) = M_col;
 115                             M_entry = M_entry + 1;
 116                             M(M_entry) = -s.area_w(ix - 1, jx, kx + 1) * ...
 117                                 0.25 * (dzdx(ix - 1, jx, kx + 1) + dzdx(ix - 1, jx + 1, kx + 1));
 118                             I(M_entry) = row_offset + (ix - 1) + (jx - 1) * nx + kx * nx * ny;
 119                             J(M_entry) = M_col;
 120                             M_entry = M_entry + 1;
 121                         end
 122                         if not(is_last_row)
 123                             M(M_entry) = -s.area_w(ix, jx, kx) * ...
 124                                 0.5 * (dzdx(ix, jx, kx) + dzdx(ix, jx + 1, kx)) * ...
 125                                 0.5 * (s.dz_at_u(ix, jx, kx) / (s.dz_at_u(ix, jx, kx - 1) + s.dz_at_u(ix, jx, kx)));
 126                             I(M_entry) = row_offset + ix + (jx - 1) * nx + (kx - 1) * nx * ny;
 127                             J(M_entry) = M_col;
 128                             M_entry = M_entry + 1;
 129                             M(M_entry) = -s.area_w(ix, jx, kx + 1) * ...
 130                                 0.25 * (dzdx(ix, jx, kx + 1) + dzdx(ix, jx + 1, kx + 1));
 131                             I(M_entry) = row_offset + ix + (jx - 1) * nx + kx * nx * ny;
 132                             J(M_entry) = M_col;
 133                             M_entry = M_entry + 1;
 134                         end
 135                     end
 136                 elseif ux == 2
 137                     % V elements
 138                     M(M_entry) = s.area_v(ix,jx,kx);
 139                     I(M_entry) = M_row;
 140                     J(M_entry) = M_col;
 141                     M_entry = M_entry + 1;
 142
 143                     % Edge cases
 144                     is_top_level = (kx == Z_len);
 145                     is_low_level = (kx == 1);
 146                     is_first_row = (jx == 1);
 147                     is_last_row = (jx == Y_len);
 148
 149                     if not(is_top_level)
 150                         if not(is_first_row)
 151                             M(M_entry) = -s.area_w(ix, jx - 1, kx + 1) * ...
 152                                 0.5 * (dzdy(ix, jx - 1, kx + 1) + dzdy(ix + 1, jx - 1, kx + 1)) * ...
 153                                 0.5 * s.dz_at_v(ix, jx - 1, kx) / (s.dz_at_v(ix, jx - 1, kx + 1) + s.dz_at_v(ix, jx - 1, kx));
 154                             I(M_entry) = row_offset + ix + (jx - 2) * nx + kx * nx * ny;
 155                             J(M_entry) = M_col;
 156                             M_entry = M_entry + 1;
 157                         end
 158                         if not(is_last_row)
 159                             M(M_entry) = -s.area_w(ix, jx, kx + 1) * ...
 160                                 0.5 * (dzdy(ix, jx, kx + 1) + dzdy(ix + 1, jx, kx + 1)) * ...
 161                                 0.5 * s.dz_at_v(ix, jx, kx) / (s.dz_at_v(ix, jx, kx + 1) + s.dz_at_v(ix, jx, kx));
 162                             I(M_entry) = row_offset + ix + (jx - 1) * nx + kx * nx * ny;
 163                             J(M_entry) = M_col;
 164                             M_entry = M_entry + 1;
 165                         end
 166
 167                         if not(is_low_level)
 168                             if not(is_first_row)
 169                                 M(M_entry) = -s.area_w(ix, jx - 1, kx) * ...
 170                                     0.5 * (dzdy(ix, jx - 1, kx) + dzdy(ix + 1, jx - 1, kx)) * ...
 171                                     0.5 * s.dz_at_v(ix, jx - 1, kx) / (s.dz_at_v(ix, jx - 1, kx - 1) + s.dz_at_v(ix, jx - 1, kx));
 172                                 I(M_entry) = row_offset + ix + (jx - 2) * nx + (kx - 1) * nx * ny;
 173                                 J(M_entry) = M_col;
 174                                 M_entry = M_entry + 1;
 175                             end
 176                             if not(is_last_row)
 177                                 M(M_entry) = -s.area_w(ix, jx, kx) * ...
 178                                     0.5 * (dzdy(ix, jx, kx) + dzdy(ix + 1, jx, kx)) * ...
 179                                     0.5 * s.dz_at_v(ix, jx, kx) / (s.dz_at_v(ix, jx, kx - 1) + s.dz_at_v(ix, jx, kx));
 180                                 I(M_entry) = row_offset + ix + (jx - 1) * nx + (kx - 1) * nx * ny;
 181                                 J(M_entry) = M_col;
 182                                 M_entry = M_entry + 1;
 183                             end
 184                         end
 185                     else
 186                         if not(is_first_row)
 187                             M(M_entry) = -s.area_w(ix, jx - 1, kx) * ...
 188                                 0.5 * (dzdy(ix, jx - 1, kx) + dzdy(ix + 1, jx - 1, kx)) * ...
 189                                 0.5 * s.dz_at_v(ix, jx - 1, kx) / (s.dz_at_v(ix, jx - 1, kx) + s.dz_at_v(ix, jx - 1, kx - 1));
 190                             I(M_entry) = row_offset + ix + (jx - 2) * nx + (kx - 1) * nx * ny;
 191                             J(M_entry) = M_col;
 192                             M_entry = M_entry + 1;
 193
 194                             M(M_entry) = -s.area_w(ix, jx - 1, kx + 1) * ...
 195                                 0.25 * (dzdy(ix, jx - 1, kx + 1) + dzdy(ix + 1, jx - 1, kx + 1));
 196                             I(M_entry) = row_offset + ix + (jx - 2) * nx + kx * nx * ny;
 197                             J(M_entry) = M_col;
 198                             M_entry = M_entry + 1;
 199                         end
 200                         if not(is_last_row)
 201                             M(M_entry) = -s.area_w(ix, jx, kx) * ...
 202                                 0.5 * (dzdy(ix, jx, kx) + dzdy(ix + 1, jx, kx)) * ...
 203                                 0.5 * (s.dz_at_v(ix, jx, kx) / (s.dz_at_v(ix, jx, kx - 1) + s.dz_at_v(ix, jx, kx)));
 204                             I(M_entry) = row_offset + ix + (jx - 1) * nx + (kx - 1) * nx * ny;
 205                             J(M_entry) = M_col;
 206                             M_entry = M_entry + 1;
 207
 208                             M(M_entry) = -s.area_w(ix, jx, kx + 1) * ...
 209                                 0.25 * (dzdy(ix, jx, kx + 1) + dzdy(ix + 1, jx, kx + 1));
 210                             I(M_entry) = row_offset + ix + (jx - 1) * nx + kx * nx * ny;
 211                             J(M_entry) = M_col;
 212                             M_entry = M_entry + 1;
 213                         end
 214                     end
 215                 elseif ux == 3 && kx ~= 1
 216                     % W elements
 217                     M(M_entry) = s.area_w(ix, jx, kx);
 218                     I(M_entry) = M_row;
 219                     J(M_entry) = M_col;
 220                     M_entry = M_entry + 1;
 221                 end
 222             M_col = M_col + 1;
 223             M_row = M_row + 1;
 224             end
 225         end
 226     end
 227 end
 228 Mmat = sparse(I,J,M,n,n);