femwind/femwind_main.m

   1 function params=femwind_main(params)
   2 % usage:
   3 % params=femwind_main    output params structure and exit
   4 % femwind_test(params)   run with params specified and copy to output
   5 %
   6 %           *** WARNING ***
   7 % DO NOT EDIT THIS FILE to set params or add outputs.
   8 % Instead, call the function to get the params, add or delete any fields you need, and call it again.
   9 % Debugging outputs can be added to the params structure as additional field.
  10 %
  11 %   *** all_test must pass fter any edits and before any merge ***
  12 %
  13 %           *** WARNING ***
  14
  15
  16 if ~exist('params','var') | isempty(params)
  17     params = params_defaults
  18     return
  19 end
  20
  21 if params.save_files >= 0
  22     diary(['femwind_',params.save_file_prefix,'_diary.txt'])
  23 end
  24 disp('femwind_main')
  25 format compact
  26
  27 if isfield(params,'mesh_top')
  28     if params.mesh_top>0
  29         disp('given params.mesh_top>0, computing params.expand')
  30         nz = params.nelem3(3); % elements in the vertical direction
  31         a = params.mesh_top/params.h(3); % desired height as multiple of first layer
  32         params.expand = findq(a,nz);
  33     end
  34
  35 end
  36
  37 params
  38
  39
  40 for sc2 = params.sc2_all
  41     for sc = params.sc_all
  42         sc,sc2
  43
  44         disp('setting up test case')
  45
  46         nel = sc*params.nelem3;  % elements in the 3 directions
  47         nel(1:2)=nel(1:2)*sc2
  48         h = params.h/sc;
  49         fprintf('mesh of %ix%ix%i cells\n',nel(1),nel(2),nel(3))
  50         params.id=sprintf('%ix%ix%i',nel); % to pass around
  51         string_diag_A=sprintf('%g %g %g',params.a); % for figure titles
  52         A = diag(params.a.^2);
  53
  54         % creating the grid
  55         expand=params.expand;
  56         X = regular_mesh(nel,h,params.expand^(1/sc));
  57         X = add_terrain_to_mesh(X,...
  58             params.terrain_shape,params.terrain_top,params.terrain_height);
  59         [CX,CH] = center_mesh(X); % get midpoints of elements
  60
  61         % initial wind at the centers of the elements
  62         rng(1);
  63         switch params.initial_wind
  64             case 'uniform'
  65                 disp('initial wind uniform in x direction')
  66                 U0={ones(nel),zeros(nel),zeros(nel)};
  67             case 'random-z'
  68                 % to test iterative methods with non-smooth initial error
  69                 disp('initial wind uniform in x direction random in z direction')
  70                 U0={ones(nel),zeros(nel),randn(nel)};
  71             case 'random-xz'
  72                 % to test iterative methods with non-smooth initial error
  73                 disp('initial wind uniform in x direction and z direction')
  74                 U0={ones(nel),zeros(nel),ones(nel)};
  75             case 'log'
  76                 disp('initial wind log profile in x direction')
  77                 U0={log(max(1,CH/params.roughness_height)),zeros(nel),zeros(nel)};
  78             otherwise
  79                 error(['unknown initial wind ',params.initial_wind])
  80         end
  81         if params.graphics>0
  82             disp('graphics: problem setup')
  83             % show mesh
  84             figure(1),clf
  85             plot_mesh_3d(X)
  86             axis equal
  87             title('The wind mesh, wind vector in centers, lambda in corners')
  88         end
  89
  90         if params.graphics>0
  91             % show initial wind
  92             figure(2),clf
  93             plot_mesh_3d(X,[1,nel(1)+1,1,nel(2)+1,1,1]), hold on,
  94             plot_wind_3d(CX,U0)
  95             hold off
  96             axis equal
  97             title('Initial wind')
  98         end
  99
 100         if params.graphics>1
 101             % show initial wind
 102             figure(3),clf
 103             plot_mesh_3d(X,[1,nel(1),1,nel(2)+1,1,1]), hold on,
 104             plot_wind_3d(CX,U0,1)
 105             hold off
 106             axis equal
 107             title('Initial wind lowest layer')
 108         end
 109
 110 %         Plot initial streamlines
 111         if params.graphics > 2
 112             figure(4), clf
 113             plot_mesh_3d(X,[1,nel(1)+1,1,nel(2)+1,2,2])
 114             hold on
 115             wind_streamlines(X, CX, U0, params);
 116             hold off
 117         end
 118         diary; diary
 119
 120         % solve
 121         if params.run_fortran && params.run_matlab
 122                 [W,rate(sc,sc2)]=femwind_fortran(A,X,U0,params);
 123                 [Wm,rate(sc,sc2)]=femwind_solve(A,X,U0,params);
 124                 for i=1:3, err_compare(i)=big(W{i}-Wm{i})/big(W{i});end
 125                 err=max(err_compare);
 126                 fprintf('relative max error matlab vs fortran %g\n',err)
 127                 if err>1e-4,
 128                     warning('error too large')
 129                 end
 130         elseif params.run_fortran
 131                 [W,rate(sc,sc2)]=femwind_fortran(A,X,U0,params);
 132         else
 133                 [W,rate(sc,sc2)]=femwind_solve(A,X,U0,params);
 134         end
 135
 136         if params.graphics>1
 137             disp('graphics: solution')
 138
 139             % plot resulting wind
 140             figure(5),crate(sc,sc2)
 141             plot_mesh_3d(X,[1,nel(1)+1,1,nel(2)+1,1,1]), hold on,
 142             plot_wind_3d(CX,W)
 143             hold off
 144             axis equal
 145             title(['Final wind a=',string_diag_A])
 146
 147             figure(6),clf
 148             plot_mesh_3d(X,[1,nel(1)+1,1,nel(2)+1,1,1]), hold on,
 149             plot_wind_3d(CX,W,1)
 150             hold off
 151             axis equal
 152             title(['Final wind lowest layer a=',string_diag_A])
 153
 154             figure(7),clf
 155             plot_mesh_3d(X,[1,nel(1)+1,1,nel(2)+1,1,1]), hold on,
 156             plot_wind_3d(CX,W,1:2)
 157             hold off
 158             axis equal
 159             title(['Final wind lowest layers a=',string_diag_A])
 160
 161         end
 162
 163         if params.graphics>0
 164             disp('graphics: wind_at_h')
 165
 166             figure(8),clf
 167             height=10;
 168             bbox=[min(X{1}(:)),max(X{1}(:)),...
 169                 min(X{2}(:)),max(X{2}(:)),...
 170                 height,height];
 171             [XH,WH]=wind_at_h(X,CX,W,[20,20,1],bbox);
 172             plot_wind_3d(XH,WH)
 173             hold on
 174             plot_mesh_3d(X,[1,nel(1)+1,1,nel(2)+1,1,1])
 175             hold off
 176             axis equal
 177             title(['Final wind with a=',string_diag_A,' at ',num2str(height),' above terrain'])
 178         end
 179         if params.graphics > 2
 180             figure(9),clf
 181             plot_mesh_3d(X,[1,nel(1)+1,1,nel(2)+1,2,2])
 182             hold on
 183             wind_streamlines(X, CX, W, params)
 184             hold off
 185         end
 186         params.rate=rate;
 187         if params.save_files>1,
 188             sfile=[params.save_file_prefix,'_workspace.mat'];
 189             disp(['saving workspace to ',sfile])
 190             save(sfile,'-v7.3')
 191         end
 192         if params.save_files>0,
 193             sfile=[params.save_file_prefix,'_params.mat'];
 194             disp(['saving params to ',sfile])
 195             save(sfile,'params','-v7.3')
 196         end
 197     end
 198 end