180 likes | 314 Views
Linear solutions of flow over mountains COSMO General Meeting WG2 'Numerics and Dynamics' 15.09.2008 Michael Baldauf Deutscher Wetterdienst, Offenbach, Germany. Development of a program to calculate a reference solution for convergence tests (task 5) Linear solutions of flow over mountain
E N D
Linear solutions of flow over mountains COSMO General Meeting WG2 'Numerics and Dynamics' 15.09.2008 Michael Baldauf Deutscher Wetterdienst, Offenbach, Germany
Development of a program to calculate a reference solution for convergence tests (task 5) • Linear solutions of flow over mountain • there exist a huge number of different linear solutions for flow over mountains: • Queney (1947, 1948), Smith (1979, 1980), ... • reason: different goals: • to get insight into the physical process: try to approximate as far as possible to get simple solutions/formulae • to test numerical models:try to reduce the approximations as far as possible (perhaps by requiring more preconditions) to get an 'analytical' solution which fits best to the PDE-system (better say: a solution which can be calculated with higher numerical confidence, e.g. calculation of integrals or FFT's instead of solving PDE's)
Starting point: compressible Euler equations COSMO-dynamical core • Preconditions: • no friction • only adiabatic processes (in particular no phase changes) • ideal gas law • cp=const., cV=const., R=const. • no Coriolis force • all movements take place on a plane (no earth curvature) • These preconditions can be fulfilled easily by a dynamical core. • Only 2 approximations will be made: • linearisation 1/Fr = h N / U <<< 1 simulate very flat mountains at high inflow velocities • the assumption that kz=const (see below; not absolutely necessary)
Compressible Euler equations (adiabatic, f=0,...) e.g. Smith (1979) 'tracer'-parameter: 1 = 0/1 : hydrostatic / non-hydrostatic approximation 2 = 0/1 : incompressible / compressible model 3 = 0/1 : shallow / deep convection/atmosphere
Base state: stationary, hydrostatic, only depends on z
Fouriertransformation: horizontal and temporal Perturbation equations: 4th 'tracer'-parameter introduced (Pichler, 1997): 4=0/1 : small / big Mach-numbers
Some abbreviations: Def.: Heterogenity (Queney, 1947) Def: Stability parameter Def: Mach-number Def:
Stationary case (=0) From perturbation equations: express u', v', ' and p' by w' equation 2nd order for w'(kx, ky, z): with coefficient functions:
Variable transformation (modified Bretherton (1960)-transformation): The whole derivation is in principle similar to many textbooks, but the only (!) approximation up to now is linearisation. There are two special cases, which can be handled easily separately: 1.) =0, kx=0, ky0 2.) =0, kx=0, ky=0
Boundary conditions: below: linearisation of the free-slip condition: above: assume kz=const. case kz² < 0 : exponentially growing solution seems unphysical choose and omit term ~B case kz² > 0 : no energy transport downwards (Smith, 1980) choose and omit term ~B
Choice of the base state: N = const. kz ~ const. (= 2nd approx.!) experience: kz=const. is fulfilled rather good max. height: ~ 35 km • Some features of 'Lin_Mountain': • written in C/C++ • Makefile for Linux available • Binary output (with an additional .ctl-File for GrADS) • Use of FFT (Numerical recipes, Press et al.) for • estimation about the accuracy of the approximation 'kz=const.'
512 nx 2 ny 160 nz 200.0 dx 200.0 dy 100.0 dz 10.0 U_0 293.0 T_00 1.0e5 p_00 0.01 N_BV 1 delta1__0=hydrostat/1=nichthydrostat 1 delta2__0=inkompress/1=kompress 1 delta3__0=flach/1=tief 1 delta4__0=kleine/1=grosse_Machzahl topo.d oro_datei topo_ft.d oro_ft_datei Gaussdamm oro_typ 10000.0 oro_breite 10.0 oro_hoehe gauss_2D Ausgabekennung ./ Ausgabepfad 100 ix_min 400 ix_max 1 ix_step 0 iy_min 0 iy_max 1 iy_step 0 iz_min 159 iz_max 1 FFT Input-File: (abc.inp)
Example 1: 2D-test case gaussian hill (for the Input-Data: slide before) w [m/s]
w [m/s] Example 3: 3D Gaussian Hill
Optimization of horizontal advection: up to COSMO 4.3: 'advection operators' = a subroutine acting on every single grid point compiler has problems to optimize loops since COSMO 4.4: advection routines using 'field operations' (and additionally the DFI modifications by Lucio Torrisi) • Efficiency gain for routine COSMO-DE at DWD (IBM): • speedup of the horizontal advection alone: ~ 3 times faster • overall reduction of model run time:~ 1 Min. / 20 Min. ~ 5% Furthermore, some inconsistencies using metrical factors could be repaired acrlat(j,1)acrlat(j,2)lent to an error of ~ -0.05% in the term v dw/dy
ngranularity: Each sample hit covers 4 bytes. Time: 16425.58 seconds called/total parents index %time self descendents called+self name index called/total children 53.46 12371.49 1/1 .__start [2] [1] 75.6 53.46 12371.49 1 .lmorg [1] 0.01 9450.25 8514/8514 .organize_dynamics [3] 0.73 2529.11 8514/8514 .organize_physics [6] 72.20 87.39 4256/4256 .initialize_loop [32] 0.01 119.25 8515/8515 .organize_diagnosis [39] 0.01 61.03 4258/8514 .organize_data [38] 0.04 36.55 12770/12770 .organize_assimilation [63] 10.19 0.82 4256/4256 .near_surface [79] 0.02 3.36 4256/4256 .exchange_runge_kutta [98] ... ----------------------------------------------- ... ----------------------------------------------- 0.01 9450.25 8514/8514 .lmorg [1] [3] 57.5 0.01 9450.25 8514 .organize_dynamics [3] 262.68 9025.25 4256/4256 .__src_runge_kutta_NMOD_org_runge_kutta [4] 113.02 49.28 4256/4256 .__src_relaxation_NMOD_sardass [31] 0.02 0.00 1/1 .init_dynamics [302] 0.00 0.00 1/1 .__src_relaxation_NMOD_init_relaxation [641] 0.00 0.00 1/1 .input_dynctl [1105] ----------------------------------------------- 262.68 9025.25 4256/4256 .organize_dynamics [3] [4] 56.5 262.68 9025.25 4256 .__src_runge_kutta_NMOD_org_runge_kutta [4] 3557.22 113.61 12768/12768 .__fast_waves_rk_NMOD_fast_waves_runge_kutta [5] ===> 602.40 1593.07 12768/12768 .__src_advection_rk_NMOD_advection_alt [7] 31.91 985.85 4256/4256 .__src_advection_rk_NMOD_advection_pd [11] 930.86 0.00 12768/12768 .__src_slow_tendencies_rk_NMOD_complete_tendencies_uvwtpp [13] ----------------------------------------------- 602.40 1593.07 12768/12768 .__src_runge_kutta_NMOD_org_runge_kutta [4] [7] 13.4 602.40 1593.07 12768 .__src_advection_rk_NMOD_advection_alt [7] 1574.75 0.00 18446744061366576640/18446744062714877440 .__numeric_utilities_rk_NMOD_udsdx [8] ... Old version with 'advection operator' (gprof - output)
ngranularity: Each sample hit covers 4 bytes. Time: 15485.44 seconds called/total parents index %time self descendents called+self name index called/total children 6.6s <spontaneous> [1] 72.9 426.78 10867.61 .lmorg [1] 0.01 8139.09 8514/8514 .organize_dynamics [2] 0.56 2496.58 8514/8514 .organize_physics [5] 0.02 120.22 8515/8515 .organize_diagnosis [42] 0.02 60.68 4258/8514 .organize_data [41] 0.00 38.48 12770/12770 .organize_assimilation [65] 10.45 0.95 4256/4256 .near_surface [81] ... ----------------------------------------------- 0.01 8139.09 8514/8514 .lmorg [1] [2] 52.6 0.01 8139.09 8514 .organize_dynamics [2] 264.96 7711.97 4256/4256 .__src_runge_kutta_NMOD_org_runge_kutta [3] ... ----------------------------------------------- 264.96 7711.97 4256/4256 .organize_dynamics [2] [3] 51.5 264.96 7711.97 4256 .__src_runge_kutta_NMOD_org_runge_kutta [3] 3592.24 180.20 12768/12768 .__fast_waves_rk_NMOD_fast_waves_runge_kutta [4] 31.61 1064.15 4256/4256 .__src_advection_rk_NMOD_advection_pd [7] 1041.89 0.00 12768/12768 .__src_slow_tendencies_rk_NMOD_complete_tendencies_uvwtpp [8] ===> 286.16 307.05 12768/12768 .__src_advection_rk_NMOD_advection [14] 65.29 394.44 4256/4256 .__hori_diffusion_NMOD_comp_hori_diff [17] ... ----------------------------------------------- ... ----------------------------------------------- 286.16 307.05 12768/12768 .__src_runge_kutta_NMOD_org_runge_kutta [3] [14] 3.8 286.16 307.05 12768 .__src_advection_rk_NMOD_advection [14] 0.04 282.22 76608/76608 .__src_advection_rk_NMOD_horiz_adv_driver [24] 0.07 18.80 25536/206951 .__environment_NMOD_exchg_boundaries [35] 5.76 0.00 56179200/1404480000 .__numeric_utilities_rk_NMOD_udsdx [38] New version with 'field operators'(gprof - output)