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                 fprintf('femwind_fortran and femwind_solve done\n')
 126                 err=max(err_compare);
 127                 fprintf('relative max error matlab vs fortran %g\n',err)
 128                 if err>1e-4,
 129                     warning('error too large')
 130                 end
 131         elseif params.run_fortran
 132                 [W,rate(sc,sc2)]=femwind_fortran(A,X,U0,params);
 133                 fprintf('femwind_fortran done ')
 134         else
 135                 [W,rate(sc,sc2)]=femwind_solve(A,X,U0,params);
 136                 fprintf('femwind_solve done ')
 137         end
 138         fprintf('rate(%i,%i)=%g\n',sc,sc2,rate(sc,sc2))
 139
 140         if params.graphics>1
 141             disp('graphics: solution')
 142
 143             % plot resulting wind
 144             figure(5),crate(sc,sc2)
 145             plot_mesh_3d(X,[1,nel(1)+1,1,nel(2)+1,1,1]), hold on,
 146             plot_wind_3d(CX,W)
 147             hold off
 148             axis equal
 149             title(['Final wind a=',string_diag_A])
 150
 151             figure(6),clf
 152             plot_mesh_3d(X,[1,nel(1)+1,1,nel(2)+1,1,1]), hold on,
 153             plot_wind_3d(CX,W,1)
 154             hold off
 155             axis equal
 156             title(['Final wind lowest layer a=',string_diag_A])
 157
 158             figure(7),clf
 159             plot_mesh_3d(X,[1,nel(1)+1,1,nel(2)+1,1,1]), hold on,
 160             plot_wind_3d(CX,W,1:2)
 161             hold off
 162             axis equal
 163             title(['Final wind lowest layers a=',string_diag_A])
 164
 165         end
 166
 167         if params.graphics>0
 168             disp('graphics: wind_at_h')
 169
 170             figure(8),clf
 171             height=10;
 172             bbox=[min(X{1}(:)),max(X{1}(:)),...
 173                 min(X{2}(:)),max(X{2}(:)),...
 174                 height,height];
 175             [XH,WH]=wind_at_h(X,CX,W,[20,20,1],bbox);
 176             plot_wind_3d(XH,WH)
 177             hold on
 178             plot_mesh_3d(X,[1,nel(1)+1,1,nel(2)+1,1,1])
 179             hold off
 180             axis equal
 181             title(['Final wind with a=',string_diag_A,' at ',num2str(height),' above terrain'])
 182         end
 183         if params.graphics > 2
 184             figure(9),clf
 185             plot_mesh_3d(X,[1,nel(1)+1,1,nel(2)+1,2,2])
 186             hold on
 187             wind_streamlines(X, CX, W, params)
 188             hold off
 189         end
 190         params.rate=rate;
 191         if params.save_files>1,
 192             sfile=[params.save_file_prefix,'_workspace.mat'];
 193             disp(['saving workspace to ',sfile])
 194             save(sfile,'-v7.3')
 195         end
 196         if params.save_files>0,
 197             sfile=[params.save_file_prefix,'_params.mat'];
 198             disp(['saving params to ',sfile])
 199             save(sfile,'params','-v7.3')
 200         end
 201     end
 202 end