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.run_fortran=1;
  18     params.run_matlab=1;
  19     params.femwind_fortran_test=0;
  20     params.test_fortran=0;
  21     params.graphics=2;  % 1=basic, 2=all
  22     params.expand=1.2;  % exponential grid expansion in the vertical
  23     params.mesh_top=1000; % if given, ignore params_expand
  24     params.sc_all=[1,2]; % mesh refinements for tests at multiple scales
  25     params.sc2_all=[1,2,4];  % additional factors for horizonal mesh extent
  26     params.nelem3=[50,50,10]; % base size in elements, horizontal=2*odd
  27     params.h=[30,30,2]; % base mesh spacing before scaling
  28     params.a=[1 1 1]; % penalty factors in x y z directions
  29     params.initial_wind='log'; % or uniform
  30     params.roughness_height=0.5;
  31     params.terrain_shape='hill'; % terrain for add_terrain_to_mesh
  32     params.terrain_top='shift'; % mesh top treatment for add_terrain_to_mesh
  33     params.terrain_height=0.1; % terrain height as part of domain height
  34     params.solver='2-level' ; % see sparse_solve.m
  35     params.maxit=50; % max iterations
  36     params.maxit_coarse=8; % 2 smoothing, coarse, 2 smoothing, coarse, 2 smoothing
  37     params.coarsest_iter=100; % 0 = direct solver n>0 number of iterations
  38     params.nsmooth=3; % smoothing iterations before correcton
  39     params.coarsening='2 linear';
  40     % params.coarse_P='variational';
  41     params.coarse_K='assembly';
  42     params.P_by_x=1;  % prolongation by geometrically linear interpolation
  43     params.smoothing='vertical sweeps';
  44     % params.smoothing='3D red-black';
  45     params.restol=1e-6;
  46     params.exact=0; % compare with exact solution to compute error
  47     params.slice=0.5; % vertical y slice of error to display, 0 to 1
  48     params.err_slice_fig=12; % figure number for residual slice
  49     params.res_slice_fig=13; % figure number for error slice
  50     params.iterations_fig=14; % figure number for iterations progress
  51     params.maxaspect=3;  % do not coarsen vertically if vertical layer is too thick
  52     params.minaspect=1/3; % do not coarsen horizontally if the layer is too thin
  53     params.levels=15;
  54     params.apply_coarse_boundary_conditions=1;
  55     params.nsmooth_coarse=2;
  56     params.save_files=0; % save progress levels=3, workspace=2 params only=1
  57     params.save_file_prefix='femwind';
  58     %Define Streamline Starting Points: Defined in terms of scale*nelem
  59     params.in_height_stream = [175];
  60     params.time_stream  = 0;
  61     params.st_contour = 1; %Produce contour planes in streamlines plot 0 off, 1 on
  62     params.st_quiver = 1; %Produce vectors along streamlines 0 off, 1 on
  63
  64
  65     return
  66 end
  67
  68 if params.save_files >= 0
  69     diary(['femwind_',params.save_file_prefix,'_diary.txt'])
  70 end
  71 disp('femwind_main')
  72 format compact
  73
  74 if isfield(params,'mesh_top')
  75     if params.mesh_top>0
  76         disp('given params.mesh_top>0, computing params.expand')
  77         nz = params.nelem3(3); % elements in the vertical direction
  78         a = params.mesh_top/params.h(3); % desired height as multiple of first layer
  79         params.expand = findq(a,nz);
  80     end
  81
  82 end
  83
  84 params
  85
  86
  87 for sc2 = params.sc2_all
  88     for sc = params.sc_all
  89         sc,sc2
  90
  91         disp('setting up test case')
  92
  93         nel = sc*params.nelem3;  % elements in the 3 directions
  94         nel(1:2)=nel(1:2)*sc2
  95         h = params.h/sc;
  96         fprintf('mesh of %ix%ix%i cells\n',nel(1),nel(2),nel(3))
  97         params.id=sprintf('%ix%ix%i',nel); % to pass around
  98         string_diag_A=sprintf('%g %g %g',params.a); % for figure titles
  99         A = diag(params.a.^2);
 100
 101         % creating the grid
 102         expand=params.expand;
 103         X = regular_mesh(nel,h,params.expand^(1/sc));
 104         X = add_terrain_to_mesh(X,...
 105             params.terrain_shape,params.terrain_top,params.terrain_height);
 106         [CX,CH] = center_mesh(X); % get midpoints of elements
 107
 108         % initial wind at the centers of the elements
 109         rng(1);
 110         switch params.initial_wind
 111             case 'uniform'
 112                 disp('initial wind uniform in x direction')
 113                 U0={ones(nel),zeros(nel),zeros(nel)};
 114             case 'random-z'
 115                 % to test iterative methods with non-smooth initial error
 116                 disp('initial wind uniform in x direction random in z direction')
 117                 U0={ones(nel),zeros(nel),randn(nel)};
 118             case 'random-xz'
 119                 % to test iterative methods with non-smooth initial error
 120                 disp('initial wind uniform in x direction and z direction')
 121                 U0={ones(nel),zeros(nel),ones(nel)};
 122             case 'log'
 123                 disp('initial wind log profile in x direction')
 124                 U0={log(max(1,CH/params.roughness_height)),zeros(nel),zeros(nel)};
 125             otherwise
 126                 error(['unknown initial wind ',params.initial_wind])
 127         end
 128         if params.graphics>0
 129             disp('graphics: problem setup')
 130             % show mesh
 131             figure(1),clf
 132             plot_mesh_3d(X)
 133             axis equal
 134             title('The wind mesh, wind vector in centers, lambda in corners')
 135         end
 136
 137         if params.graphics>0
 138             % show initial wind
 139             figure(2),clf
 140             plot_mesh_3d(X,[1,nel(1)+1,1,nel(2)+1,1,1]), hold on,
 141             plot_wind_3d(CX,U0)
 142             hold off
 143             axis equal
 144             title('Initial wind')
 145         end
 146
 147         if params.graphics>1
 148             % show initial wind
 149             figure(3),clf
 150             plot_mesh_3d(X,[1,nel(1),1,nel(2)+1,1,1]), hold on,
 151             plot_wind_3d(CX,U0,1)
 152             hold off
 153             axis equal
 154             title('Initial wind lowest layer')
 155         end
 156
 157 %         Plot initial streamlines
 158         if params.graphics > 2
 159             figure(4), clf
 160             plot_mesh_3d(X,[1,nel(1)+1,1,nel(2)+1,2,2])
 161             hold on
 162             wind_streamlines(X, CX, U0, params);
 163             hold off
 164         end
 165         diary; diary
 166
 167         % solve
 168         [W,rate(sc,sc2)]=femwind_fortran(A,X,U0,params);
 169
 170         if params.graphics>1
 171             disp('graphics: solution')
 172
 173             % plot resulting wind
 174             figure(5),crate(sc,sc2)
 175             plot_mesh_3d(X,[1,nel(1)+1,1,nel(2)+1,1,1]), hold on,
 176             plot_wind_3d(CX,W)
 177             hold off
 178             axis equal
 179             title(['Final wind a=',string_diag_A])
 180
 181             figure(6),clf
 182             plot_mesh_3d(X,[1,nel(1)+1,1,nel(2)+1,1,1]), hold on,
 183             plot_wind_3d(CX,W,1)
 184             hold off
 185             axis equal
 186             title(['Final wind lowest layer a=',string_diag_A])
 187
 188             figure(7),clf
 189             plot_mesh_3d(X,[1,nel(1)+1,1,nel(2)+1,1,1]), hold on,
 190             plot_wind_3d(CX,W,1:2)
 191             hold off
 192             axis equal
 193             title(['Final wind lowest layers a=',string_diag_A])
 194
 195         end
 196
 197         if params.graphics>0
 198             disp('graphics: wind_at_h')
 199
 200             figure(8),clf
 201             height=10;
 202             bbox=[min(X{1}(:)),max(X{1}(:)),...
 203                 min(X{2}(:)),max(X{2}(:)),...
 204                 height,height];
 205             [XH,WH]=wind_at_h(X,CX,W,[20,20,1],bbox);
 206             plot_wind_3d(XH,WH)
 207             hold on
 208             plot_mesh_3d(X,[1,nel(1)+1,1,nel(2)+1,1,1])
 209             hold off
 210             axis equal
 211             title(['Final wind with a=',string_diag_A,' at ',num2str(height),' above terrain'])
 212         end
 213         if params.graphics > 2
 214             figure(9),clf
 215             plot_mesh_3d(X,[1,nel(1)+1,1,nel(2)+1,2,2])
 216             hold on
 217             wind_streamlines(X, CX, W, params)
 218             hold off
 219         end
 220         params.rate=rate;
 221         if params.save_files>1,
 222             sfile=[params.save_file_prefix,'_workspace.mat'];
 223             disp(['saving workspace to ',sfile])
 224             save(sfile,'-v7.3')
 225         end
 226         if params.save_files>0,
 227             sfile=[params.save_file_prefix,'_params.mat'];
 228             disp(['saving params to ',sfile])
 229             save(sfile,'params','-v7.3')
 230         end
 231     end
 232 end