quicwind/Mmul_v.m

   1 function res = Mmul_v(X,trans,v)
   2 % multiplies the transformation of variable flux matrix
   3 % (or its transpose) and a vector
   4 % the operation is res = op(M) * v, where op(M) = M or M^T
   5 % set:
   6 %    X   mesh
   7 %    t   indicates whether the matrix is transposed or not:
   8 %           trans = 'T' or trans = 't': op(M) = M^T
   9 %           trans = 'N' or trans = 'n': op(M) = M
  10 %        default is 'N'
  11 %    v   vector to multiply
  12
  13 % Error checking
  14 check_mesh(X);
  15 if ~exist('trans','var')
  16     trans = 'N';
  17 end
  18 if (trans ~= 'T' && trans ~= 't' && trans ~= 'N' && trans ~= 'n')
  19     error('trans must be N, T, n, or t')
  20 end
  21
  22 % Initialize result
  23 res = zeros(size(v));
  24
  25 % Mesh sizes and vars
  26 nx = size(X{1},1)-1;
  27 ny = size(X{1},2)-1;
  28 nz = size(X{1},3)-1;
  29 x = X{1}; y = X{2}; z = X{3};
  30
  31 % This needs to be removed and replaced in code
  32 s = cell_sizes(X);
  33 % continue one layer up
  34 s.dz_at_u(:,:,nz+1)=s.dz_at_u(:,:,nz);
  35 s.dz_at_v(:,:,nz+1)=s.dz_at_v(:,:,nz);
  36
  37 % slope in x direction
  38 dzdx = zeros(nx,ny+1,nz+1);
  39 for k=1:nz+1
  40     for j=1:ny+1
  41         for i=1:nx
  42             dzdx(i,j,k) = ((z(i+1,j,k)-z(i,j,k))/(x(i+1,j,k)-x(i,j,k)));
  43         end
  44     end
  45 end
  46
  47 % slope in y direction
  48 dzdy = zeros(nx+1,ny,nz+1);
  49 for k=1:nz+1
  50     for j=1:ny
  51         for i=1:nx+1
  52             dzdy(i,j,k) = ((z(i,j+1,k)-z(i,j,k))/(y(i,j+1,k)-y(i,j,k)));
  53         end
  54     end
  55 end
  56
  57 if (trans == 't' || trans == 'T')
  58
  59 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% M^T * v %%%%%%%%%%%%%%%%%%%%%%%%%%%%
  60 z_factor = (nx+1)*ny*nz + nx*(ny+1)*nz;
  61 for k=1:nz
  62     for j=1:ny
  63         for i=1:(nx+1)
  64             area = s.area_u(i,j,k);
  65             res(i + (j-1) * (nx+1) + (k-1) * ny * (nx+1)) = ...
  66                 area * v(i + (j-1) * (nx+1) + (k-1) * ny * (nx+1));
  67             is_low_level = (k == 1);
  68             is_top_level = ((k+1) == (nz+1));
  69             is_first_row = (i == 1);
  70             is_last_row = (i == (nx + 1));
  71             if (is_low_level)
  72                 % Upper level
  73                 if (~is_last_row)
  74                     area_w = s.area_w(i,j,k+1);
  75                     dzdx_at_k = 0.5*(dzdx(i,j,k+1) + dzdx(i,j+1,k+1));
  76                     mul_fact = 0.5*(s.dz_at_u(i,j,k) / (s.dz_at_u(i,j,k) + s.dz_at_u(i,j,k+1)));
  77                     res(i + (j-1) * (nx+1) + (k-1) * ny * (nx+1)) = ...
  78                         res(i + (j-1) * (nx+1) + (k-1) * ny * (nx+1)) - ...
  79                         area_w * dzdx_at_k * mul_fact * v(z_factor + i + (j-1)*nx + k*nx*ny);
  80                 end
  81                 if (~is_first_row)
  82                     area_w = s.area_w(i-1,j,k+1);
  83                     dzdx_at_k = 0.5*(dzdx(i-1,j,k+1) + dzdx(i-1,j+1,k+1));
  84                     mul_fact = 0.5*(s.dz_at_u(i-1,j,k) / (s.dz_at_u(i-1,j,k) + s.dz_at_u(i-1,j,k+1)));
  85                     res(i + (j-1) * (nx+1) + (k-1) * ny * (nx+1)) = ...
  86                         res(i + (j-1) * (nx+1) + (k-1) * ny * (nx+1)) - ...
  87                         area_w * dzdx_at_k * mul_fact * v(z_factor + i - 1 + (j-1)*nx + k*nx*ny);
  88                 end
  89             elseif (is_top_level)
  90
  91                 % Current level
  92                 if (~is_last_row)
  93                     area_w = s.area_w(i,j,k);
  94                     dzdx_at_k = 0.5*(dzdx(i,j,k) + dzdx(i,j+1,k));
  95                     mul_fact = 0.5*(s.dz_at_u(i,j,k) / (s.dz_at_u(i,j,k) + s.dz_at_u(i,j,k-1)));
  96                     res(i + (j-1) * (nx+1) + (k-1) * ny * (nx+1)) = ...
  97                         res(i + (j-1) * (nx+1) + (k-1) * ny * (nx+1)) - ...
  98                         area_w * dzdx_at_k * mul_fact * v(z_factor + i + (j-1)*nx + (k-1)*nx*ny);
  99                 end
 100                 if (~is_first_row)
 101                     area_w = s.area_w(i-1,j,k);
 102                     dzdx_at_k = 0.5*(dzdx(i-1,j,k) + dzdx(i-1,j+1,k));
 103                     mul_fact = 0.5*(s.dz_at_u(i-1,j,k) / (s.dz_at_u(i-1,j,k) + s.dz_at_u(i-1,j,k-1)));
 104                     res(i + (j-1) * (nx+1) + (k-1) * ny * (nx+1)) = ...
 105                         res(i + (j-1) * (nx+1) + (k-1) * ny * (nx+1)) - ...
 106                         area_w * dzdx_at_k * mul_fact * v(z_factor + i - 1 + (j-1)*nx + (k-1)*nx*ny);
 107                 end
 108
 109                 % Upper and current level
 110                 if (~is_last_row)
 111                     area_w = s.area_w(i,j,k+1);
 112                     dzdx_at_k = 0.5*(dzdx(i,j,k+1) + dzdx(i,j+1,k+1));
 113                     mul_fact = s.dz_at_u(i,j,k) / (s.dz_at_u(i,j,k) + s.dz_at_u(i,j,k));
 114                     res(i + (j-1) * (nx+1) + (k-1) * ny * (nx+1)) = ...
 115                         res(i + (j-1) * (nx+1) + (k-1) * ny * (nx+1)) - ...
 116                         area_w * dzdx_at_k * mul_fact * v(z_factor + i + (j-1)*nx + k*nx*ny);
 117                 end
 118                 if (~is_first_row)
 119                     area_w = s.area_w(i-1,j,k+1);
 120                     dzdx_at_k = 0.5*(dzdx(i-1,j,k+1) + dzdx(i-1,j+1,k+1));
 121                     mul_fact = s.dz_at_u(i-1,j,k) / (s.dz_at_u(i-1,j,k) + s.dz_at_u(i-1,j,k));
 122                     res(i + (j-1) * (nx+1) + (k-1) * ny * (nx+1)) = ...
 123                         res(i + (j-1) * (nx+1) + (k-1) * ny * (nx+1)) - ...
 124                         area_w * dzdx_at_k * mul_fact * v(z_factor + i - 1 + (j-1)*nx + k*nx*ny);
 125                 end
 126             else
 127                 % Upper level
 128                 if (~is_last_row)
 129                     area_w = s.area_w(i,j,k);
 130                     dzdx_at_k = 0.5*(dzdx(i,j,k) + dzdx(i,j+1,k));
 131                     mul_fact = 0.5*(s.dz_at_u(i,j,k) / (s.dz_at_u(i,j,k) + s.dz_at_u(i,j,k-1)));
 132                     res(i + (j-1) * (nx+1) + (k-1) * ny * (nx+1)) = ...
 133                         res(i + (j-1) * (nx+1) + (k-1) * ny * (nx+1)) - ...
 134                         area_w * dzdx_at_k * mul_fact * v(z_factor + i + (j-1)*nx + (k-1)*nx*ny);
 135                 end
 136                 if (~is_first_row)
 137                     area_w = s.area_w(i-1,j,k);
 138                     dzdx_at_k = 0.5*(dzdx(i-1,j,k) + dzdx(i-1,j+1,k));
 139                     mul_fact = 0.5*(s.dz_at_u(i-1,j,k) / (s.dz_at_u(i-1,j,k) + s.dz_at_u(i-1,j,k-1)));
 140                     res(i + (j-1) * (nx+1) + (k-1) * ny * (nx+1)) = ...
 141                         res(i + (j-1) * (nx+1) + (k-1) * ny * (nx+1)) - ...
 142                         area_w * dzdx_at_k * mul_fact * v(z_factor + i - 1 + (j-1)*nx + (k-1)*nx*ny);
 143                 end
 144
 145                 % Current level
 146                 if (~is_last_row)
 147                     area_w = s.area_w(i,j,k+1);
 148                     dzdx_at_k = 0.5*(dzdx(i,j,k+1) + dzdx(i,j+1,k+1));
 149                     mul_fact = 0.5 * s.dz_at_u(i,j,k) / (s.dz_at_u(i,j,k) + s.dz_at_u(i,j,k+1));
 150                     res(i + (j-1) * (nx+1) + (k-1) * ny * (nx+1)) = ...
 151                         res(i + (j-1) * (nx+1) + (k-1) * ny * (nx+1)) - ...
 152                         area_w * dzdx_at_k * mul_fact * v(z_factor + i + (j-1)*nx + k*nx*ny);
 153                 end
 154                 if (~is_first_row)
 155                     area_w = s.area_w(i-1,j,k+1);
 156                     dzdx_at_k = 0.5*(dzdx(i-1,j,k+1) + dzdx(i-1,j+1,k+1));
 157                     mul_fact = 0.5 * (s.dz_at_u(i-1,j,k) / (s.dz_at_u(i-1,j,k) + s.dz_at_u(i-1,j,k+1)));
 158                     res(i + (j-1) * (nx+1) + (k-1) * ny * (nx+1)) = ...
 159                         res(i + (j-1) * (nx+1) + (k-1) * ny * (nx+1)) - ...
 160                         area_w * dzdx_at_k * mul_fact * v(z_factor + i - 1 + (j-1)*nx + k*nx*ny);
 161                 end
 162             end
 163         end
 164     end
 165 end
 166 add_factor = (nx+1)*ny*nz;
 167 for k=1:nz
 168     for j=1:(ny+1)
 169         for i=1:nx
 170             area = s.area_v(i,j,k);
 171             res(add_factor + i + (j-1) * nx + (k-1) * (ny+1) * nx) = ...
 172                 area * v(add_factor + i + (j-1) * nx + (k-1) * (ny+1) * nx);
 173             is_low_level = (k == 1);
 174             is_top_level = ((k+1) == (nz+1));
 175             is_first_col = (j == 1);
 176             is_last_col = (j == (ny + 1));
 177             if (is_low_level)
 178                 % Upper level
 179                 if (~is_last_col)
 180                     area_w = s.area_w(i,j,k+1);
 181                     dzdy_at_k = 0.5*(dzdy(i,j,k+1) + dzdy(i+1,j,k+1));
 182                     mul_fact = 0.5*(s.dz_at_v(i,j,k) / (s.dz_at_v(i,j,k) + s.dz_at_v(i,j,k+1)));
 183                     res(add_factor + i + (j-1) * nx + (k-1) * (ny+1) * nx) = ...
 184                         res(add_factor + i + (j-1) * nx + (k-1) * (ny+1) * nx) - ...
 185                         area_w * dzdy_at_k * mul_fact * v(z_factor + i + (j-1)*nx + k*nx*ny);
 186                 end
 187                 if (~is_first_col)
 188                     area_w = s.area_w(i,j-1,k+1);
 189                     dzdy_at_k = 0.5*(dzdy(i,j-1,k+1) + dzdy(i+1,j-1,k+1));
 190                     mul_fact = 0.5*(s.dz_at_v(i,j-1,k) / (s.dz_at_v(i,j-1,k) + s.dz_at_v(i,j-1,k+1)));
 191                     res(add_factor + i + (j-1) * nx + (k-1) * (ny+1) * nx) = ...
 192                         res(add_factor + i + (j-1) * nx + (k-1) * (ny + 1) * nx) - ...
 193                         area_w * dzdy_at_k * mul_fact * v(z_factor + i + (j-2)*nx + k*nx*ny);
 194                 end
 195             elseif (is_top_level)
 196
 197                 % Current level
 198                 if (~is_last_col)
 199                     area_w = s.area_w(i,j,k);
 200                     dzdy_at_k = 0.5*(dzdy(i,j,k) + dzdy(i+1,j,k));
 201                     mul_fact = 0.5*(s.dz_at_v(i,j,k) / (s.dz_at_v(i,j,k) + s.dz_at_v(i,j,k-1)));
 202                     res(add_factor + i + (j-1) * nx + (k-1) * (ny + 1) * nx) = ...
 203                         res(add_factor + i + (j-1) * nx + (k-1) * (ny + 1) * nx) - ...
 204                         area_w * dzdy_at_k * mul_fact * v(z_factor + i + (j-1)*nx + (k-1)*nx*ny);
 205                 end
 206                 if (~is_first_col)
 207                     area_w = s.area_w(i,j-1,k);
 208                     dzdy_at_k = 0.5*(dzdy(i,j-1,k) + dzdy(i+1,j-1,k));
 209                     mul_fact = 0.5*(s.dz_at_v(i,j-1,k) / (s.dz_at_v(i,j-1,k) + s.dz_at_v(i,j-1,k-1)));
 210                     res(add_factor + i + (j-1) * nx + (k-1) * (ny + 1) * nx) = ...
 211                         res(add_factor + i + (j-1) * nx + (k-1) * (ny + 1) * nx) - ...
 212                         area_w * dzdy_at_k * mul_fact * v(z_factor + i + (j-2)*nx + (k-1)*nx*ny);
 213                 end
 214
 215                 % Cuurent and upper level
 216                 if (~is_last_col)
 217                     area_w = s.area_w(i,j,k+1);
 218                     dzdy_at_k = 0.5*(dzdy(i,j,k+1) + dzdy(i+1,j,k+1));
 219                     mul_fact = s.dz_at_v(i,j,k) / (s.dz_at_v(i,j,k) + s.dz_at_v(i,j,k));
 220                     res(add_factor + i + (j-1) * nx + (k-1) * (ny + 1) * nx) = ...
 221                         res(add_factor + i + (j-1) * nx + (k-1) * (ny + 1) * nx) - ...
 222                         area_w * dzdy_at_k * mul_fact * v(z_factor + i + (j-1)*nx + k*nx*ny);
 223                 end
 224                 if (~is_first_col)
 225                     area_w = s.area_w(i,j-1,k+1);
 226                     dzdy_at_k = 0.5*(dzdy(i,j-1,k+1) + dzdy(i+1,j-1,k+1));
 227                     mul_fact = s.dz_at_u(i,j-1,k) / (s.dz_at_u(i,j-1,k) + s.dz_at_u(i,j-1,k));
 228                     res(add_factor + i + (j-1) * nx + (k-1) * (ny + 1) * nx) = ...
 229                         res(add_factor + i + (j-1) * nx + (k-1) * (ny + 1) * nx) - ...
 230                         area_w * dzdy_at_k * mul_fact * v(z_factor + i + (j-2)*nx + k*nx*ny);
 231                 end
 232             else
 233                 % Upper level
 234                 if (~is_last_col)
 235                     area_w = s.area_w(i,j,k);
 236                     dzdy_at_k = 0.5*(dzdy(i,j,k) + dzdy(i+1,j,k));
 237                     mul_fact = 0.5*(s.dz_at_v(i,j,k) / (s.dz_at_v(i,j,k) + s.dz_at_v(i,j,k-1)));
 238                     res(add_factor + i + (j-1) * nx + (k-1) * (ny + 1) * nx) = ...
 239                         res(add_factor + i + (j-1) * nx + (k-1) * (ny + 1) * nx) - ...
 240                         area_w * dzdy_at_k * mul_fact * v(z_factor + i + (j-1)*nx + (k-1)*nx*ny);
 241                 end
 242                 if (~is_first_col)
 243                     area_w = s.area_w(i,j-1,k);
 244                     dzdy_at_k = 0.5*(dzdy(i,j-1,k) + dzdy(i+1,j-1,k));
 245                     mul_fact = 0.5*(s.dz_at_v(i,j-1,k) / (s.dz_at_v(i,j-1,k) + s.dz_at_v(i,j-1,k-1)));
 246                     res(add_factor + i + (j-1) * nx + (k-1) * (ny + 1) * nx) = ...
 247                         res(add_factor + i + (j-1) * nx + (k-1) * (ny + 1) * nx) - ...
 248                         area_w * dzdy_at_k * mul_fact * v(z_factor + i + (j-2)*nx + (k-1)*nx*ny);
 249                 end
 250
 251                 % Current level
 252                 if (~is_last_col)
 253                     area_w = s.area_w(i,j,k+1);
 254                     dzdy_at_k = 0.5*(dzdy(i,j,k+1) + dzdy(i+1,j,k+1));
 255                     mul_fact = 0.5 * s.dz_at_v(i,j,k) / (s.dz_at_v(i,j,k) + s.dz_at_v(i,j,k+1));
 256                     res(add_factor + i + (j-1) * nx + (k-1) * (ny + 1) * nx) = ...
 257                         res(add_factor + i + (j-1) * nx + (k-1) * (ny + 1) * nx) - ...
 258                         area_w * dzdy_at_k * mul_fact * v(z_factor + i + (j-1)*nx + k*nx*ny);
 259                 end
 260                 if (~is_first_col)
 261                     area_w = s.area_w(i,j-1,k+1);
 262                     dzdy_at_k = 0.5*(dzdy(i,j-1,k+1) + dzdy(i+1,j-1,k+1));
 263                     mul_fact = 0.5 * (s.dz_at_v(i,j-1,k) / (s.dz_at_v(i,j-1,k) + s.dz_at_v(i,j-1,k+1)));
 264                     res(add_factor + i + (j-1) * nx + (k-1) * (ny + 1) * nx) = ...
 265                         res(add_factor + i + (j-1) * nx + (k-1) * (ny + 1) * nx) - ...
 266                         area_w * dzdy_at_k * mul_fact * v(z_factor + i + (j-2)*nx + k*nx*ny);
 267                 end
 268             end
 269         end
 270     end
 271 end
 272 add_factor = add_factor + nx*(ny+1)*nz;
 273 for k=1:(nz+1)
 274     for j=1:ny
 275         for i=1:nx
 276             is_low_level = (k == 1);
 277             % Zero on the ground
 278             if (is_low_level)
 279                 res(add_factor + i + (j-1)*nx + (k-1)*nx*ny) = 0;
 280             else
 281                 res(add_factor + i + (j-1)*nx + (k-1)*nx*ny) = ...
 282                     s.area_w(i,j,k) * v(add_factor + i + (j-1)*nx + (k-1)*nx*ny);
 283             end
 284         end
 285     end
 286 end
 287
 288 else % trans = 'N' || trans = 'n'
 289
 290 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% M * v %%%%%%%%%%%%%%%%%%%%%%%%%%%%
 291
 292 for k=1:nz
 293     for j=1:ny
 294         for i=1:(nx+1)
 295             area = s.area_u(i,j,k);
 296             res(i + (j-1) * (nx+1) + (k-1) * ny * (nx+1)) = ...
 297                 area * v(i + (j-1) * (nx+1) + (k-1) * ny * (nx+1));
 298         end
 299     end
 300 end
 301 add_factor = (nx+1)*ny*nz;
 302 for k=1:nz
 303     for j=1:(ny+1)
 304         for i=1:nx
 305             area = s.area_v(i,j,k);
 306             res(add_factor + i + (j-1) * nx + (k-1) * (ny+1) * nx) = ...
 307                 area * v(add_factor + i + (j-1) * nx + (k-1) * (ny+1) * nx);
 308         end
 309     end
 310 end
 311 add_factor_y = add_factor;
 312 add_factor = add_factor + nx*(ny+1)*nz;
 313 for k=1:(nz+1)
 314     for j=1:ny
 315         for i=1:nx
 316
 317             is_top_level = (k == (nz+1));
 318             is_low_level = (k == 1);
 319
 320             % Zero on the ground
 321             if (is_low_level)
 322                 res(add_factor + i + (j-1)*nx + (k-1)*nx*ny) = 0;
 323             else
 324                 area = s.area_w(i,j,k);
 325                 dzdx_at_k = 0.5*(dzdx(i,j,k) + dzdx(i,j+1,k));
 326                 dzdy_at_k = 0.5*(dzdy(i,j,k) + dzdy(i+1,j,k));
 327                 if (is_top_level)
 328                     u_at_k = 0.5 * ((v(i + (j-1)*(nx+1) + (k-2)*(nx+1)*ny) + v(i + 1 + (j-1)*(nx+1) + (k-2)*(nx+1)*ny)) * s.dz_at_u(i,j,k) + ...
 329                                     (v(i + (j-1)*(nx+1) + (k-2)*(nx+1)*ny) + v(i + 1 + (j-1)*(nx+1) + (k-2)*(nx+1)*ny)) * s.dz_at_u(i,j,k-1)) / ...
 330                                     (s.dz_at_u(i,j,k) + s.dz_at_u(i,j,k-1));
 331                     v_at_k = 0.5 * ((v(add_factor_y + i + (j-1)*nx + (k-2)*nx*(ny+1)) + v(add_factor_y + i + j*nx + (k-2)*nx*(ny+1))) * s.dz_at_v(i,j,k) + ...
 332                                     (v(add_factor_y + i + (j-1)*nx + (k-2)*nx*(ny+1)) + v(add_factor_y + i + j*nx + (k-2)*nx*(ny+1))) * s.dz_at_v(i,j,k-1)) / ...
 333                                     (s.dz_at_v(i,j,k) + s.dz_at_v(i,j,k-1));
 334                     res(add_factor + i + (j-1)*nx + (k-1)*nx*ny) = ...
 335                         area * (v(add_factor + i + (j-1)*nx + (k-1)*nx*ny) - ...
 336                         dzdx_at_k * u_at_k - dzdy_at_k * v_at_k);
 337                 else
 338                     u_at_k = 0.5 * ((v(i + (j-1)*(nx+1) + (k-1)*(nx+1)*ny) + v(i + 1 + (j-1)*(nx+1) + (k-1)*(nx+1)*ny)) * s.dz_at_u(i,j,k) + ...
 339                                     (v(i + (j-1)*(nx+1) + (k-2)*(nx+1)*ny) + v(i + 1 + (j-1)*(nx+1) + (k-2)*(nx+1)*ny)) * s.dz_at_u(i,j,k-1)) / ...
 340                                     (s.dz_at_u(i,j,k) + s.dz_at_u(i,j,k-1));
 341                     v_at_k = 0.5 * ((v(add_factor_y + i + (j-1)*nx + (k-1)*nx*(ny+1)) + v(add_factor_y + i + j*nx + (k-1)*nx*(ny+1))) * s.dz_at_v(i,j,k) + ...
 342                                     (v(add_factor_y + i + (j-1)*nx + (k-2)*nx*(ny+1)) + v(add_factor_y + i + j*nx + (k-2)*nx*(ny+1))) * s.dz_at_v(i,j,k-1)) / ...
 343                                     (s.dz_at_v(i,j,k) + s.dz_at_v(i,j,k-1));
 344                     res(add_factor + i + (j-1)*nx + (k-1)*nx*ny) = ...
 345                         area * (v(add_factor + i + (j-1)*nx + (k-1)*nx*ny) - ...
 346                         dzdx_at_k * u_at_k - dzdy_at_k * v_at_k);
 347                 end
 348             end
 349         end
 350     end
 351 end
 352 end
 353
 354
 355