diff --git a/CONTRIBUTORS.md b/CONTRIBUTORS.md index 10e9ced2..c18d22d2 100644 --- a/CONTRIBUTORS.md +++ b/CONTRIBUTORS.md @@ -12,4 +12,5 @@ | DrTVockerodtMO | Terence Vockerodt | Met Office | 2026-01-08 | | MetBenjaminWent | Benjamin Went | Met Office | 2026-01-15 | | timgraham-Met | Tim Graham | Met Office | 2026-01-15 | -| mo-alistairp | Alistair Pirrie | Met Office | 2026-01-19 | \ No newline at end of file +| mo-alistairp | Alistair Pirrie | Met Office | 2026-01-19 | +| cjohnson-pi | Christine Johnson | Met Office | 2026-01-21 | \ No newline at end of file diff --git a/applications/linear_model/plot_convergence/plot_convergence.sh b/applications/linear_model/plot_convergence/plot_convergence.sh deleted file mode 100755 index ab6e366a..00000000 --- a/applications/linear_model/plot_convergence/plot_convergence.sh +++ /dev/null @@ -1,83 +0,0 @@ -############################################################################## -# (c) Crown copyright 2022 Met Office. All rights reserved. -# The file LICENCE, distributed with this code, contains details of the terms -# under which the code may be used. -############################################################################## - -#---------------------------------------------------------------------- -# Plots the convergence rate of the tangent linear model. -#---------------------------------------------------------------------- - -# INSTRUCTIONS TO RUN -# 1. Specify CONFIG -# 2. Run using . plot_convergence.sh, from the plot_convergence directory - -# SCIENCE DETAILS -# The relative linearisation error is -# E = || N(x+ gamma x') - N(x) - L(x) gamma x' || / || L(x) gamma x' || -# where N=nonlinear model, L=linear model, x=linearisation state -# x'=perturbation, gamma=scalar. -# From the Taylor series expansion, E(gamma) = O(gamma) i.e. of the order gamma -# So the relative error should be a linear function of gamma - -# SCRIPT STEPS -# 1. Produce the data: The integration test tl_test_timesteps is extended by -# running over 10 values of gamma, rather than 2 values of gamma. -# 2. Plot the data: The data is plotted for each prognostic variable. - -# EXTENSION -# The plot_configuration.nml can also be extended to other configurations e.g -# * increase the number of timesteps (timesteps_end) -# * increase the number of timesteps between updating the linearisation state -# (update_ls_frequency) - -#-------------------------------------------------------------------------- - -# CONFIG can be specified as either runge_kutta or semi_implicit -CONFIG=semi_implicit - -# Define directories using the current working directory -Working_dir=$PWD -Linear_dir="$(dirname "$PWD")" - -# Integration tests executable name -exe=$Linear_dir/test/$CONFIG - -# Build the integration tests, unless that has already been completed -if [ -f $exe ] ; then - echo "Do not need to build the executable as $exe exists" -else - echo "$exe does not exist, so now building the executable" - cd $Linear_dir - make integration-tests - - if [$? -ne 0 ]; then - echo "Error building the executable" - return - fi -fi - -# Setup the configuration - to test with 10 values of gamma -cd $Linear_dir/test/test_files/$CONFIG -cp ${CONFIG}_configuration.nml plot_configuration.nml -sed -i 's/number_gamma_values=2/number_gamma_values=10/g' plot_configuration.nml -if [ $? -ne 0 ]; then - echo "Error in creating plot_configuration.nml" - return -fi - -# Run the tl_test_timesteps integration test -echo "Running the integration test" -../../$CONFIG plot_configuration.nml test_timesteps > outfile -if [ $? -ne 0 ]; then - echo "Error in creating outfile data" - return -else - echo "Data created successfully" -fi - -# Plot the data, together with the expected gradient -echo "Plotting the data" -python $Working_dir/plot_convergence.py - - diff --git a/science/gungho/source/driver/gungho_step_mod.x90 b/science/gungho/source/driver/gungho_step_mod.x90 index cd516978..82b47579 100644 --- a/science/gungho/source/driver/gungho_step_mod.x90 +++ b/science/gungho/source/driver/gungho_step_mod.x90 @@ -122,11 +122,13 @@ module gungho_step_mod '(A,I0)' ) 'Start of timestep ', model_clock%get_step() call log_event( log_scratch_space, LOG_LEVEL_INFO ) - temp_corr_io_value => get_io_value( modeldb%values, 'temperature_correction_io_value') - call modeldb%values%get_value( 'total_dry_mass', total_dry_mass ) - call modeldb%values%get_value( 'total_energy', total_energy ) - call modeldb%values%get_value( 'total_energy_previous', & - total_energy_previous ) + if ( encorr_usage /= encorr_usage_none .or. write_conservation_diag ) then + temp_corr_io_value => get_io_value( modeldb%values, 'temperature_correction_io_value') + call modeldb%values%get_value( 'total_dry_mass', total_dry_mass ) + call modeldb%values%get_value( 'total_energy', total_energy ) + call modeldb%values%get_value( 'total_energy_previous', & + total_energy_previous ) + end if use_moisture = ( moisture_formulation /= moisture_formulation_dry ) diff --git a/science/linear/integration-test/nwp_gal9/ReadMe b/science/linear/integration-test/nwp_gal9/ReadMe new file mode 100644 index 00000000..911e9288 --- /dev/null +++ b/science/linear/integration-test/nwp_gal9/ReadMe @@ -0,0 +1 @@ +The nwp_gal9.nml matches the configuration from the default rose-app.conf, in linear_model. diff --git a/science/linear/integration-test/nwp_gal9/iodef.xml b/science/linear/integration-test/nwp_gal9/iodef.xml new file mode 100644 index 00000000..35ae0ee0 --- /dev/null +++ b/science/linear/integration-test/nwp_gal9/iodef.xml @@ -0,0 +1,124 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + performance + 1.0 + + + + true + 50 + true + + + + + diff --git a/science/linear/integration-test/nwp_gal9/nwp_gal9.f90 b/science/linear/integration-test/nwp_gal9/nwp_gal9.f90 new file mode 100644 index 00000000..2280e549 --- /dev/null +++ b/science/linear/integration-test/nwp_gal9/nwp_gal9.f90 @@ -0,0 +1,144 @@ +!----------------------------------------------------------------------------- +! (C) Crown copyright 2026 Met Office. All rights reserved. +! The file LICENCE, distributed with this code, contains details of the terms +! under which the code may be used. +!----------------------------------------------------------------------------- + +!>@brief The top level program for the tangent linear tests. +!>@details The default is to run all available tests - which +!! test whether the linear code is tangent linear to the +!! corresponding nonlinear code. +program nwp_gal9 + + use driver_collections_mod, only: init_collections, final_collections + use driver_time_mod, only: init_time, final_time + use driver_comm_mod, only: init_comm, final_comm + use driver_log_mod, only: init_logger, final_logger + use driver_config_mod, only: init_config, final_config + use driver_modeldb_mod, only: modeldb_type + use lfric_mpi_mod, only: global_mpi + use gungho_mod, only: gungho_required_namelists + use log_mod, only: log_event, & + LOG_LEVEL_ERROR, & + LOG_LEVEL_INFO + use linear_driver_mod, only: initialise, finalise + use tl_test_driver_mod, only: run_timesteps, & + run_transport_control, & + run_semi_imp_alg, & + run_rhs_alg + + implicit none + + ! Model run working data set + type(modeldb_type) :: modeldb + character(*), parameter :: application_name = 'nwp_gal9' + character(:), allocatable :: filename + + ! Variables used for parsing command line arguments + integer :: length, status, nargs + character(len=0) :: dummy + character(len=:), allocatable :: program_name, test_flag + + ! Flags which determine the tests that will be carried out + logical :: do_test_timesteps = .false. + logical :: do_test_transport_control = .false. + logical :: do_test_semi_imp_alg = .false. + logical :: do_test_rhs_alg = .false. + + ! Usage message to print + character(len=256) :: usage_message + + modeldb%mpi => global_mpi + + call modeldb%configuration%initialise( application_name, table_len=10 ) + call modeldb%values%initialise('values', 5) + + call log_event( 'TL testing running ...', LOG_LEVEL_INFO ) + + ! Create the depository, prognostics and diagnostics field collections + call modeldb%fields%add_empty_field_collection("depository", table_len = 100) + call modeldb%fields%add_empty_field_collection("prognostic_fields", & + table_len = 100) + call modeldb%fields%add_empty_field_collection("diagnostic_fields", & + table_len = 100) + call modeldb%fields%add_empty_field_collection("lbc_fields", & + table_len = 100) + call modeldb%fields%add_empty_field_collection("radiation_fields", & + table_len = 100) + call modeldb%fields%add_empty_field_collection("fd_fields", & + table_len = 100) + + + call modeldb%io_contexts%initialise(application_name, 100) + + ! Parse command line parameters + call get_command_argument( 0, dummy, length, status ) + allocate(character(length)::program_name) + call get_command_argument( 0, program_name, length, status ) + nargs = command_argument_count() + + ! Print out usage message if wrong number of arguments is specified + if (nargs /= 2) then + write(usage_message,*) "Usage: ",trim(program_name), & + " " // & + " test_XXX with XXX in { " // & + " timesteps, " // & + " transport_control, " // & + " semi_imp_alg, " // & + " rhs_alg, " // & + " } " + call log_event( trim(usage_message), LOG_LEVEL_ERROR ) + end if + + call get_command_argument( 1, dummy, length, status ) + allocate( character(length) :: filename ) + call get_command_argument( 1, filename, length, status ) + + call get_command_argument( 2, dummy, length, status ) + allocate(character(length)::test_flag) + call get_command_argument( 2, test_flag, length, status ) + + ! Choose test case depending on flag provided in the first command + ! line argument + select case (trim(test_flag)) + case ("test_timesteps") + do_test_timesteps = .true. + case ("test_transport_control") + do_test_transport_control = .true. + case ("test_semi_imp_alg") + do_test_semi_imp_alg = .true. + case ("test_rhs_alg") + do_test_rhs_alg = .true. + case default + call log_event( "Unknown test", LOG_LEVEL_ERROR ) + end select + + call init_comm( application_name, modeldb ) + call init_config( filename, gungho_required_namelists, & + modeldb%configuration ) + call init_logger( modeldb%mpi%get_comm(), application_name ) + call init_collections() + call init_time( modeldb ) + call initialise( application_name, modeldb ) + + if (do_test_timesteps) then + call run_timesteps(modeldb) + endif + if (do_test_transport_control) then + call run_transport_control(modeldb) + endif + if (do_test_rhs_alg) then + call run_rhs_alg(modeldb) + endif + if (do_test_semi_imp_alg) then + call run_semi_imp_alg(modeldb) + endif + + call finalise( application_name, modeldb ) + call final_time( modeldb ) + call final_collections() + call final_logger( application_name ) + call final_config() + call final_comm( modeldb ) + +end program nwp_gal9 diff --git a/science/linear/integration-test/nwp_gal9/nwp_gal9.py b/science/linear/integration-test/nwp_gal9/nwp_gal9.py new file mode 100755 index 00000000..aa141417 --- /dev/null +++ b/science/linear/integration-test/nwp_gal9/nwp_gal9.py @@ -0,0 +1,105 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +############################################################################## +# (c) Crown copyright 2026 Met Office. All rights reserved. +# The file LICENCE, distributed with this code, contains details of the terms +# under which the code may be used. +############################################################################## +''' +Run the linear model integration tests for the default (nwp_gal9) configuration + +''' +import os +import re +import sys + + +from testframework import LFRicLoggingTest, TestEngine, TestFailed + + +class TLTest(LFRicLoggingTest): + ''' + Run the linear model integration tests + ''' + + def __init__(self, flag): + self._flag = flag + if 'MPIEXEC_BROKEN' in os.environ: + TLTest.set_mpiexec_broken() + super(TLTest, self).__init__([sys.argv[1], + 'resources/nwp_gal9_configuration.nml', + 'test_' + self._flag], + processes=1, + name='tl_test.Log') + + def test(self, return_code, out, err): + ''' + Error messages if the test failed to run + ''' + if return_code != 0: + message = 'Test program failed with exit code: {code}' + raise TestFailed(message.format(code=return_code), + stdout=out, stderr=err, + log=self.getLFRicLoggingLog()) + + # "out" becomes self.getLFRicLoggingLog() when PE>1 + if not self.test_passed(out): + message = 'Test {} failed' + raise TestFailed(message.format(self._flag), + stdout=out, stderr=err, + log=self.getLFRicLoggingLog()) + + return 'TL test : '+self._flag + + def test_passed(self, out): + ''' + Examine the output to see if the validity test passed + ''' + success = False + pattern = re.compile(r'\s+test\s+.*?:\s*PASS\s*$') + for line in out.split("\n"): + match = pattern.search(line) + if match: + success = True + return success + +class tl_test_semi_imp_alg(TLTest): + ''' + Test the semi implicit timestep + ''' + def __init__(self): + flag = "semi_imp_alg" + super(tl_test_semi_imp_alg, self).__init__(flag) + + +class tl_test_rhs_alg(TLTest): + ''' + Test the right hand side forcing for the mixed solver + ''' + def __init__(self): + flag = "rhs_alg" + super(tl_test_rhs_alg, self).__init__(flag) + +class tl_test_transport_control(TLTest): + ''' + Test the transport + ''' + def __init__(self): + flag = "transport_control" + super(tl_test_transport_control, self).__init__(flag) + + +class tl_test_timesteps(TLTest): + ''' + Test running over multiple timesteps + ''' + def __init__(self): + flag = "timesteps" + super(tl_test_timesteps, self).__init__(flag) + +if __name__ == '__main__': + TestEngine.run(tl_test_rhs_alg()) + TestEngine.run(tl_test_transport_control()) + TestEngine.run(tl_test_semi_imp_alg()) + TestEngine.run(tl_test_timesteps()) + diff --git a/science/linear/integration-test/nwp_gal9/resources/mesh_C12_MG.nc b/science/linear/integration-test/nwp_gal9/resources/mesh_C12_MG.nc new file mode 100644 index 00000000..ebe51890 Binary files /dev/null and b/science/linear/integration-test/nwp_gal9/resources/mesh_C12_MG.nc differ diff --git a/science/linear/integration-test/nwp_gal9/resources/nwp_gal9_configuration.nml b/science/linear/integration-test/nwp_gal9/resources/nwp_gal9_configuration.nml new file mode 100644 index 00000000..c687adfe --- /dev/null +++ b/science/linear/integration-test/nwp_gal9/resources/nwp_gal9_configuration.nml @@ -0,0 +1,412 @@ +&base_mesh +file_prefix='resources/mesh_C12_MG', +geometry='spherical', +prepartitioned=.false., +prime_mesh_name='dynamics', +topology='fully_periodic', +/ +&boundaries +limited_area=.false., +transport_overwrite_freq='final', +/ +&checks +limit_cfl=.false., +/ +§ion_choice +aerosol='none', +boundary_layer='none', +chemistry='none', +cloud='none', +dynamics='gungho', +external_forcing=.false., +iau=.false., +iau_sst=.false., +iau_surf=.false., +methane_oxidation=.false., +orographic_drag='none', +radiation='none', +spectral_gwd='none', +stochastic_physics='none', +surface='none', +/ +&convection +dx_ref=50000.0, +l_cvdiag_ctop_qmax=.false., +qlmin=4.0e-4, +resdep_precipramp=.false., +/ +&cosp +l_cosp=.false., +/ +&damping_layer +dl_base=40000.0, +dl_str=0.05, +dl_type='standard', +/ +&departure_points +horizontal_limit='cap', +horizontal_method='ffsl', +n_dep_pt_iterations=1, +share_stencil_extent=.true., +vertical_limit='exponential', +vertical_method='timeaverage', +vertical_sorting=.false., +/ +&energy_correction +encorr_usage='none', +integral_method='fd', +/ +&extrusion +domain_height=80000.0, +method='um_L70_50t_20s_80km', +number_of_layers=70, +planet_radius=6371229.0, +stretching_height=17507.0, +stretching_method='smooth', +/ +&files +ancil_directory='/data/users/lfricadmin/data/ancils/basic-gal/yak/C12', +checkpoint_stem_name='', +diag_stem_name='diagGungho', +ls_directory='/data/users/lfricadmin/data/tangent-linear/Ticket354', +ls_filename='final_ls', +orography_mean_ancil_path='orography/gmted_ramp2/qrparm.orog', +start_dump_directory='/data/users/lfricadmin/data/tangent-linear/Ticket354', +start_dump_filename='final_pert', +/ +&finite_element +cellshape='quadrilateral', +coord_order=1, +coord_system='native', +element_order_h=0, +element_order_v=0, +rehabilitate=.true., +vorticity_in_w1=.false., +/ +&formulation +dlayer_on=.false., +dry_static_adjust=.true., +eos_method='sampled', +exner_from_eos=.false., +horizontal_physics_predictor=.false., +horizontal_transport_predictor=.false., +init_exner_bt=.true., +l_multigrid=.true., +lagged_orog=.true., +moisture_formulation='traditional', +moisture_in_solver=.true., +p2theta_vert=.true., +rotating=.true., +shallow=.true., +si_momentum_equation=.false., +theta_moist_source=.false., +use_multires_coupling=.false., +use_physics=.true., +use_wavedynamics=.true., +vector_invariant=.false., +/ +&helmholtz_solver +gcrk=8, +method='prec_only', +monitor_convergence=.false., +normalise=.true., +preconditioner='multigrid', +si_pressure_a_tol=1.0e-8, +si_pressure_maximum_iterations=400, +si_pressure_tolerance=1.0e-4, +/ +&iau_addinf_io +/ +&iau_addinf_io +/ +&iau_ainc_io +/ +&iau_ainc_io +/ +&iau_bcorr_io +/ +&iau +/ +&idealised +f_lon_deg=0.0, +perturb_init=.false., +test='none', +/ +&ideal_surface +canopy_height=19.01,16.38,0.79,1.26,1.0, +leaf_area_index=5.0,4.0,1.5,1.5,1.5, +n_snow_layers=11*0, +snow_depth=11*0.0, +snow_layer_ice_mass=27*0.0, +snow_layer_temp=27*273.0, +snow_layer_thickness=27*0.0, +soil_moisture=15.86,98.861,274.35,862.27, +soil_temperature=284.508,286.537,289.512,293.066, +surf_tile_fracs=9*0.0,1.0,0.0, +surf_tile_temps=9*295.0,300.0,265.0, +tile_snow_mass=11*0.0, +/ +&initialization +ancil_option='none', +coarse_aerosol_ancil=.false., +coarse_orography_ancil=.false., +coarse_ozone_ancil=.false., +init_option='fd_start_dump', +lbc_option='none', +ls_option='file', +model_eos_height=100, +n_orog_smooth=0, +read_w2h_wind=.true., +sea_ice_source='ancillary', +snow_source='start_dump', +w0_orography_mapping=.false., +zero_w2v_wind=.false., +/ +&initial_density +density_background=0.1, +density_max=2.0, +r1=0.4, +r2=0.4, +x1=0.4, +x2=-0.4, +y1=0.0, +y2=0.0, +z1=0.0, +z2=0.0, +/ +&initial_pressure +method='balanced', +surface_pressure=1000.0e2, +/ +&initial_temperature +bvf_square=0.0001, +pert_centre=60.0, +pert_width_scaling=1.0, +perturb='none', +theta_surf=300.0, +/ +&initial_vapour +/ +&initial_wind +nl_constant=0.0, +profile='constant_uv', +sbr_angle_lat=0.0, +sbr_angle_lon=0.0, +smp_init_wind=.true., +u0=2.0, +v0=0.0, +wind_time_period=0.0, +/ +&io +checkpoint_read=.false., +checkpoint_write=.false., +counter_output_suffix='counter.txt', +diag_active_files='lfric_diag', +diag_always_on_sampling=.false., +diagnostic_frequency=8, +file_convention='UGRID', +multifile_io=.false., +nodal_output_on_w3=.false., +subroutine_counters=.false., +subroutine_timers=.true., +timer_output_path='timer.txt', +use_xios_io=.true., +write_conservation_diag=.false., +write_diag=.false., +write_dump=.false., +write_fluxes=.false., +write_minmax_tseries=.false., +/ +&linear +fixed_ls=.false., +l_stabilise_bl=.false., +ls_read_w2h=.false., +max_bl_stabilisation=0.75, +n_bl_levels_to_stabilise=15, +pert_option='file', +/ +&logging +log_to_rank_zero_only=.false., +run_log_level='debug', +/ +&mixed_solver +eliminate_variables='discrete', +fail_on_non_converged=.true., +gcrk=4, +guess_np1=.false., +mixed_solver_a_tol=1.0e-3, +monitor_convergence=.true., +normalise=.true., +reference_reset_time=1800.0, +si_maximum_iterations=10, +si_method='block_gcr', +si_preconditioner='pressure', +si_tolerance=1.0e-1, +split_w=.true., +/ +&mixing +leonard_term=.false., +smagorinsky=.false., +viscosity=.false., +viscosity_mu=0.0, +/ +&multigrid +chain_mesh_tags='dynamics','multigrid_l1','multigrid_l2', +multigrid_chain_nitems=3, +n_coarsesmooth=4, +n_postsmooth=2, +n_presmooth=2, +smooth_relaxation=0.8, +/ +&esm_couple +l_esm_couple_test=.false., +/ +&orography +orog_init_option='ancil', +/ +&partitioning +generate_inner_halos=.false., +panel_decomposition='auto', +panel_xproc=6, +panel_yproc=1, +partitioner='cubedsphere', +/ +&physics +configure_segments=.false., +limit_drag_incs=.false., +sample_physics_scalars=.true., +sample_physics_winds=.true., +sample_physics_winds_correction=.false., +/ +&planet +cp=1005.0, +gravity=9.80665, +omega=7.292116E-5, +p_zero=100000.0, +rd=287.05, +scaling_factor=1.0, +/ +&radiative_gases +cfc113_rad_opt='off', +cfc11_mix_ratio=1.110e-09, +cfc11_rad_opt='constant', +cfc12_mix_ratio=2.187e-09, +cfc12_rad_opt='constant', +ch4_mix_ratio=1.006e-06, +ch4_rad_opt='constant', +co2_mix_ratio=6.002e-04, +co2_rad_opt='constant', +co_rad_opt='off', +cs_rad_opt='off', +h2_rad_opt='off', +h2o_rad_opt='prognostic', +hcfc22_rad_opt='off', +hcn_rad_opt='off', +he_rad_opt='off', +hfc134a_rad_opt='off', +k_rad_opt='off', +l_cts_fcg_rates=.false., +li_rad_opt='off', +n2_rad_opt='off', +n2o_mix_ratio=4.945e-07, +n2o_rad_opt='constant', +na_rad_opt='off', +nh3_rad_opt='off', +o2_mix_ratio=0.2314, +o2_rad_opt='constant', +o3_rad_opt='ancil', +rb_rad_opt='off', +so2_rad_opt='off', +tio_rad_opt='off', +vo_rad_opt='off', +/ +&solver +gcrk=18, +maximum_iterations=7, +method='chebyshev', +monitor_convergence=.false., +preconditioner='diagonal', +tolerance=1.0e-6, +/ +&specified_surface +/ +&time +calendar='timestep', +calendar_origin='2016-01-01 15:00:00', +calendar_start='2016-01-01 15:00:00', +calendar_type='gregorian', +timestep_end='3', +timestep_start='1', +/ +×tepping +alpha=0.55, +dt=1800, +inner_iterations=2, +method='semi_implicit', +outer_iterations=2, +runge_kutta_method='forward_euler', +spinup_alpha=.false.tau, +tau_r=1.0, +tau_t=1.0, +tau_u=0.55, +/ +&transport +adjust_theta=.false., +adjust_vhv_wind=.false., +ageofair_reset_level=10, +broken_w2_projection=.false., +calculate_detj='upwind', +cap_density_predictor=0.5, +cfl_mol_1d_stab=1.0, +cfl_mol_2d_stab=1.0, +cfl_mol_3d_stab=1.0, +cheap_update=.false., +consistent_metric=.false., +dep_pt_stencil_extent=3, +dry_field_name='density', +enforce_min_value=5*.false., +equation_form=1,2,2,2,2, +ffsl_inner_order=0, +ffsl_outer_order=1, +ffsl_splitting=5*1, +ffsl_unity_3d=.false., +ffsl_vertical_order=2,2,1,2,2, +field_names='density','potential_temperature','wind','moisture', +'con_tracer', +fv_horizontal_order=2, +fv_vertical_order=2, +horizontal_method=5*1, +horizontal_monotone=5*1, +log_space=5*.false., +max_vert_cfl_calc='dep_point', +min_val_abs_tol=-1.0e-12, +min_val_max_iterations=10, +min_val_method='iterative', +min_value=0.0,0.0,-99999999.0,0.0,0.0, +oned_reconstruction=.false., +operators='fv', +panel_edge_high_order=.true., +panel_edge_treatment='none', +profile_size=5, +reversible=5*.false., +runge_kutta_method='ssp3', +scheme=5*1, +si_outer_transport='none', +slice_order='parabola', +special_edges_monotone=5*1, +splitting=5*1, +substep_transport='off', +theta_dispersion_correction=.false., +theta_variable='dry', +transport_ageofair=.false., +use_density_predictor=.false., +vertical_method=5*1, +vertical_monotone=5*1, +vertical_monotone_order=5*3, +vertical_sl_order='cubic', +wind_mono_top=.false., +/ +&validity_test +number_gamma_values=2, +update_ls_frequency=1, +/ diff --git a/science/linear/integration-test/runge_kutta/resources/runge_kutta_configuration.nml b/science/linear/integration-test/runge_kutta/resources/runge_kutta_configuration.nml index 9298f96f..c396ffe9 100644 --- a/science/linear/integration-test/runge_kutta/resources/runge_kutta_configuration.nml +++ b/science/linear/integration-test/runge_kutta/resources/runge_kutta_configuration.nml @@ -131,7 +131,7 @@ subroutine_counters=.false., subroutine_timers=.true., use_xios_io=.false., write_conservation_diag=.false., -write_diag=.true., +write_diag=.false., write_dump=.false., write_fluxes=.false., write_minmax_tseries=.false., @@ -194,7 +194,7 @@ tolerance=1.0e-6, / &time calendar='timestep', -timestep_end='1', +timestep_end='6', timestep_start='1', calendar_type='gregorian', calendar_start='2016-01-01 15:00:00', diff --git a/science/linear/integration-test/runge_kutta/runge_kutta.f90 b/science/linear/integration-test/runge_kutta/runge_kutta.f90 index f460201a..301f0b7f 100644 --- a/science/linear/integration-test/runge_kutta/runge_kutta.f90 +++ b/science/linear/integration-test/runge_kutta/runge_kutta.f90 @@ -10,21 +10,19 @@ !! corresponding nonlinear code. program runge_kutta - use configuration_mod, only: read_configuration, final_configuration use driver_collections_mod, only: init_collections, final_collections use driver_time_mod, only: init_time, final_time + use driver_comm_mod, only: init_comm, final_comm + use driver_log_mod, only: init_logger, final_logger + use driver_config_mod, only: init_config, final_config use driver_modeldb_mod, only: modeldb_type - use halo_comms_mod, only: initialise_halo_comms, & - finalise_halo_comms - use lfric_mpi_mod, only: create_comm, destroy_comm, global_mpi, & - lfric_comm_type - use log_mod, only: initialise_logging, finalise_logging, & - log_event, & + use lfric_mpi_mod, only: global_mpi + use gungho_mod, only: gungho_required_namelists + use log_mod, only: log_event, & LOG_LEVEL_ERROR, & LOG_LEVEL_INFO - use tl_test_driver_mod, only: initialise, & - finalise, & - run_timesteps, & + use linear_driver_mod, only: initialise, finalise + use tl_test_driver_mod, only: run_timesteps_random, & run_kinetic_energy_gradient, & run_advect_density_field, & run_advect_theta_field, & @@ -39,7 +37,6 @@ program runge_kutta ! Model run working data set type(modeldb_type) :: modeldb character(*), parameter :: application_name = "runge_kutta" - character(:), allocatable :: filename ! Variables used for parsing command line arguments @@ -47,8 +44,6 @@ program runge_kutta character(len=0) :: dummy character(len=:), allocatable :: program_name, test_flag - type(lfric_comm_type) :: communicator - ! Flags which determine the tests that will be carried out logical :: do_test_timesteps = .false. logical :: do_test_kinetic_energy_gradient = .false. @@ -65,11 +60,8 @@ program runge_kutta modeldb%mpi => global_mpi - call create_comm( communicator ) - call modeldb%mpi%initialise( communicator ) - call initialise_logging( communicator%get_comm_mpi_val(), & - "linear_integration-runge_kutta-test" ) - call initialise_halo_comms( communicator ) + call modeldb%configuration%initialise( application_name, table_len=10 ) + call modeldb%values%initialise('values', 5) call log_event( 'TL testing running ...', LOG_LEVEL_INFO ) @@ -86,7 +78,7 @@ program runge_kutta call modeldb%fields%add_empty_field_collection("fd_fields", & table_len = 100) - call modeldb%io_contexts%initialise(program_name, 100) + call modeldb%io_contexts%initialise(application_name, 100) ! Parse command line parameters call get_command_argument( 0, dummy, length, status ) @@ -145,16 +137,16 @@ program runge_kutta call log_event( "Unknown test", LOG_LEVEL_ERROR ) end select - call modeldb%configuration%initialise( program_name, table_len=10 ) - call read_configuration( filename, modeldb%configuration ) - deallocate( filename ) - + call init_comm( application_name, modeldb ) + call init_config( filename, gungho_required_namelists, & + modeldb%configuration ) + call init_logger( modeldb%mpi%get_comm(), application_name ) call init_collections() call init_time( modeldb ) - call initialise( application_name, modeldb, modeldb%calendar ) + call initialise( application_name, modeldb ) if (do_test_timesteps) then - call run_timesteps(modeldb) + call run_timesteps_random(modeldb) endif if (do_test_kinetic_energy_gradient) then call run_kinetic_energy_gradient(modeldb) @@ -184,9 +176,8 @@ program runge_kutta call finalise( application_name, modeldb ) call final_time( modeldb ) call final_collections() - call final_configuration() - call finalise_halo_comms() - call finalise_logging() - call destroy_comm() + call final_logger( application_name ) + call final_config() + call final_comm( modeldb ) end program runge_kutta diff --git a/science/linear/integration-test/semi_implicit/resources/semi_implicit_configuration.nml b/science/linear/integration-test/semi_implicit/resources/semi_implicit_configuration.nml index 2b208a45..6be19971 100644 --- a/science/linear/integration-test/semi_implicit/resources/semi_implicit_configuration.nml +++ b/science/linear/integration-test/semi_implicit/resources/semi_implicit_configuration.nml @@ -132,7 +132,7 @@ subroutine_counters=.false., subroutine_timers=.true., use_xios_io=.false., write_conservation_diag=.false., -write_diag=.true., +write_diag=.false., write_dump=.false., write_fluxes=.false., write_minmax_tseries=.false., diff --git a/science/linear/integration-test/semi_implicit/semi_implicit.f90 b/science/linear/integration-test/semi_implicit/semi_implicit.f90 index af625729..57ac5a80 100644 --- a/science/linear/integration-test/semi_implicit/semi_implicit.f90 +++ b/science/linear/integration-test/semi_implicit/semi_implicit.f90 @@ -10,35 +10,25 @@ !! corresponding nonlinear code. program semi_implicit - use configuration_mod, only: read_configuration, final_configuration use driver_collections_mod, only: init_collections, final_collections use driver_time_mod, only: init_time, final_time + use driver_comm_mod, only: init_comm, final_comm + use driver_log_mod, only: init_logger, final_logger + use driver_config_mod, only: init_config, final_config use driver_modeldb_mod, only: modeldb_type - use halo_comms_mod, only: initialise_halo_comms, finalise_halo_comms - use lfric_mpi_mod, only: global_mpi, & - create_comm, destroy_comm, & - lfric_comm_type - use log_mod, only: initialise_logging, finalise_logging, & - log_event, & + use lfric_mpi_mod, only: global_mpi + use gungho_mod, only: gungho_required_namelists + use log_mod, only: log_event, & LOG_LEVEL_ERROR, & LOG_LEVEL_INFO - use namelist_collection_mod, only: namelist_collection_type - use tl_test_driver_mod, only: initialise, & - finalise, & - run_timesteps, & - run_transport_control, & - run_semi_imp_alg, & - run_rhs_sample_eos, & - run_rhs_project_eos, & - run_rhs_alg + use linear_driver_mod, only: initialise, finalise + use tl_test_driver_mod, only: run_timesteps_random implicit none ! Model run working data set type(modeldb_type) :: modeldb - character(*), parameter :: application_name = 'semi_implicit' - character(:), allocatable :: filename ! Variables used for parsing command line arguments @@ -46,26 +36,16 @@ program semi_implicit character(len=0) :: dummy character(len=:), allocatable :: program_name, test_flag - type(lfric_comm_type) :: communicator - ! Flags which determine the tests that will be carried out logical :: do_test_timesteps = .false. - logical :: do_test_transport_control = .false. - logical :: do_test_semi_imp_alg = .false. - logical :: do_test_rhs_alg = .false. - logical :: do_test_rhs_project_eos = .false. - logical :: do_test_rhs_sample_eos = .false. ! Usage message to print character(len=256) :: usage_message modeldb%mpi => global_mpi - call create_comm( communicator ) - call modeldb%mpi%initialise( communicator ) - call initialise_logging( communicator%get_comm_mpi_val(), & - "linear_interface-semi_implicit-test" ) - call initialise_halo_comms( communicator ) + call modeldb%configuration%initialise( application_name, table_len=10 ) + call modeldb%values%initialise('values', 5) call log_event( 'TL testing running ...', LOG_LEVEL_INFO ) @@ -82,7 +62,8 @@ program semi_implicit call modeldb%fields%add_empty_field_collection("fd_fields", & table_len = 100) - call modeldb%io_contexts%initialise(program_name, 100) + + call modeldb%io_contexts%initialise(application_name, 100) ! Parse command line parameters call get_command_argument( 0, dummy, length, status ) @@ -99,8 +80,6 @@ program semi_implicit " transport_control, " // & " semi_imp_alg, " // & " rhs_alg, " // & - " rhs_project_eos, " // & - " rhs_sample_eos, " // & " } " call log_event( trim(usage_message), LOG_LEVEL_ERROR ) end if @@ -118,53 +97,27 @@ program semi_implicit select case (trim(test_flag)) case ("test_timesteps") do_test_timesteps = .true. - case ("test_transport_control") - do_test_transport_control = .true. - case ("test_semi_imp_alg") - do_test_semi_imp_alg = .true. - case ("test_rhs_alg") - do_test_rhs_alg = .true. - case ("test_rhs_project_eos") - do_test_rhs_project_eos = .true. - case ("test_rhs_sample_eos") - do_test_rhs_sample_eos = .true. case default call log_event( "Unknown test", LOG_LEVEL_ERROR ) end select - call modeldb%configuration%initialise( program_name, table_len=10 ) - call read_configuration( filename, modeldb%configuration ) - deallocate( filename ) - + call init_comm( application_name, modeldb ) + call init_config( filename, gungho_required_namelists, & + modeldb%configuration ) + call init_logger( modeldb%mpi%get_comm(), application_name ) call init_collections() call init_time( modeldb ) - call initialise( application_name, modeldb, modeldb%calendar ) + call initialise( application_name, modeldb ) if (do_test_timesteps) then - call run_timesteps(modeldb) - endif - if (do_test_transport_control) then - call run_transport_control(modeldb) - endif - if (do_test_rhs_alg) then - call run_rhs_alg(modeldb) - endif - if (do_test_rhs_project_eos) then - call run_rhs_project_eos(modeldb) - endif - if (do_test_rhs_sample_eos) then - call run_rhs_sample_eos(modeldb) - endif - if (do_test_semi_imp_alg) then - call run_semi_imp_alg(modeldb) + call run_timesteps_random(modeldb) endif call finalise( application_name, modeldb ) call final_time( modeldb ) call final_collections() - call final_configuration() - call finalise_halo_comms() - call finalise_logging() - call destroy_comm() + call final_logger( application_name ) + call final_config() + call final_comm( modeldb ) end program semi_implicit diff --git a/science/linear/integration-test/semi_implicit/semi_implicit.py b/science/linear/integration-test/semi_implicit/semi_implicit.py index af210aff..534eb989 100755 --- a/science/linear/integration-test/semi_implicit/semi_implicit.py +++ b/science/linear/integration-test/semi_implicit/semi_implicit.py @@ -63,43 +63,6 @@ def test_passed(self, out): success = True return success - -class tl_test_semi_imp_alg(TLTest): - ''' - Test the semi implicit timestep - ''' - def __init__(self): - flag = "semi_imp_alg" - super(tl_test_semi_imp_alg, self).__init__(flag) - - -class tl_test_rhs_alg(TLTest): - ''' - Test the right hand side forcing for the mixed solver - ''' - def __init__(self): - flag = "rhs_alg" - super(tl_test_rhs_alg, self).__init__(flag) - -class tl_test_rhs_sample_eos(TLTest): - def __init__(self): - flag = "rhs_sample_eos" - super(tl_test_rhs_sample_eos, self).__init__(flag) - -class tl_test_rhs_project_eos(TLTest): - def __init__(self): - flag = "rhs_project_eos" - super(tl_test_rhs_project_eos, self).__init__(flag) - -class tl_test_transport_control(TLTest): - ''' - Test the transport - ''' - def __init__(self): - flag = "transport_control" - super(tl_test_transport_control, self).__init__(flag) - - class tl_test_timesteps(TLTest): ''' Test running over multiple timesteps @@ -109,9 +72,4 @@ def __init__(self): super(tl_test_timesteps, self).__init__(flag) if __name__ == '__main__': - TestEngine.run( tl_test_rhs_sample_eos() ) - TestEngine.run( tl_test_rhs_project_eos() ) - TestEngine.run(tl_test_transport_control()) - TestEngine.run( tl_test_semi_imp_alg() ) - TestEngine.run( tl_test_rhs_alg() ) - TestEngine.run(tl_test_timesteps()) + TestEngine.run( tl_test_timesteps() ) diff --git a/science/linear/integration-test/tl_test/tl_test_advect_density_field_mod.x90 b/science/linear/integration-test/tl_test/tl_test_advect_density_field_mod.x90 index c7ee7136..b704f347 100644 --- a/science/linear/integration-test/tl_test/tl_test_advect_density_field_mod.x90 +++ b/science/linear/integration-test/tl_test/tl_test_advect_density_field_mod.x90 @@ -185,7 +185,6 @@ module tl_test_advect_density_field_mod call invoke( X_minus_Y( diff, n2_advection_inc, n1_advection_inc ), & inc_X_minus_Y( diff, p_advection_inc ), & X_innerproduct_X( norm_diff, diff ) ) - norm_diff = sqrt(norm_diff) write( log_scratch_space, '(A, E32.12, A, E32.12)' ) & 'gamma = ', gamma_u, ' norm = ', norm_diff diff --git a/science/linear/integration-test/tl_test/tl_test_advect_theta_field_mod.x90 b/science/linear/integration-test/tl_test/tl_test_advect_theta_field_mod.x90 index f712ba30..4630bdf1 100644 --- a/science/linear/integration-test/tl_test/tl_test_advect_theta_field_mod.x90 +++ b/science/linear/integration-test/tl_test/tl_test_advect_theta_field_mod.x90 @@ -188,8 +188,6 @@ module tl_test_advect_theta_field_mod inc_X_minus_Y( diff, p_advection_inc ), & X_innerproduct_X( norm_diff, diff ) ) - norm_diff = sqrt(norm_diff) - write( log_scratch_space, '(A, E32.12, A, E32.12)' ) & 'gamma = ' , gamma_u , ' norm = ' , norm_diff diff --git a/science/linear/integration-test/tl_test/tl_test_convergence_rate_check.f90 b/science/linear/integration-test/tl_test/tl_test_convergence_rate_check.f90 index a13655df..b090ba7f 100644 --- a/science/linear/integration-test/tl_test/tl_test_convergence_rate_check.f90 +++ b/science/linear/integration-test/tl_test/tl_test_convergence_rate_check.f90 @@ -6,7 +6,7 @@ !>@brief Test the convergence rate (Taylor remainder convergence). module tl_test_convergence_rate_check - use constants_mod, only: r_def, str_def + use constants_mod, only: r_def, str_def, i_def use log_mod, only: log_event, & log_scratch_space, & LOG_LEVEL_INFO @@ -14,17 +14,129 @@ module tl_test_convergence_rate_check implicit none private - public convergence_rate_check + public convergence_pass_string, & + array_convergence_rate_check, & + convergence_rate_check + +contains + + !> @brief Test the convergence rate, with application of square root. + !> @details If the convergence rate is not close to 4, within the + !! specified tolerance, set the pass string to FAIL. + !> @param[in] name Variable or test name + !> @param[in] norm Current norm + !> @param[in] norm_prev Previous norm + !> @param[in,out] pass_str Pass string (either PASS or FAIL) + !> @param[in] tol Tolerance + subroutine convergence_pass_string( name, norm, norm_prev, pass_str, tol ) + implicit none + + real(r_def), intent(in) :: norm, norm_prev + character(len=4), intent(inout) :: pass_str + real(r_def), intent(in) :: tol + character(str_def), intent(in) :: name + + real(r_def) :: conv_rate + + ! Let the error between the nonlinear (N) difference and the linear (L) be + ! E(g) = || N(x + g dx) - N(x) - g Ldx || = O(g^2) + ! where g is a scalar, x and dx are vectors, O is the order + ! and || . || = (x^T x)^1/2 is the L2 norm. + ! + ! Then the ratio + ! E(2g) / E(g) = O(4 g^2) / O(g^2) = 4 + ! i.e. we need to check whether the convergence rate is close to 4. + + ! The norms have not had a square root applied yet. + conv_rate = sqrt(norm_prev / norm) + + if ( abs( conv_rate - 4.0_r_def ) >= tol ) then + pass_str = "FAIL" + else + pass_str = "PASS" + end if + + write( log_scratch_space, '(A, A, A, E16.8, A, E16.8)') & + name, pass_str, " Convergence rate: ", conv_rate, " Tolerance: ", tol + call log_event( log_scratch_space, LOG_LEVEL_INFO ) + + end subroutine convergence_pass_string + + !> @brief Test the convergence rate for individual variables and the sum. + !> @details Calculate the convergence rate based on the norms + !> at two different iterations, and compare with the + !> expected value. Print out either PASS or FAIL, which + !> will then be used by the top level integration test. + !> @param[in] array_norm Norm at second iteration + !> @param[in] array_norm_prev Norm at first iteration + !> @param[in] array_names Name of the variable being tested + !> @param[in] n_variables Array lengths (number of variables) + !> @param[in] label Test name + !> @param[in] tol Tolerance value + subroutine array_convergence_rate_check( array_norm, array_norm_prev, array_names, n_variables, label, tol, indiv_tol) + integer(i_def), intent(in) :: n_variables + real(r_def), intent(in) :: array_norm(n_variables) + real(r_def), intent(in) :: array_norm_prev(n_variables) + character(str_def), intent(in) :: array_names(n_variables) + character(str_def), intent(in) :: label + real(r_def), optional, intent(in) :: tol + real(r_def), optional, intent(in) :: indiv_tol + real(r_def) :: tolerance, individual_tolerance + character(len=4) :: pass_str_arr(n_variables) + character(len=4) :: sum_pass_str, pass_str + character(str_def), parameter :: sum_name = "sum" + integer(i_def) :: i + real(r_def) :: sum, sum_prev - contains + call log_event( "Checking convergence rate", LOG_LEVEL_INFO ) + + if ( present(tol) ) then + tolerance = tol + else + tolerance = 1.0E-8_r_def + end if + + if ( present(indiv_tol) ) then + individual_tolerance = indiv_tol + else + individual_tolerance = 1.0E-8_r_def + end if + + sum = 0.0_r_def + sum_prev = 0.0_r_def + do i= 1, n_variables + ! Check individual convergence rates + call convergence_pass_string( & + array_names(i), array_norm(i), & + array_norm_prev(i), pass_str_arr(i), & + individual_tolerance ) + + ! Weighted sum + sum = sum + array_norm(i) / array_norm_prev(i) + sum_prev = sum_prev + array_norm_prev(i) / array_norm_prev(i) + end do + + call convergence_pass_string( & + sum_name, sum, sum_prev, sum_pass_str, tolerance) + + pass_str = "PASS" + do i= 1, n_variables + if ( pass_str_arr(i) == "FAIL" ) pass_str = "FAIL" + end do + if ( sum_pass_str == "FAIL" ) pass_str = "FAIL" + + write(log_scratch_space,'(" test",A32," : ",A4)') trim(label), pass_str + call log_event( log_scratch_space, LOG_LEVEL_INFO ) + + end subroutine array_convergence_rate_check !> @brief Calculate and test the convergence rate. !> @details Calculate the convergence rate based on the norms !> at two different iterations, and compare with the !> expected value. Print out either PASS or FAIL, which !> will then be used by the top level integration test. - !> @param[in] norm_diff Norm at first iteration - !> @param[in] norm_diff_prev Norm at second iteration + !> @param[in] norm_diff Norm at second iteration + !> @param[in] norm_diff_prev Norm at first iteration !> @param[in] label Name of the code being tested !> @param[in] tol Tolerance value subroutine convergence_rate_check( norm_diff, norm_diff_prev, label, tol ) @@ -35,7 +147,6 @@ subroutine convergence_rate_check( norm_diff, norm_diff_prev, label, tol ) character(str_def), intent(in) :: label real(r_def), optional, intent(in) :: tol real(r_def) :: tolerance - real(r_def) :: conv_rate character(len=4) :: pass_str if ( present(tol) ) then @@ -44,23 +155,8 @@ subroutine convergence_rate_check( norm_diff, norm_diff_prev, label, tol ) tolerance = 1.0E-8_r_def end if - conv_rate = norm_diff_prev/ norm_diff - - - - write( log_scratch_space, '(A)' ) & - "TL Test: " // trim(label) - call log_event( log_scratch_space, LOG_LEVEL_INFO ) - write( log_scratch_space, '(A, E16.8)') & - "Convergence rate: ", conv_rate - call log_event( log_scratch_space, LOG_LEVEL_INFO ) - - - if ( abs(conv_rate - 4.0_r_def ) < tolerance ) then - pass_str = "PASS" - else - pass_str = "FAIL" - end if + call convergence_pass_string( & + label, norm_diff, norm_diff_prev, pass_str, tolerance) write(log_scratch_space,'(" test",A32," : ",A4)') trim(label), pass_str call log_event( log_scratch_space, LOG_LEVEL_INFO ) diff --git a/science/linear/integration-test/tl_test/tl_test_driver_mod.f90 b/science/linear/integration-test/tl_test/tl_test_driver_mod.f90 index b93d850d..0d1c7e5e 100644 --- a/science/linear/integration-test/tl_test/tl_test_driver_mod.f90 +++ b/science/linear/integration-test/tl_test/tl_test_driver_mod.f90 @@ -48,13 +48,13 @@ module tl_test_driver_mod use tl_test_rhs_alg_mod, only : test_rhs_alg use tl_test_semi_imp_alg_mod, only : test_semi_imp_alg use tl_test_timesteps_alg_mod, only : test_timesteps + use tl_test_timesteps_random_alg_mod, only : test_timesteps_random implicit none private - public initialise, & - finalise, & - run_timesteps, & + public run_timesteps, & + run_timesteps_random, & run_kinetic_energy_gradient, & run_advect_density_field, & run_advect_theta_field, & @@ -72,109 +72,47 @@ module tl_test_driver_mod type(mesh_type), pointer :: mesh => null() type(mesh_type), pointer :: twod_mesh => null() - type(mesh_type), pointer :: aerosol_mesh => null() - type(mesh_type), pointer :: aerosol_twod_mesh => null() contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !>@brief Sets up the required state in preparation for run. - !>@param [in] program_name An identifier given to the model being run - !>@param [in,out] modeldb The structure that holds model state - !> - subroutine initialise( program_name, modeldb, calendar ) + !>@brief Tests the tangent linear model for multiple timesteps + !>@param [in,out] modeldb The structure that holds model state + subroutine run_timesteps(modeldb) implicit none - character(*), intent(in) :: program_name - type(modeldb_type), intent(inout) :: modeldb - class(calendar_type), intent(in) :: calendar - - type(gungho_time_axes_type) :: model_axes - type(io_value_type) :: temp_corr_io_value - type(io_value_type) :: random_seed_io_value - integer(i_def) :: random_seed_size - real(r_def), allocatable :: real_array(:) - - call modeldb%values%initialise( 'values', 5 ) - - ! Initialise infrastructure and setup constants - ! - call initialise_infrastructure( program_name, modeldb ) - - ! Add a place to store time axes in modeldb - call modeldb%values%add_key_value('model_axes', model_axes) + type(modeldb_type), intent(inout) :: modeldb - ! Get primary and 2D meshes for initialising model data mesh => mesh_collection%get_mesh(prime_mesh_name) twod_mesh => mesh_collection%get_mesh(mesh, TWOD) - ! Assume aerosol mesh is the same as dynamics mesh - aerosol_mesh => mesh - aerosol_twod_mesh => twod_mesh - - ! gungho_init_field() expects these values to exist. The dependency of - ! the linear application tests on this procedure will hopefully be resolved - ! in the future, at which point this initialisation may be removed. - ! - call temp_corr_io_value%init('temperature_correction_rate', [0.0_r_def]) - call modeldb%values%add_key_value( 'temperature_correction_io_value', & - temp_corr_io_value ) - call modeldb%values%add_key_value( 'total_dry_mass', 0.0_r_def ) - call modeldb%values%add_key_value( 'total_energy', 0.0_r_def ) - call modeldb%values%add_key_value( 'total_energy_previous', 0.0_r_def ) - if ( stochastic_physics == stochastic_physics_um ) then - ! Random seed for stochastic physics - call random_seed(size = random_seed_size) - allocate(real_array(random_seed_size)) - real_array(1:random_seed_size) = 0.0_r_def - call random_seed_io_value%init("random_seed", real_array) - call modeldb%values%add_key_value( 'random_seed_io_value', & - random_seed_io_value ) - deallocate(real_array) - end if - - ! Instantiate the fields stored in model_data - call create_model_data( modeldb, & - mesh, & - twod_mesh, & - aerosol_mesh, & - aerosol_twod_mesh ) - - ! Instantiate the linearisation state - call linear_create_ls( modeldb, mesh ) - - ! Initialise the fields stored in the model_data prognostics. This needs - ! to be done before initialise_model. - call initialise_model_data( modeldb, mesh, twod_mesh ) - - ! Model configuration initialisation - call initialise_model( mesh, & - modeldb ) - - ! Initialise the linearisation state - call linear_init_ls( mesh, twod_mesh, modeldb ) - - ! Finalise model - call finalise_model(modeldb) - - end subroutine initialise + call test_timesteps( modeldb, & + mesh, & + twod_mesh, & + modeldb%clock ) + + end subroutine run_timesteps !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !>@brief Tests the tangent linear model for multiple timesteps + !! using prescribed random data for the initial conditions !>@param [in,out] modeldb The structure that holds model state - subroutine run_timesteps(modeldb) + subroutine run_timesteps_random(modeldb) implicit none type(modeldb_type), intent(inout) :: modeldb - call test_timesteps( modeldb, & - mesh, & - twod_mesh, & - modeldb%clock ) + mesh => mesh_collection%get_mesh(prime_mesh_name) + twod_mesh => mesh_collection%get_mesh(mesh, TWOD) - end subroutine run_timesteps + call test_timesteps_random( modeldb, & + mesh, & + twod_mesh, & + modeldb%clock ) + + end subroutine run_timesteps_random !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !>@brief Tests the tangent linear model kinetic energy gradient kernel @@ -185,6 +123,9 @@ subroutine run_kinetic_energy_gradient(modeldb) type(modeldb_type), intent(inout) :: modeldb + mesh => mesh_collection%get_mesh(prime_mesh_name) + twod_mesh => mesh_collection%get_mesh(mesh, TWOD) + call test_kinetic_energy_gradient( modeldb, & mesh, & twod_mesh ) @@ -200,6 +141,9 @@ subroutine run_advect_density_field(modeldb) type(modeldb_type), intent(inout) :: modeldb + mesh => mesh_collection%get_mesh(prime_mesh_name) + twod_mesh => mesh_collection%get_mesh(mesh, TWOD) + call test_advect_density_field( modeldb, & mesh, & twod_mesh ) @@ -215,6 +159,9 @@ subroutine run_advect_theta_field(modeldb) type(modeldb_type), intent(inout) :: modeldb + mesh => mesh_collection%get_mesh(prime_mesh_name) + twod_mesh => mesh_collection%get_mesh(mesh, TWOD) + call test_advect_theta_field( modeldb, & mesh, & twod_mesh ) @@ -230,6 +177,9 @@ subroutine run_vorticity_advection(modeldb) type(modeldb_type), intent(inout) :: modeldb + mesh => mesh_collection%get_mesh(prime_mesh_name) + twod_mesh => mesh_collection%get_mesh(mesh, TWOD) + call test_vorticity_advection( modeldb, & mesh, & twod_mesh ) @@ -245,6 +195,9 @@ subroutine run_project_eos_pressure(modeldb) type(modeldb_type), intent(inout) :: modeldb + mesh => mesh_collection%get_mesh(prime_mesh_name) + twod_mesh => mesh_collection%get_mesh(mesh, TWOD) + call test_project_eos_pressure( modeldb, & mesh, & twod_mesh ) @@ -260,6 +213,9 @@ subroutine run_sample_eos_pressure(modeldb) type(modeldb_type), intent(inout) :: modeldb + mesh => mesh_collection%get_mesh(prime_mesh_name) + twod_mesh => mesh_collection%get_mesh(mesh, TWOD) + call test_sample_eos_pressure( modeldb, & mesh, & twod_mesh ) @@ -275,6 +231,9 @@ subroutine run_hydrostatic(modeldb) type(modeldb_type), intent(inout) :: modeldb + mesh => mesh_collection%get_mesh(prime_mesh_name) + twod_mesh => mesh_collection%get_mesh(mesh, TWOD) + call test_hydrostatic( modeldb, & mesh, & twod_mesh ) @@ -290,6 +249,9 @@ subroutine run_pressure_gradient_bd(modeldb) type(modeldb_type), intent(inout) :: modeldb + mesh => mesh_collection%get_mesh(prime_mesh_name) + twod_mesh => mesh_collection%get_mesh(mesh, TWOD) + call test_pressure_gradient_bd( modeldb, & mesh, & twod_mesh ) @@ -305,6 +267,8 @@ subroutine run_rk_alg(modeldb) type(modeldb_type), intent(inout) :: modeldb + mesh => mesh_collection%get_mesh(prime_mesh_name) + call test_rk_alg( modeldb, & mesh) @@ -319,6 +283,9 @@ subroutine run_transport_control(modeldb) type(modeldb_type), intent(inout) :: modeldb + mesh => mesh_collection%get_mesh(prime_mesh_name) + twod_mesh => mesh_collection%get_mesh(mesh, TWOD) + call test_transport_control( modeldb, & mesh, & twod_mesh ) @@ -334,6 +301,9 @@ subroutine run_semi_imp_alg(modeldb) type(modeldb_type), intent(inout) :: modeldb + mesh => mesh_collection%get_mesh(prime_mesh_name) + twod_mesh => mesh_collection%get_mesh(mesh, TWOD) + call test_semi_imp_alg( modeldb, & mesh, & twod_mesh ) @@ -349,6 +319,9 @@ subroutine run_rhs_alg(modeldb) type(modeldb_type), intent(inout) :: modeldb + mesh => mesh_collection%get_mesh(prime_mesh_name) + twod_mesh => mesh_collection%get_mesh(mesh, TWOD) + call test_rhs_alg( modeldb, & mesh, & twod_mesh ) @@ -364,6 +337,9 @@ subroutine run_rhs_project_eos(modeldb) type(modeldb_type), intent(inout) :: modeldb + mesh => mesh_collection%get_mesh(prime_mesh_name) + twod_mesh => mesh_collection%get_mesh(mesh, TWOD) + call test_rhs_project_eos( modeldb, & mesh, & twod_mesh ) @@ -379,31 +355,13 @@ subroutine run_rhs_sample_eos(modeldb) type(modeldb_type), intent(inout) :: modeldb + mesh => mesh_collection%get_mesh(prime_mesh_name) + twod_mesh => mesh_collection%get_mesh(mesh, TWOD) + call test_rhs_sample_eos( modeldb, & mesh, & twod_mesh ) end subroutine run_rhs_sample_eos - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !>@brief Tidies up after a run. - !>@param [in] program_name An identifier given to the model being run - !>@param [in,out] modeldb The structure that holds model state - subroutine finalise( program_name, modeldb ) - - implicit none - - character(*), intent(in) :: program_name - type(modeldb_type), intent(inout) :: modeldb - - call log_event( 'Finalising '//program_name//' ...', LOG_LEVEL_ALWAYS ) - - ! Destroy the fields stored in model_data - call finalise_model_data( modeldb ) - - ! Finalise infrastructure and constants - call finalise_infrastructure(modeldb) - - end subroutine finalise - end module tl_test_driver_mod diff --git a/science/linear/integration-test/tl_test/tl_test_hydrostatic_mod.x90 b/science/linear/integration-test/tl_test/tl_test_hydrostatic_mod.x90 index fbcb7baa..5aa4f198 100644 --- a/science/linear/integration-test/tl_test/tl_test_hydrostatic_mod.x90 +++ b/science/linear/integration-test/tl_test/tl_test_hydrostatic_mod.x90 @@ -189,8 +189,6 @@ module tl_test_hydrostatic_mod inc_X_minus_Y( diff, p_rhs_u ), & X_innerproduct_X( norm_diff, diff ) ) - norm_diff = sqrt(norm_diff) - write( log_scratch_space, '(A, E32.12, A, E32.12)' ) & 'gamma = ', gamma, ' norm = ', norm_diff diff --git a/science/linear/integration-test/tl_test/tl_test_kinetic_energy_gradient_mod.x90 b/science/linear/integration-test/tl_test/tl_test_kinetic_energy_gradient_mod.x90 index 2138662a..f2910b1b 100644 --- a/science/linear/integration-test/tl_test/tl_test_kinetic_energy_gradient_mod.x90 +++ b/science/linear/integration-test/tl_test/tl_test_kinetic_energy_gradient_mod.x90 @@ -135,7 +135,6 @@ module tl_test_kinetic_energy_gradient_mod call invoke( X_minus_Y( diff, n2_rhs_u, n1_rhs_u ), & inc_X_minus_Y( diff, p_rhs_u ), & X_innerproduct_X( norm_diff, diff ) ) - norm_diff = sqrt(norm_diff) write( log_scratch_space, '(A, E32.12, A, E32.12)' ) & 'gamma = ', gamma, ' norm = ', norm_diff diff --git a/science/linear/integration-test/tl_test/tl_test_pressure_grad_bd_mod.x90 b/science/linear/integration-test/tl_test/tl_test_pressure_grad_bd_mod.x90 index 8bb6117a..5a9f0210 100644 --- a/science/linear/integration-test/tl_test/tl_test_pressure_grad_bd_mod.x90 +++ b/science/linear/integration-test/tl_test/tl_test_pressure_grad_bd_mod.x90 @@ -213,8 +213,6 @@ module tl_test_pressure_grad_bd_mod inc_X_minus_Y( diff, p_rhs_u ), & X_innerproduct_X( norm_diff, diff ) ) - norm_diff = sqrt(norm_diff) - write( log_scratch_space, '(A, E32.12, A, E32.12)' ) & 'gamma = ', gamma, ' norm = ', norm_diff diff --git a/science/linear/integration-test/tl_test/tl_test_project_eos_pressure_mod.x90 b/science/linear/integration-test/tl_test/tl_test_project_eos_pressure_mod.x90 index 01ef9b69..9fc26e06 100644 --- a/science/linear/integration-test/tl_test/tl_test_project_eos_pressure_mod.x90 +++ b/science/linear/integration-test/tl_test/tl_test_project_eos_pressure_mod.x90 @@ -195,7 +195,6 @@ module tl_test_project_eos_pressure_mod call invoke( X_minus_Y( diff, n2_exner, n1_exner ), & inc_X_minus_Y( diff, p_exner ), & X_innerproduct_X( norm_diff, diff ) ) - norm_diff = sqrt(norm_diff) write( log_scratch_space, '(A, E32.12, A, E32.12)' ) & 'gamma = ', gamma , ' norm = ', norm_diff diff --git a/science/linear/integration-test/tl_test/tl_test_rhs_alg_mod.x90 b/science/linear/integration-test/tl_test/tl_test_rhs_alg_mod.x90 index 1c483372..6207afcd 100644 --- a/science/linear/integration-test/tl_test/tl_test_rhs_alg_mod.x90 +++ b/science/linear/integration-test/tl_test/tl_test_rhs_alg_mod.x90 @@ -18,7 +18,7 @@ module tl_test_rhs_alg_mod use function_space_mod, only: function_space_type use derived_config_mod, only: bundle_size use moist_dyn_mod, only: num_moist_factors - use sci_assign_field_random_kernel_mod, only: assign_field_random_kernel_type + use mr_indices_mod, only: nummr use field_indices_mod, only: igh_u, igh_t, igh_d, igh_p use sci_field_bundle_builtins_mod, only: clone_bundle, & add_bundle, & @@ -29,7 +29,9 @@ module tl_test_rhs_alg_mod use tl_rhs_alg_mod, only: tl_rhs_alg use log_mod, only: log_event, LOG_LEVEL_INFO, & LOG_LEVEL_ERROR - use tl_test_convergence_rate_check, only: convergence_rate_check + use moist_dyn_factors_alg_mod, only: moist_dyn_factors_alg + use tl_moist_dyn_factors_alg_mod, only: tl_moist_dyn_factors_alg + use tl_test_convergence_rate_check, only: array_convergence_rate_check implicit none @@ -56,8 +58,8 @@ module tl_test_rhs_alg_mod character(str_def) :: label = "rhs_alg" - type(field_collection_type ), pointer :: ls_fields - + type( field_collection_type ), pointer :: ls_fields + type( field_collection_type ), pointer :: prognostic_fields => null() type( field_type ) :: state(bundle_size) type( field_type ) :: ls_state(bundle_size) @@ -73,37 +75,62 @@ module tl_test_rhs_alg_mod type(field_type), pointer :: ls_theta => null() type(field_type), pointer :: ls_exner => null() type(field_type), pointer :: ls_moist_dyn(:) => null() + type(field_type), pointer :: ls_mr(:) => null() + type(field_type), pointer :: u => null() + type(field_type), pointer :: rho => null() + type(field_type), pointer :: theta => null() + type(field_type), pointer :: exner => null() + type(field_type), pointer :: moist_dyn(:) => null() + type(field_type), pointer :: mr(:) => null() type(field_type), dimension(num_moist_factors) :: p_moist_dyn type(field_type), dimension(num_moist_factors) :: n_moist_dyn type(field_type), dimension(num_moist_factors) :: r_moist_dyn type(field_collection_type), pointer :: moisture_fields => null() + type(field_array_type), pointer :: mr_array => null() + type(field_array_type), pointer :: moist_dyn_array => null() + type(field_array_type), pointer :: ls_mr_array => null() type(field_array_type), pointer :: ls_moist_dyn_array => null() real(r_def) :: gamma_u, gamma_rho, gamma_exner, gamma_theta, & gamma_moist_dyn - real(r_def) :: norm_u, norm_rho, norm_exner, norm_theta, norm_moist_dyn - real(r_def) :: norm_moist_dyn_tmp - real(r_def) :: norm_diff, norm_diff_prev + real(r_def) :: norm_u, norm_exner + integer(i_def), parameter :: n_variables = 2 + real(r_def) :: array_norm(n_variables), array_norm_prev(n_variables) + character(str_def) :: array_names(n_variables) - real(r_def), parameter :: tol = 1.e-2_r_def + real(r_def), parameter :: tol = 5.e-2_r_def + real(r_def), parameter :: indiv_tol = 5.e-1_r_def integer :: i, n call log_event( "TL Test: " // trim(label), & LOG_LEVEL_INFO ) + prognostic_fields => modeldb%fields%get_field_collection( & + "prognostic_fields") ls_fields => modeldb%fields%get_field_collection("ls_fields") moisture_fields => modeldb%fields%get_field_collection("moisture_fields") + call moisture_fields%get_field("mr", mr_array) + call moisture_fields%get_field("moist_dyn", moist_dyn_array) + mr => mr_array%bundle + moist_dyn => moist_dyn_array%bundle + call moisture_fields%get_field("ls_mr", ls_mr_array) call moisture_fields%get_field("ls_moist_dyn", ls_moist_dyn_array) + ls_mr => ls_mr_array%bundle ls_moist_dyn => ls_moist_dyn_array%bundle - ! Input + ! Linearisation state call ls_fields%get_field('ls_u', ls_u) call ls_fields%get_field('ls_rho', ls_rho) call ls_fields%get_field('ls_theta', ls_theta) call ls_fields%get_field('ls_exner', ls_exner) + ! Perturbation + call prognostic_fields%get_field('u', u) + call prognostic_fields%get_field('rho', rho) + call prognostic_fields%get_field('theta', theta) + call prognostic_fields%get_field('exner', exner) call ls_state(igh_u)%initialise( vector_space = ls_u%get_function_space() ) call ls_state(igh_t)%initialise( vector_space = ls_theta%get_function_space() ) @@ -127,24 +154,27 @@ module tl_test_rhs_alg_mod setval_X(ls_state(igh_d), ls_rho ), & setval_X(ls_state(igh_p), ls_exner) ) - call invoke( assign_field_random_kernel_type( random(igh_u), 1.0_r_def ) , & - assign_field_random_kernel_type( random(igh_d), 1.0_r_def ) , & - assign_field_random_kernel_type( random(igh_t), 1.0_r_def ) , & - assign_field_random_kernel_type( random(igh_p), 1.0_r_def ) ) + call moist_dyn_factors_alg(ls_moist_dyn, ls_mr) - do i = 1, num_moist_factors - call invoke( assign_field_random_kernel_type( r_moist_dyn(i), 1.0_r_def ) ) - enddo + ! Set the 'random data' to be the perturbation + call invoke( name = "copy_fields_to_state", & + setval_X(random(igh_u), u ), & + setval_X(random(igh_p), exner ), & + setval_X(random(igh_t), theta), & + setval_X(random(igh_d), rho ) ) - gamma_u = 1.e2_r_def - gamma_theta = 1.e2_r_def - gamma_rho = 1.e2_r_def - gamma_exner = 1.e1_r_def - gamma_moist_dyn = 1.e-1_r_def + call tl_moist_dyn_factors_alg(r_moist_dyn, mr ) + + gamma_u = 1.0_r_def + gamma_theta = 1.0_r_def + gamma_rho = 1.0_r_def + gamma_exner = 1.0_r_def + gamma_moist_dyn = 1.0_r_def call set_bundle_scalar( 0.0_r_def, n1_rhs, bundle_size ) call rhs_alg( n1_rhs, dt, ls_state, ls_state, ls_moist_dyn, & - .true., .true., .false., modeldb%clock ) + compute_eos=.true., compute_rhs_t_d=.true., & + dlayer_rhs=.true., model_clock=modeldb%clock ) do n=1,2 gamma_u = gamma_u / 2.0_r_def @@ -166,42 +196,45 @@ module tl_test_rhs_alg_mod do i = 1, num_moist_factors call invoke( & a_times_X( p_moist_dyn(i), gamma_moist_dyn, r_moist_dyn(i) ), & - setval_X( n_moist_dyn(i), ls_moist_dyn(i) ), & + setval_X( n_moist_dyn(i), ls_moist_dyn(i) ), & inc_X_plus_Y( n_moist_dyn(i), p_moist_dyn(i) ) ) enddo call rhs_alg( n2_rhs, dt, state, state, n_moist_dyn, & - .true., .true., .false., modeldb%clock ) + compute_eos=.true., compute_rhs_t_d=.true., & + dlayer_rhs=.true., model_clock=modeldb%clock ) call tl_rhs_alg(p_rhs, dt, p_state, p_state, p_moist_dyn, & ls_state, ls_moist_dyn, & - .true., .false., modeldb%clock) + compute_eos=.true., dlayer_rhs=.true., & + model_clock=modeldb%clock) ! diff = n2_rhs - n1_rhs call minus_bundle( n2_rhs, n1_rhs, diff, bundle_size ) call invoke( & inc_X_minus_Y( diff(igh_u), p_rhs(igh_u) ), & - inc_X_minus_Y( diff(igh_t), p_rhs(igh_t) ), & - inc_X_minus_Y( diff(igh_d), p_rhs(igh_d) ), & inc_X_minus_Y( diff(igh_p), p_rhs(igh_p) ) ) call invoke( X_innerproduct_X( norm_u, diff(igh_u) ) , & - X_innerproduct_X( norm_rho, diff(igh_d) ) , & - X_innerproduct_X( norm_theta, diff(igh_t) ) , & X_innerproduct_X( norm_exner, diff(igh_p) ) ) - print*, norm_u, norm_rho, norm_theta, norm_exner - norm_diff = norm_u + norm_rho + norm_theta + norm_exner - norm_diff = sqrt(norm_diff) - print*,'norm', norm_diff + ! The equations for rho and theta are already linear + ! so only evaluate u and exner + array_norm(1) = norm_u + array_norm(2) = norm_exner + array_names(1) = 'norm_u' + array_names(2) = 'norm_exner' if (n == 2) then - call convergence_rate_check( norm_diff, norm_diff_prev, & - label, tol=tol ) - endif + call array_convergence_rate_check( & + array_norm, array_norm_prev, & + array_names, n_variables, label, tol=tol, & + indiv_tol=indiv_tol ) + end if + + array_norm_prev = array_norm - norm_diff_prev = norm_diff enddo end subroutine test_rhs_alg diff --git a/science/linear/integration-test/tl_test/tl_test_rhs_project_eos_mod.x90 b/science/linear/integration-test/tl_test/tl_test_rhs_project_eos_mod.x90 index 9504122f..8b39ecb5 100644 --- a/science/linear/integration-test/tl_test/tl_test_rhs_project_eos_mod.x90 +++ b/science/linear/integration-test/tl_test/tl_test_rhs_project_eos_mod.x90 @@ -209,7 +209,6 @@ contains call invoke( X_minus_Y( diff, n2_rhs, n1_rhs ), & inc_X_minus_Y( diff, p_rhs ), & X_innerproduct_X( norm_diff, diff ) ) - norm_diff = sqrt(norm_diff) if (n == 2) then call convergence_rate_check( norm_diff, norm_diff_prev, & diff --git a/science/linear/integration-test/tl_test/tl_test_rhs_sample_eos_mod.x90 b/science/linear/integration-test/tl_test/tl_test_rhs_sample_eos_mod.x90 index 49d69273..0d1cd26d 100644 --- a/science/linear/integration-test/tl_test/tl_test_rhs_sample_eos_mod.x90 +++ b/science/linear/integration-test/tl_test/tl_test_rhs_sample_eos_mod.x90 @@ -195,7 +195,6 @@ contains call invoke( X_minus_Y( diff, n2_rhs, n1_rhs ), & inc_X_minus_Y( diff, p_rhs ), & X_innerproduct_X( norm_diff, diff ) ) - norm_diff = sqrt(norm_diff) if (n == 2) then call convergence_rate_check( norm_diff, norm_diff_prev, & diff --git a/science/linear/integration-test/tl_test/tl_test_rk_alg_mod.x90 b/science/linear/integration-test/tl_test/tl_test_rk_alg_mod.x90 index 9b99f0fc..d16f4c14 100644 --- a/science/linear/integration-test/tl_test/tl_test_rk_alg_mod.x90 +++ b/science/linear/integration-test/tl_test/tl_test_rk_alg_mod.x90 @@ -251,6 +251,7 @@ module tl_test_rk_alg_mod ! Evolve the perturbation fields with the linear timestep. This ! gives p_x=L(gamma r_x) + call tl_rk_alg_final() call tl_rk_alg_init(mesh, p_u, p_rho, p_theta, p_exner, & ls_u, ls_rho, ls_theta, ls_exner) @@ -300,7 +301,6 @@ module tl_test_rk_alg_mod call log_event( log_scratch_space, LOG_LEVEL_INFO ) norm_diff = norm_rho + norm_theta + norm_exner + norm_u + norm_moist_dyn - norm_diff = sqrt(norm_diff) if (n == 2) then ! Compare results from first and second iteration diff --git a/science/linear/integration-test/tl_test/tl_test_sample_eos_pressure_mod.x90 b/science/linear/integration-test/tl_test/tl_test_sample_eos_pressure_mod.x90 index b6f56add..1915214c 100644 --- a/science/linear/integration-test/tl_test/tl_test_sample_eos_pressure_mod.x90 +++ b/science/linear/integration-test/tl_test/tl_test_sample_eos_pressure_mod.x90 @@ -175,7 +175,6 @@ module tl_test_sample_eos_pressure_mod call invoke( X_minus_Y( diff, n2_exner, n1_exner ), & inc_X_minus_Y( diff, p_exner ), & X_innerproduct_X( norm_diff, diff ) ) - norm_diff = sqrt(norm_diff) write( log_scratch_space, '(A, E32.12, A, E32.12)' ) & 'gamma = ', gamma , ' norm = ', norm_diff diff --git a/science/linear/integration-test/tl_test/tl_test_semi_imp_alg_mod.x90 b/science/linear/integration-test/tl_test/tl_test_semi_imp_alg_mod.x90 index 2b060eac..30bf004e 100644 --- a/science/linear/integration-test/tl_test/tl_test_semi_imp_alg_mod.x90 +++ b/science/linear/integration-test/tl_test/tl_test_semi_imp_alg_mod.x90 @@ -10,8 +10,6 @@ !! implicit algorithm, using the Taylor remainder convergence test. module tl_test_semi_imp_alg_mod - use sci_assign_field_random_kernel_mod, & - only: assign_field_random_kernel_type use constants_mod, only: i_def, r_def, str_def use field_array_mod, only: field_array_type use field_mod, only: field_type @@ -39,7 +37,7 @@ module tl_test_semi_imp_alg_mod use timestep_method_mod, only: timestep_method_type use moist_dyn_factors_alg_mod, only: moist_dyn_factors_alg use tl_moist_dyn_factors_alg_mod, only: tl_moist_dyn_factors_alg - use tl_test_convergence_rate_check, only: convergence_rate_check + use tl_test_convergence_rate_check, only: array_convergence_rate_check implicit none @@ -128,9 +126,12 @@ module tl_test_semi_imp_alg_mod real(r_def) :: gamma_mr real(r_def) :: norm_u, norm_rho, norm_exner, norm_theta real(r_def) :: norm_mr, norm_mr_tmp - real(r_def) :: norm_diff, norm_diff_prev + integer(i_def), parameter :: n_variables = 5 + real(r_def) :: array_norm(n_variables), array_norm_prev(n_variables) + character(str_def) :: array_names(n_variables) - real(r_def), parameter :: tol = 5.e-1_r_def + real(r_def), parameter :: tol = 5.e-2_r_def + real(r_def), parameter :: indiv_tol = 5.e-1_r_def real(r_def), parameter :: dtemp_encorr = 0.0_r_def integer :: n, i @@ -173,12 +174,14 @@ module tl_test_semi_imp_alg_mod aerosol_fields => modeldb%fields%get_field_collection("aerosol_fields") stph_fields => modeldb%fields%get_field_collection("stph_fields") lbc_fields => modeldb%fields%get_field_collection("lbc_fields") - ! Input + + ! Linearisation State call ls_fields%get_field('ls_u', ls_u) call ls_fields%get_field('ls_rho', ls_rho) call ls_fields%get_field('ls_theta', ls_theta) call ls_fields%get_field('ls_exner', ls_exner) + ! Perturbation call prognostic_fields%get_field('u', u) call prognostic_fields%get_field('rho', rho) call prognostic_fields%get_field('theta', theta) @@ -228,15 +231,29 @@ module tl_test_semi_imp_alg_mod call clone_bundle(mr, b_mr, nummr) call clone_bundle(mr, diff_mr, nummr) - call invoke( assign_field_random_kernel_type( r_u, 1.0_r_def ) , & - assign_field_random_kernel_type( r_rho, 1.0_r_def ) , & - assign_field_random_kernel_type( r_theta, 1.0_r_def ) , & - assign_field_random_kernel_type( r_exner, 1.0_r_def ) ) - + ! Set the 'random data' to be the perturbation + call invoke( name = "copy_fields_to_state", & + setval_X(r_u, u ), & + setval_X(r_exner, exner ), & + setval_X(r_theta, theta), & + setval_X(r_rho, rho ) ) + + ! Overwrite the moisture fields with potential temperature + ! + ! Note: The perturbation fields are read in from file. However the + ! moisture perturbation fields have values close to zero (especially + ! at upper levels) and this gives an inaccurate linearisation test. Using + ! potential temperature is an easy fix to create a realistic 'random' + ! field using real data that has mostly non-zero values. do i = 1, nummr - call invoke( assign_field_random_kernel_type( r_mr(i), 1.0_r_def ) ) + call invoke( & + setval_X( r_mr(i), theta ), & + inc_a_times_X( 0.000001_r_def, r_mr(i) ), & + setval_X( ls_mr(i), ls_theta ), & + inc_a_times_X( 0.000001_r_def, ls_mr(i) )) end do + ! Set the prognostic fields to be the linearisation state call invoke( setval_X( u, ls_u ), & setval_X( theta, ls_theta ), & setval_X( rho, ls_rho ), & @@ -247,7 +264,9 @@ module tl_test_semi_imp_alg_mod end do do i = 1, nummr call invoke( setval_X( mr(i), ls_mr(i) ) ) - end do + end do + + call moist_dyn_factors_alg(moist_dyn, mr) ! Clean up solvers left over from previous tests, as semi-implicit ! method will create new solvers which cause double-allocates @@ -274,15 +293,15 @@ module tl_test_semi_imp_alg_mod call invoke( setval_X( b_mr(i), mr(i) ) ) end do - gamma_u = 1.e4_r_def - gamma_theta = 1.e2_r_def - gamma_rho = 8.0e-1_r_def - gamma_exner = 6.5e-1_r_def + gamma_u = 6.5_r_def + gamma_theta = 6.5_r_def + gamma_rho = 6.5_r_def + gamma_exner = 6.5_r_def if ( moisture_formulation == moisture_formulation_dry ) then gamma_mr = 0.0_r_def else - gamma_mr = 1.e-2_r_def + gamma_mr = 6.5_r_def end if do n=1,2 @@ -328,6 +347,8 @@ module tl_test_semi_imp_alg_mod ! Evolve the perturbation fields with the linear timestep. This ! gives p_x=L(gamma r_x) + call final_si_operators() + call tl_semi_implicit_alg_final() call tl_semi_implicit_alg_init(mesh, p_u, p_rho, p_theta, p_exner, & p_mr, ls_u, ls_rho, ls_theta, ls_exner, & ls_mr, ls_moist_dyn) @@ -365,35 +386,24 @@ module tl_test_semi_imp_alg_mod norm_mr = norm_mr + norm_mr_tmp end do - write( log_scratch_space, & - '(A, E32.12, A, E32.12 )' ) & - ' norm_u = ' , norm_u, & - ' norm_rho = ' , norm_rho - - call log_event( log_scratch_space, LOG_LEVEL_INFO ) - - write( log_scratch_space, & - '(A, E32.12, A, E32.12 )' ) & - ' norm_mr = ' , norm_mr, & - ' norm_exner = ' , norm_exner - - call log_event( log_scratch_space, LOG_LEVEL_INFO ) - - write( log_scratch_space, & - '(A, E32.12 )' ) & - ' norm_theta = ' , norm_theta - - call log_event( log_scratch_space, LOG_LEVEL_INFO ) - - norm_diff = norm_rho + norm_theta + norm_exner + norm_u + norm_mr - norm_diff = sqrt(norm_diff) + array_norm(1) = norm_u + array_norm(2) = norm_rho + array_norm(3) = norm_theta + array_norm(4) = norm_mr + array_norm(5) = norm_exner + array_names(1) = 'norm_u' + array_names(2) = 'norm_rho' + array_names(3) = 'norm_theta' + array_names(4) = 'norm_mr' + array_names(5) = 'norm_exner' if (n == 2) then - call convergence_rate_check( norm_diff, norm_diff_prev, & - label, tol=tol ) + call array_convergence_rate_check( array_norm, array_norm_prev, & + array_names, n_variables, label, tol=tol, indiv_tol=indiv_tol ) end if - norm_diff_prev = norm_diff + array_norm_prev = array_norm + end do end subroutine test_semi_imp_alg diff --git a/science/linear/integration-test/tl_test/tl_test_timesteps_alg_mod.x90 b/science/linear/integration-test/tl_test/tl_test_timesteps_alg_mod.x90 index 65759abf..686e0cad 100644 --- a/science/linear/integration-test/tl_test/tl_test_timesteps_alg_mod.x90 +++ b/science/linear/integration-test/tl_test/tl_test_timesteps_alg_mod.x90 @@ -44,7 +44,7 @@ module tl_test_timesteps_alg_mod method, & method_rk use time_config_mod, only: timestep_start, timestep_end - use tl_test_convergence_rate_check, only: convergence_rate_check + use tl_test_convergence_rate_check, only: array_convergence_rate_check use validity_test_config_mod, only: number_gamma_values, & update_ls_frequency use moist_dyn_factors_alg_mod, only: moist_dyn_factors_alg @@ -52,7 +52,6 @@ module tl_test_timesteps_alg_mod implicit none - private convergence_pass_string public test_timesteps contains @@ -120,8 +119,12 @@ module tl_test_timesteps_alg_mod real(r_def) :: denom_u, denom_rho, denom_exner, denom_theta, denom_mr real(r_def) :: denom_mr_tmp real(r_def) :: denom_total + integer(i_def), parameter :: n_variables = 5 + real(r_def) :: array_norm(n_variables), array_norm_prev(n_variables) + character(str_def) :: array_names(n_variables) - real(r_def), parameter :: tol = 9.e-1_r_def + real(r_def), parameter :: tol = 3.e-1_r_def + real(r_def), parameter :: indiv_tol = 1.2_r_def integer(i_def) :: n, i, t character(len=4) :: pass_str @@ -212,24 +215,32 @@ module tl_test_timesteps_alg_mod function_space_collection%get_fs(mesh, element_order_h, element_order_v, Wtheta) ) call clone_bundle(mr, diff_mr, nummr) - call invoke( assign_field_random_kernel_type( r_u, 1.0_r_def ) , & - assign_field_random_kernel_type( r_rho, 1.0_r_def ) , & - assign_field_random_kernel_type( r_theta, 1.0_r_def ) , & - assign_field_random_kernel_type( r_exner, 1.0_r_def ) ) + ! Set the 'random data' to be the perturbation + call invoke( name = "copy_fields_to_state", & + setval_X(r_u, u ), & + setval_X(r_exner, exner ), & + setval_X(r_theta, theta), & + setval_X(r_rho, rho ) ) + ! Overwrite the moisture fields with potential temperature do i = 1, nummr - call invoke( assign_field_random_kernel_type( r_mr(i), 1.0_r_def ) ) + call invoke( & + setval_X( r_mr(i), theta ), & + inc_a_times_X( 0.000001_r_def, r_mr(i) ), & + setval_X( ls_mr(i), ls_theta ), & + inc_a_times_X( 0.000001_r_def, ls_mr(i) )) end do + call moist_dyn_factors_alg( ls_moist_dyn, ls_mr ) - gamma_u = 1.e4_r_def - gamma_theta = 1.e2_r_def - gamma_rho = 1.e0_r_def - gamma_exner = 6.5e-1_r_def + gamma_u = 6.5_r_def + gamma_theta = 6.5_r_def + gamma_rho = 6.5_r_def + gamma_exner = 6.5_r_def - if ( use_moisture ) then - gamma_mr = 1.e-1_r_def + if ( use_moisture ) then + gamma_mr = 6.5_r_def else - gamma_mr = 0.0_r_def + gamma_mr = 0._r_def end if do n = 1, number_gamma_values @@ -277,6 +288,10 @@ module tl_test_timesteps_alg_mod !------------------------------------------------------------------------ ! Nonlinear model run with perturbed state !------------------------------------------------------------------------ + write(log_scratch_space,'(A, I2)') & + "Perturbed Nonlinear, timestep: ", t + call log_event(log_scratch_space, LOG_LEVEL_INFO) + call invoke( setval_X( u, n1_u ), & setval_X( theta, n1_theta ), & setval_X( rho, n1_rho ), & @@ -314,6 +329,10 @@ module tl_test_timesteps_alg_mod !------------------------------------------------------------------------ ! Linear model run !------------------------------------------------------------------------ + write(log_scratch_space,'(A, I2)') & + "Linear, timestep: ", t + call log_event(log_scratch_space, LOG_LEVEL_INFO) + call invoke( setval_X( u, p_u ), & setval_X( theta, p_theta ), & setval_X( rho, p_rho ), & @@ -337,6 +356,7 @@ module tl_test_timesteps_alg_mod call moist_dyn_factors_alg( ls_moist_dyn, ls_mr ) endif + call finalise_linear_model() call initialise_linear_model( mesh, & modeldb ) @@ -359,6 +379,10 @@ module tl_test_timesteps_alg_mod !------------------------------------------------------------------------ ! Nonlinear model run with un-perturbed state !------------------------------------------------------------------------ + write(log_scratch_space,'(A, I2)') & + "Unperturbed NonLinear, timestep: ", t + call log_event(log_scratch_space, LOG_LEVEL_INFO) + call invoke( setval_X( u, n2_u ), & setval_X( theta, n2_theta ), & setval_X( rho, n2_rho ), & @@ -386,6 +410,10 @@ module tl_test_timesteps_alg_mod setval_X( n2_rho, rho ), & setval_X( n2_exner, exner ) ) + do i = 1, nummr + call invoke( setval_X( n2_mr(i), mr(i) ) ) + end do + do i = 1, nummr call invoke( setval_X( ls_mr(i), mr(i) ) ) end do @@ -417,160 +445,67 @@ module tl_test_timesteps_alg_mod norm_mr = norm_mr + norm_mr_tmp end do - norm_diff = norm_rho + norm_theta + norm_exner + norm_u + norm_mr - norm_diff = sqrt( norm_diff ) - norm_u = sqrt( norm_u ) - norm_rho= sqrt( norm_rho ) - norm_theta = sqrt( norm_theta ) - norm_exner = sqrt( norm_exner ) - norm_mr = sqrt( norm_mr ) - - call invoke( X_innerproduct_X( denom_u, p_u ) ) - call invoke( X_innerproduct_X( denom_rho, p_rho ) ) - call invoke( X_innerproduct_X( denom_exner, p_exner ) ) - call invoke( X_innerproduct_X( denom_theta, p_theta ) ) + array_norm(1) = norm_u + array_norm(2) = norm_rho + array_norm(3) = norm_theta + array_norm(4) = norm_mr + array_norm(5) = norm_exner + array_names(1) = 'norm_u' + array_names(2) = 'norm_rho' + array_names(3) = 'norm_theta' + array_names(4) = 'norm_mr' + array_names(5) = 'norm_exner' - denom_mr=0.0_r_def - do i = 1, nummr - call invoke( & - X_innerproduct_X( denom_mr_tmp, p_mr(i) ) ) - denom_mr = denom_mr + denom_mr_tmp - end do + if (n == 2) then + call array_convergence_rate_check( array_norm, array_norm_prev, & + array_names, n_variables, label, tol=tol, indiv_tol=indiv_tol ) + end if - denom_total = denom_rho + denom_theta + denom_exner + denom_u + denom_mr - denom_total = sqrt( denom_total ) - denom_u = sqrt( denom_u ) - denom_rho = sqrt( denom_rho ) - denom_theta = sqrt( denom_theta ) - denom_exner = sqrt( denom_exner ) - denom_mr = sqrt( denom_mr ) + array_norm_prev = array_norm !---------------------------------------------------------------------------- ! Print values out to enable plotting a graph !---------------------------------------------------------------------------- - write( log_scratch_space, & - '( A, E32.12, A, E32.12 )' ) & - ' gamma_total = ' , gamma_u, & - ' norm = ' , norm_diff / denom_total - - call log_event( log_scratch_space, LOG_LEVEL_INFO ) write( log_scratch_space, & '(A, E32.12, A, E32.12)' ) & - ' gamma_rho = ' , gamma_u, & - ' norm = ' , norm_rho / denom_rho + ' gamma_rho = ' , gamma_rho, & + ' norm = ' , norm_rho call log_event( log_scratch_space, LOG_LEVEL_INFO ) write( log_scratch_space, & '(A, E32.12, A, E32.12)' ) & - ' gamma_exner = ' , gamma_u, & - ' norm = ' , norm_exner / denom_exner + ' gamma_exner = ' , gamma_exner, & + ' norm = ' , norm_exner call log_event( log_scratch_space, LOG_LEVEL_INFO ) write( log_scratch_space, & '(A, E32.12, A, E32.12)' ) & - ' gamma_theta = ' , gamma_u, & - ' norm = ' , norm_theta / denom_theta + ' gamma_theta = ' , gamma_theta, & + ' norm = ' , norm_theta call log_event( log_scratch_space, LOG_LEVEL_INFO ) write( log_scratch_space, & '(A, E32.12, A, E32.12)' ) & ' gamma_u = ' , gamma_u, & - ' norm = ' , norm_u / denom_u + ' norm = ' , norm_u call log_event( log_scratch_space, LOG_LEVEL_INFO ) write( log_scratch_space, & '(A, E32.12, A, E32.12)' ) & - ' gamma_mr = ' , gamma_u, & - ' norm = ' , norm_mr / denom_mr + ' gamma_mr = ' , gamma_mr, & + ' norm = ' , norm_mr call log_event( log_scratch_space, LOG_LEVEL_INFO ) -!---------------------------------------------------------------------------- -! Test that gives a PASS or FAIL -!---------------------------------------------------------------------------- - if (n == 2) then - - pass_str = "PASS" - - write( log_scratch_space, '(A)' ) & - "TL Test: " // trim(label) - call log_event( log_scratch_space, LOG_LEVEL_INFO ) - - write( log_scratch_space, '(A)' ) "Total " - call log_event( log_scratch_space, LOG_LEVEL_INFO ) - call convergence_pass_string( norm_diff, norm_diff_prev, pass_str, tol ) - - write( log_scratch_space, '(A)' ) "Theta " - call log_event( log_scratch_space, LOG_LEVEL_INFO ) - call convergence_pass_string( norm_theta, norm_theta_prev, pass_str, tol ) - - write( log_scratch_space, '(A)' ) "Rho " - call log_event( log_scratch_space, LOG_LEVEL_INFO ) - call convergence_pass_string( norm_rho, norm_rho_prev, pass_str, tol ) - - write( log_scratch_space, '(A)' ) "Exner " - call log_event( log_scratch_space, LOG_LEVEL_INFO ) - call convergence_pass_string( norm_exner, norm_exner_prev, pass_str, tol ) - - write( log_scratch_space, '(A)' ) "u " - call log_event( log_scratch_space, LOG_LEVEL_INFO ) - call convergence_pass_string( norm_u, norm_u_prev, pass_str, tol ) - - if ( use_moisture ) then - write( log_scratch_space, '(A)' ) "mr " - call log_event( log_scratch_space, LOG_LEVEL_INFO ) - call convergence_pass_string( norm_mr, norm_mr_prev, pass_str, tol ) - end if - - write( log_scratch_space,'(" test",A32," : ",A4)' ) trim(label), pass_str - call log_event( log_scratch_space, LOG_LEVEL_INFO ) - - end if - - norm_diff_prev = norm_diff - norm_rho_prev = norm_rho - norm_theta_prev = norm_theta - norm_exner_prev = norm_exner - norm_u_prev = norm_u - norm_mr_prev = norm_mr - end do ! Loop over values of gamma end subroutine test_timesteps - !> @brief Test the convergence rate. - !> @details If the convergence rate is not close to 4, within the - !! specified tolerance, then change the pass string to FAIL. - !> @param[in] norm_diff Current norm - !> @param[in] norm_diff_prev Previous norm - !> @param[inout] pass_str Pass string (either PASS or FAIL) - !> @param[in] tol Test Tolerance - subroutine convergence_pass_string( norm_diff, norm_diff_prev, pass_str, tol ) - implicit none - - real(r_def), intent(in) :: norm_diff, norm_diff_prev - character(len=4), intent(inout) :: pass_str - real(r_def), intent(in) :: tol - - real(r_def) :: conv_rate - conv_rate = norm_diff_prev / norm_diff - - conv_rate = norm_diff_prev / norm_diff - write( log_scratch_space, '(A, E16.8)') & - "Convergence rate: ", conv_rate - call log_event( log_scratch_space, LOG_LEVEL_INFO ) - - if ( abs( conv_rate - 4.0_r_def ) >= tol ) then - pass_str = "FAIL" - endif - - end subroutine convergence_pass_string - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Converts a string to a timestep number. !> diff --git a/science/linear/integration-test/tl_test/tl_test_timesteps_random_alg_mod.x90 b/science/linear/integration-test/tl_test/tl_test_timesteps_random_alg_mod.x90 new file mode 100644 index 00000000..6ea6ff29 --- /dev/null +++ b/science/linear/integration-test/tl_test/tl_test_timesteps_random_alg_mod.x90 @@ -0,0 +1,597 @@ +!----------------------------------------------------------------------------- +! (C) Crown copyright 2026 Met Office. All rights reserved. +! The file LICENCE, distributed with this code, contains details of the terms +! under which the code may be used. +!----------------------------------------------------------------------------- + +!>@brief The tangent linear test for multiple timesteps +!! (Taylor remainder convergence). +!>@details Test whether the tangent linear code is tangent linear +!! to the corresponding nonlinear code, for multiple timesteps, +!! using the Taylor remainder convergence test. +module tl_test_timesteps_random_alg_mod + use sci_assign_field_random_kernel_mod, & + only: assign_field_random_kernel_type + use constants_mod, only: i_timestep, i_def, r_def, l_def, & + str_def + use field_array_mod, only: field_array_type + use field_mod, only: field_type + use sci_field_bundle_builtins_mod, only: clone_bundle + use field_collection_mod, only: field_collection_type + use finite_element_config_mod, only: element_order_h, element_order_v + use formulation_config_mod, only: moisture_formulation, & + moisture_formulation_dry + use fs_continuity_mod, only: W2, W3, Wtheta + use function_space_collection_mod, only: function_space_collection + use gungho_step_mod, only: gungho_step + use gungho_model_mod, only: initialise_model, & + finalise_model + use driver_modeldb_mod, only: modeldb_type + use log_mod, only: log_event, & + log_scratch_space, & + LOG_LEVEL_INFO, & + LOG_LEVEL_ERROR + use linear_step_mod, only: linear_step + use linear_model_mod, only: initialise_linear_model, & + finalise_linear_model + use mesh_mod, only: mesh_type + use model_clock_mod, only: model_clock_type + use moist_dyn_mod, only: num_moist_factors + use mr_indices_mod, only: nummr + use semi_implicit_solver_alg_mod, only: semi_implicit_solver_alg_final + use si_operators_alg_mod, only: final_si_operators + use timestepping_config_mod, only: dt, & + method, & + method_rk + use time_config_mod, only: timestep_start, timestep_end + use tl_test_convergence_rate_check, only: convergence_rate_check + use validity_test_config_mod, only: number_gamma_values, & + update_ls_frequency + use moist_dyn_factors_alg_mod, only: moist_dyn_factors_alg + use tl_moist_dyn_factors_alg_mod, only: tl_moist_dyn_factors_alg + + implicit none + + private convergence_pass_string + public test_timesteps_random + + contains + + !> @brief Test the tangent linear runge-kutta algorithm. + !> @param[in] modeldb The working data set for a model run + !> @param[in] mesh The current 3d mesh + !> @param[in] twod_mesh The current 2d mesh + !> @param[in] model_clock Time within the model. + !> + subroutine test_timesteps_random( modeldb, & + mesh, & + twod_mesh, & + model_clock ) + + implicit none + + type(modeldb_type), target, intent(inout) :: modeldb + type(mesh_type), pointer, intent(in) :: mesh + type(mesh_type), pointer, intent(in) :: twod_mesh + class(model_clock_type), intent(in) :: model_clock + + character(str_def) :: label = "timesteps" + + type(field_collection_type), pointer :: prognostic_fields + type(field_collection_type), pointer :: ls_fields + + type(field_type), pointer :: ls_u + type(field_type), pointer :: ls_rho + type(field_type), pointer :: ls_theta + type(field_type), pointer :: ls_exner + type(field_type), pointer :: ls_moist_dyn(:) + type(field_type), pointer :: ls_mr(:) + + type(field_type), pointer :: u + type(field_type), pointer :: rho + type(field_type), pointer :: theta + type(field_type), pointer :: exner + type(field_type), pointer :: moist_dyn(:) + type(field_type), pointer :: mr(:) + + type(field_type), dimension(nummr) :: p_mr + type(field_type), dimension(nummr) :: n1_mr + type(field_type), dimension(nummr) :: n2_mr + type(field_type), dimension(nummr) :: r_mr + type(field_type), dimension(nummr) :: diff_mr + type(field_type) :: p_u, p_rho, p_theta, p_exner + type(field_type) :: n1_u, n1_rho, n1_theta, n1_exner + type(field_type) :: n2_u, n2_rho, n2_theta, n2_exner + type(field_type) :: r_u, r_rho, r_theta, r_exner + type(field_type) :: diff_u, diff_rho, diff_theta, diff_exner + + type(field_collection_type), pointer :: moisture_fields + type(field_array_type), pointer :: mr_array + type(field_array_type), pointer :: moist_dyn_array + type(field_array_type), pointer :: ls_mr_array + type(field_array_type), pointer :: ls_moist_dyn_array + + real(r_def) :: gamma_u, gamma_rho, gamma_exner, gamma_theta, gamma_mr + real(r_def) :: norm_u, norm_rho, norm_exner, norm_theta, norm_mr + real(r_def) :: norm_mr_tmp + real(r_def) :: norm_diff, norm_diff_prev + real(r_def) :: norm_u_prev, norm_rho_prev, norm_exner_prev + real(r_def) :: norm_theta_prev, norm_mr_prev + + real(r_def), parameter :: tol = 9.e-1_r_def + + integer(i_def) :: n, i, t + character(len=4) :: pass_str + logical(l_def) :: use_moisture + + call log_event( "TL Test: " // trim(label), & + LOG_LEVEL_INFO ) + + use_moisture = ( moisture_formulation /= moisture_formulation_dry ) + + if ( method == method_rk .and. use_moisture ) then + call log_event( & + 'Not possible to use Runge Kutta time stepping with moisture', & + log_level_error ) + end if + + prognostic_fields => modeldb%fields%get_field_collection( & + "prognostic_fields") + ls_fields => modeldb%fields%get_field_collection("ls_fields") + moisture_fields => modeldb%fields%get_field_collection("moisture_fields") + call moisture_fields%get_field("mr", mr_array) + call moisture_fields%get_field("moist_dyn", moist_dyn_array) + mr => mr_array%bundle + moist_dyn => moist_dyn_array%bundle + call moisture_fields%get_field("ls_mr", ls_mr_array) + call moisture_fields%get_field("ls_moist_dyn", ls_moist_dyn_array) + ls_mr => ls_mr_array%bundle + ls_moist_dyn => ls_moist_dyn_array%bundle + + ! Input + call ls_fields%get_field('ls_u', ls_u) + call ls_fields%get_field('ls_rho', ls_rho) + call ls_fields%get_field('ls_theta', ls_theta) + call ls_fields%get_field('ls_exner', ls_exner) + + call prognostic_fields%get_field('u', u) + call prognostic_fields%get_field('rho', rho) + call prognostic_fields%get_field('theta', theta) + call prognostic_fields%get_field('exner', exner) + + call r_u%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, W2) ) + call r_rho%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, W3) ) + call r_exner%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, W3) ) + call r_theta%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, Wtheta) ) + call clone_bundle(mr, r_mr, nummr) + + call p_u%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, W2) ) + call p_rho%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, W3) ) + call p_exner%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, W3) ) + call p_theta%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, Wtheta) ) + call clone_bundle(mr, p_mr, nummr) + + call n1_u%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, W2) ) + call n1_rho%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, W3) ) + call n1_exner%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, W3) ) + call n1_theta%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, Wtheta) ) + call clone_bundle(mr, n1_mr, nummr) + + call n2_u%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, W2) ) + call n2_rho%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, W3) ) + call n2_exner%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, W3) ) + call n2_theta%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, Wtheta) ) + call clone_bundle(mr, n2_mr, nummr) + + call diff_u%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, W2) ) + call diff_rho%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, W3) ) + call diff_exner%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, W3) ) + call diff_theta%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, Wtheta) ) + call clone_bundle(mr, diff_mr, nummr) + + call invoke( assign_field_random_kernel_type( r_u, 1.0_r_def ) , & + assign_field_random_kernel_type( r_rho, 1.0_r_def ) , & + assign_field_random_kernel_type( r_theta, 1.0_r_def ) , & + assign_field_random_kernel_type( r_exner, 1.0_r_def ) ) + + do i = 1, nummr + call invoke( assign_field_random_kernel_type( r_mr(i), 1.0_r_def ) ) + end do + + gamma_u = 1.e4_r_def + gamma_theta = 1.0_r_def + gamma_rho = 1.e0_r_def + gamma_exner = 6.5e-1_r_def + + if ( use_moisture ) then + gamma_mr = 1.e-3_r_def + else + gamma_mr = 0.0_r_def + end if + + gamma_u = gamma_u /24.0_r_def + gamma_theta = gamma_theta / 24.0_r_def + gamma_rho = gamma_rho / 24.0_r_def + gamma_exner = gamma_exner / 24.0_r_def + gamma_mr = gamma_mr / 24.0_r_def + + do n = 1, number_gamma_values + gamma_u = gamma_u / 2.0_r_def + gamma_theta = gamma_theta / 2.0_r_def + gamma_rho = gamma_rho / 2.0_r_def + gamma_exner = gamma_exner / 2.0_r_def + gamma_mr = gamma_mr / 2.0_r_def + + ! Initial conditions for nonlinear perturbed + call invoke( aX_plus_Y( n1_u, gamma_u, r_u, ls_u ), & + aX_plus_Y( n1_theta, gamma_theta, r_theta, ls_theta ), & + aX_plus_Y( n1_rho, gamma_rho, r_rho, ls_rho ), & + aX_plus_Y( n1_exner, gamma_exner, r_exner, ls_exner ) ) + + do i = 1, nummr + call invoke( aX_plus_Y( n1_mr(i), gamma_mr, & + r_mr(i), ls_mr(i) ) ) + end do + + ! Initial conditions for nonlinear unperturbed + call invoke( setval_X( n2_u, ls_u ), & + setval_X( n2_theta, ls_theta ), & + setval_X( n2_rho, ls_rho ), & + setval_X( n2_exner, ls_exner ) ) + + do i = 1, nummr + call invoke( setval_X( n2_mr(i), ls_mr(i) ) ) + end do + + ! Initial conditions for the linear + call invoke( a_times_X( p_u, gamma_u, r_u ), & + a_times_X( p_theta, gamma_theta, r_theta ), & + a_times_X( p_rho, gamma_rho, r_rho ), & + a_times_X( p_exner, gamma_exner, r_exner ) ) + + do i = 1, nummr + call invoke( a_times_X( p_mr(i), gamma_mr, & + r_mr(i) ) ) + end do + + ! Multiple timesteps + do t = parse_instance(timestep_start), parse_instance(timestep_end) + +!------------------------------------------------------------------------ +! Nonlinear model run with perturbed state +!------------------------------------------------------------------------ + call invoke( setval_X( u, n1_u ), & + setval_X( theta, n1_theta ), & + setval_X( rho, n1_rho ), & + setval_X( exner, n1_exner ) ) + do i = 1, nummr + call invoke( setval_X( mr(i), n1_mr(i) ) ) + enddo + call moist_dyn_factors_alg( moist_dyn, mr ) + + ! Clean up solvers and timestep method left over from running + ! the model to set up model state and linear state + call final_si_operators() + call semi_implicit_solver_alg_final() + call finalise_model(modeldb) + + call initialise_model( mesh, & + modeldb ) + + call gungho_step( mesh, & + twod_mesh, & + modeldb, & + model_clock ) + + call finalise_model(modeldb) + + call invoke( setval_X( n1_u, u ), & + setval_X( n1_theta, theta ), & + setval_X( n1_rho, rho ), & + setval_X( n1_exner, exner) ) + + do i = 1, nummr + call invoke( setval_X( n1_mr(i), mr(i) ) ) + end do + +!------------------------------------------------------------------------ +! Linear model run +!------------------------------------------------------------------------ + call invoke( setval_X( u, p_u ), & + setval_X( theta, p_theta ), & + setval_X( rho, p_rho ), & + setval_X( exner, p_exner) ) + do i = 1, nummr + call invoke( setval_X( mr(i), p_mr(i) ) ) + end do + call tl_moist_dyn_factors_alg( moist_dyn, mr ) + + ! Set the linearisation state + ! The linearisation state is copied from the nonlinear run, + ! every few timesteps. + if ( mod(t, update_ls_frequency) == 0 ) then + call invoke( setval_X( ls_u, n2_u ), & + setval_X( ls_theta, n2_theta ), & + setval_X( ls_rho, n2_rho ), & + setval_X( ls_exner, n2_exner) ) + do i = 1, nummr + call invoke( setval_X( ls_mr(i), n2_mr(i) ) ) + enddo + call moist_dyn_factors_alg( ls_moist_dyn, ls_mr ) + endif + + call finalise_linear_model() + call initialise_linear_model( mesh, & + modeldb ) + + call linear_step( mesh, & + twod_mesh, & + modeldb, & + model_clock ) + + call finalise_linear_model() + + + call invoke( setval_X( p_u, u ), & + setval_X( p_theta, theta ), & + setval_X( p_rho, rho ), & + setval_X( p_exner, exner) ) + do i = 1, nummr + call invoke( setval_X( p_mr(i), mr(i) ) ) + enddo + +!------------------------------------------------------------------------ +! Nonlinear model run with un-perturbed state +!------------------------------------------------------------------------ + call invoke( setval_X( u, n2_u ), & + setval_X( theta, n2_theta ), & + setval_X( rho, n2_rho ), & + setval_X( exner, n2_exner ) ) + + do i = 1, nummr + call invoke( setval_X( mr(i), n2_mr(i) ) ) + end do + + call moist_dyn_factors_alg( moist_dyn, mr ) + + + call initialise_model( mesh, & + modeldb ) + + call gungho_step( mesh, & + twod_mesh, & + modeldb, & + model_clock ) + + call finalise_model(modeldb) + + call invoke( setval_X( n2_u, u ), & + setval_X( n2_theta, theta ), & + setval_X( n2_rho, rho ), & + setval_X( n2_exner, exner ) ) + + do i = 1, nummr + call invoke( setval_X( n2_mr(i), mr(i) ) ) + call invoke( setval_X( ls_mr(i), mr(i) ) ) + end do + + end do ! timesteps loop + +!------------------------------------------------------------------------ +! Calculate norms +!------------------------------------------------------------------------ + call invoke( X_minus_Y( diff_u, n1_u, n2_u ), & + inc_X_minus_Y( diff_u, p_u ), & + X_innerproduct_X( norm_u, diff_u ) ) + call invoke( X_minus_Y( diff_rho, n1_rho, n2_rho ), & + inc_X_minus_Y( diff_rho, p_rho ), & + X_innerproduct_X( norm_rho, diff_rho ) ) + call invoke( X_minus_Y( diff_exner, n1_exner, n2_exner ), & + inc_X_minus_Y( diff_exner, p_exner ), & + X_innerproduct_X( norm_exner, diff_exner ) ) + call invoke( X_minus_Y( diff_theta, n1_theta, n2_theta ), & + inc_X_minus_Y( diff_theta, p_theta ), & + X_innerproduct_X( norm_theta, diff_theta ) ) + + norm_mr=0.0_r_def + do i = 1, nummr + call invoke( & + X_minus_Y( diff_mr(i), n1_mr(i), n2_mr(i) ), & + inc_X_minus_Y( diff_mr(i), p_mr(i) ), & + X_innerproduct_X( norm_mr_tmp, diff_mr(i) ) ) + norm_mr = norm_mr + norm_mr_tmp + end do + + norm_diff = norm_rho + norm_theta + norm_exner + norm_u + norm_mr + +!---------------------------------------------------------------------------- +! Print values out to enable plotting a graph +!---------------------------------------------------------------------------- + write( log_scratch_space, & + '( A, E32.12, A, E32.12 )' ) & + ' gamma_total = ' , gamma_u, & + ' norm = ' , norm_diff + + call log_event( log_scratch_space, LOG_LEVEL_INFO ) + + write( log_scratch_space, & + '(A, E32.12, A, E32.12)' ) & + ' gamma_rho = ' , gamma_u, & + ' norm = ' , norm_rho + + call log_event( log_scratch_space, LOG_LEVEL_INFO ) + + write( log_scratch_space, & + '(A, E32.12, A, E32.12)' ) & + ' gamma_exner = ' , gamma_u, & + ' norm = ' , norm_exner + + call log_event( log_scratch_space, LOG_LEVEL_INFO ) + + write( log_scratch_space, & + '(A, E32.12, A, E32.12)' ) & + ' gamma_theta = ' , gamma_u, & + ' norm = ' , norm_theta + + call log_event( log_scratch_space, LOG_LEVEL_INFO ) + + write( log_scratch_space, & + '(A, E32.12, A, E32.12)' ) & + ' gamma_u = ' , gamma_u, & + ' norm = ' , norm_u + + call log_event( log_scratch_space, LOG_LEVEL_INFO ) + + write( log_scratch_space, & + '(A, E32.12, A, E32.12)' ) & + ' gamma_mr = ' , gamma_u, & + ' norm = ' , norm_mr + + call log_event( log_scratch_space, LOG_LEVEL_INFO ) + +!---------------------------------------------------------------------------- +! Test that gives a PASS or FAIL +!---------------------------------------------------------------------------- + if (n == 2) then + + pass_str = "PASS" + + write( log_scratch_space, '(A)' ) & + "TL Test: " // trim(label) + call log_event( log_scratch_space, LOG_LEVEL_INFO ) + + write( log_scratch_space, '(A)' ) "Total " + call log_event( log_scratch_space, LOG_LEVEL_INFO ) + call convergence_pass_string( norm_diff, norm_diff_prev, pass_str, tol ) + + write( log_scratch_space, '(A)' ) "Theta " + call log_event( log_scratch_space, LOG_LEVEL_INFO ) + call convergence_pass_string( norm_theta, norm_theta_prev, pass_str, tol ) + + write( log_scratch_space, '(A)' ) "Rho " + call log_event( log_scratch_space, LOG_LEVEL_INFO ) + call convergence_pass_string( norm_rho, norm_rho_prev, pass_str, tol ) + + write( log_scratch_space, '(A)' ) "Exner " + call log_event( log_scratch_space, LOG_LEVEL_INFO ) + call convergence_pass_string( norm_exner, norm_exner_prev, pass_str, tol ) + + write( log_scratch_space, '(A)' ) "u " + call log_event( log_scratch_space, LOG_LEVEL_INFO ) + call convergence_pass_string( norm_u, norm_u_prev, pass_str, tol ) + + if ( use_moisture ) then + write( log_scratch_space, '(A)' ) "mr " + call log_event( log_scratch_space, LOG_LEVEL_INFO ) + call convergence_pass_string( norm_mr, norm_mr_prev, pass_str, tol ) + end if + + write( log_scratch_space,'(" test",A32," : ",A4)' ) trim(label), pass_str + call log_event( log_scratch_space, LOG_LEVEL_INFO ) + + end if + + norm_diff_prev = norm_diff + norm_rho_prev = norm_rho + norm_theta_prev = norm_theta + norm_exner_prev = norm_exner + norm_u_prev = norm_u + norm_mr_prev = norm_mr + + end do ! Loop over values of gamma + + end subroutine test_timesteps_random + + !> @brief Test the convergence rate. + !> @details If the convergence rate is not close to 4, within the + !! specified tolerance, then change the pass string to FAIL. + !! This is a special routine for this particular test, to allow + !! for norms that have already had the square root applied. + !> @param[in] norm_diff Current norm + !> @param[in] norm_diff_prev Previous norm + !> @param[inout] pass_str Pass string (either PASS or FAIL) + !> @param[in] tol Test Tolerance + subroutine convergence_pass_string( norm_diff, norm_diff_prev, pass_str, tol ) + implicit none + + real(r_def), intent(in) :: norm_diff, norm_diff_prev + character(len=4), intent(inout) :: pass_str + real(r_def), intent(in) :: tol + + real(r_def) :: conv_rate + conv_rate = norm_diff_prev / norm_diff + + conv_rate = norm_diff_prev / norm_diff + write( log_scratch_space, '(A, E16.8)') & + "Convergence rate: ", conv_rate + call log_event( log_scratch_space, LOG_LEVEL_INFO ) + + if ( abs( conv_rate - 4.0_r_def ) >= tol ) then + pass_str = "FAIL" + endif + + end subroutine convergence_pass_string + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Converts a string to a timestep number. + !> + !> @param[in] string Correctly formatted timestep number. + !> @return Timestep number. + !> + function parse_instance( string ) result(instance) + + implicit none + + character(*), intent(in) :: string + integer(i_timestep) :: instance + + integer :: string_size + integer :: format_size + integer :: status + character(:), allocatable :: fmt + + ! The "I0" formatting string does not work for input. Instead the exact + ! number of digits to be read must be requested. Thus we need to construct + ! a format string from the size of the string. + ! + string_size = len(string) + format_size = string_size + 3 + allocate( character(format_size) :: fmt, stat=status ) + if ( status /= 0 ) then + call log_event( & + 'TL test parse_instance: Unable to allocate format', & + log_level_error ) + end if + write( fmt, '("(I", I0, ")")' ) string_size + read( string, fmt ) instance + deallocate( fmt ) + + if ( instance < 0 ) then + call log_event( & + 'TL test parse_instance: Instances may not be negative', & + log_level_error ) + end if + + end function parse_instance + +end module tl_test_timesteps_random_alg_mod diff --git a/science/linear/integration-test/tl_test/tl_test_transport_control_mod.x90 b/science/linear/integration-test/tl_test/tl_test_transport_control_mod.x90 index 02d96d2e..d74fd3df 100644 --- a/science/linear/integration-test/tl_test/tl_test_transport_control_mod.x90 +++ b/science/linear/integration-test/tl_test/tl_test_transport_control_mod.x90 @@ -21,7 +21,6 @@ module tl_test_transport_control_mod use function_space_collection_mod, only: function_space_collection use derived_config_mod, only: bundle_size use driver_modeldb_mod, only: modeldb_type - use moist_dyn_mod, only: num_moist_factors use mr_indices_mod, only: nummr use log_mod, only: log_event, & log_scratch_space, & @@ -31,7 +30,7 @@ module tl_test_transport_control_mod only: gungho_transport_control_alg_init, & gungho_transport_control_alg use tl_transport_control_alg_mod, only: tl_transport_control_alg - use tl_test_convergence_rate_check, only: convergence_rate_check + use tl_test_convergence_rate_check, only: array_convergence_rate_check implicit none @@ -61,10 +60,11 @@ module tl_test_transport_control_mod type(field_type), pointer :: ls_u => null() type(field_type), pointer :: ls_rho => null() type(field_type), pointer :: ls_theta => null() - type(field_type), pointer :: ls_exner => null() - type(field_type), pointer :: ls_moist_dyn(:) => null() type(field_type), pointer :: ls_mr(:) => null() - type(field_type), pointer :: moist_dyn(:) => null() + type(field_type), pointer :: ls_exner => null() + type(field_type), pointer :: u => null() + type(field_type), pointer :: rho => null() + type(field_type), pointer :: theta => null() type(field_type), pointer :: mr(:) => null() type(field_type) :: state(bundle_size) @@ -75,10 +75,7 @@ module tl_test_transport_control_mod type(field_type) :: n1_rhs(bundle_size) type(field_type) :: n2_rhs(bundle_size) type(field_type) :: diff(bundle_size) - type(field_type), dimension(num_moist_factors) :: p_moist_dyn - type(field_type), dimension(num_moist_factors) :: n_moist_dyn - type(field_type), dimension(num_moist_factors) :: r_moist_dyn - type(field_type), dimension(num_moist_factors) :: diff_moist_dyn + type(field_type), dimension(nummr) :: p_mr_in type(field_type), dimension(nummr) :: n_mr_in type(field_type), dimension(nummr) :: ls_mr_out @@ -91,36 +88,32 @@ module tl_test_transport_control_mod type(field_collection_type), pointer :: moisture_fields => null() type(field_array_type), pointer :: mr_array => null() - type(field_array_type), pointer :: moist_dyn_array => null() type(field_array_type), pointer :: ls_mr_array => null() - type(field_array_type), pointer :: ls_moist_dyn_array => null() - - real(r_def) :: gamma_u, gamma_rho, gamma_exner, gamma_theta - real(r_def) :: gamma_moist_dyn, gamma_mr - real(r_def) :: norm_u, norm_rho, norm_exner, norm_theta - real(r_def) :: norm_moist_dyn, norm_mr - real(r_def) :: norm_diff, norm_diff_prev - real(r_def) :: norm_moist_dyn_tmp, norm_mr_tmp - real(r_def), parameter :: tol = 1.e-2_r_def + real(r_def) :: gamma_u, gamma_rho, gamma_theta + real(r_def) :: gamma_mr + real(r_def) :: norm_u, norm_rho, norm_theta + real(r_def) :: norm_mr + real(r_def) :: norm_mr_tmp + integer(i_def), parameter :: n_variables = 4 + real(r_def) :: array_norm(n_variables), array_norm_prev(n_variables) + character(str_def) :: array_names(n_variables) + real(r_def), parameter :: tol = 5.e-2_r_def + real(r_def), parameter :: indiv_tol = 5.e-1_r_def integer :: n, outer, i call log_event( "TL Test: " // trim(label), & - LOG_LEVEL_INFO ) + LOG_LEVEL_INFO ) prognostic_fields => modeldb%fields%get_field_collection( & "prognostic_fields") ls_fields => modeldb%fields%get_field_collection("ls_fields") moisture_fields => modeldb%fields%get_field_collection("moisture_fields") call moisture_fields%get_field("mr", mr_array) - call moisture_fields%get_field("moist_dyn", moist_dyn_array) mr => mr_array%bundle - moist_dyn => moist_dyn_array%bundle call moisture_fields%get_field("ls_mr", ls_mr_array) - call moisture_fields%get_field("ls_moist_dyn", ls_moist_dyn_array) ls_mr => ls_mr_array%bundle - ls_moist_dyn => ls_moist_dyn_array%bundle ! Input call ls_fields%get_field('ls_u', ls_u) @@ -128,6 +121,12 @@ module tl_test_transport_control_mod call ls_fields%get_field('ls_theta', ls_theta) call ls_fields%get_field('ls_exner', ls_exner) + ! Overwrite the moisture fields with potential temperature + do i = 1, nummr + call invoke( & + setval_X( ls_mr(i), ls_theta)) + end do + call ls_state(igh_u)%initialise( vector_space = ls_u%get_function_space() ) call ls_state(igh_t)%initialise( vector_space = ls_theta%get_function_space() ) call ls_state(igh_d)%initialise( vector_space = ls_rho%get_function_space() ) @@ -144,11 +143,6 @@ module tl_test_transport_control_mod call clone_bundle(ls_state, n2_rhs, bundle_size) call clone_bundle(ls_state, diff, bundle_size) - call clone_bundle(moist_dyn, p_moist_dyn, num_moist_factors) - call clone_bundle(moist_dyn, r_moist_dyn, num_moist_factors) - call clone_bundle(moist_dyn, n_moist_dyn, num_moist_factors) - call clone_bundle(moist_dyn, diff_moist_dyn, num_moist_factors) - call clone_bundle(mr, r_mr, nummr) call clone_bundle(mr, p_mr_in, nummr) call clone_bundle(mr, n_mr_in, nummr) @@ -161,23 +155,36 @@ module tl_test_transport_control_mod setval_X(ls_state(igh_u), ls_u ), & setval_X(ls_state(igh_t), ls_theta), & setval_X(ls_state(igh_d), ls_rho ), & - setval_X(ls_state(igh_p), ls_exner), & + setval_X(ls_state(igh_p), ls_exner ), & setval_X(ls_advected_u, ls_u ) ) call gungho_transport_control_alg_init( mesh ) - call invoke( assign_field_random_kernel_type( random(igh_u), 1.0_r_def ) , & - assign_field_random_kernel_type( random(igh_d), 1.0_r_def ) , & - assign_field_random_kernel_type( random(igh_t), 1.0_r_def ) , & - assign_field_random_kernel_type( random(igh_p), 1.0_r_def ) ) + ! Perturbation + call prognostic_fields%get_field('u', u) + call prognostic_fields%get_field('rho', rho) + call prognostic_fields%get_field('theta', theta) + + ! Set the 'random perturbation as the ls_state plus pert' + call invoke( name = "copy_fields_to_state", & + setval_X(random(igh_u), ls_u ), & + setval_X(random(igh_t), ls_theta), & + setval_X(random(igh_d), ls_rho ) ) + + do i = 1, nummr + call invoke( & + setval_X( r_mr(i), ls_mr(i) )) + end do - do i = 1, num_moist_factors - call invoke( assign_field_random_kernel_type( r_moist_dyn(i), 1.0_r_def ) ) - enddo + call invoke( name = "inc", & + inc_X_plus_Y(random(igh_u), u ), & + inc_X_plus_Y(random(igh_t), theta), & + inc_X_plus_Y(random(igh_d), rho ) ) do i = 1, nummr - call invoke( assign_field_random_kernel_type( r_mr(i), 1.0_r_def ) ) - enddo + call invoke( & + inc_X_plus_Y( r_mr(i), theta )) + end do call set_bundle_scalar( 0.0_r_def, n1_rhs, bundle_size ) @@ -192,19 +199,15 @@ module tl_test_transport_control_mod outer, & cheap_update = .false.) - gamma_u = 7.e8_r_def - gamma_theta = 1.e5_r_def - gamma_rho = 1.e5_r_def - gamma_exner = 1.e5_r_def - gamma_moist_dyn=1.e0_r_def - gamma_mr = 1.e0_r_def + gamma_u = 1._r_def + gamma_theta = 1._r_def + gamma_rho = 1._r_def + gamma_mr = 1._r_def do n=1,2 gamma_u = gamma_u/2.0_r_def gamma_theta = gamma_theta/2.0_r_def gamma_rho = gamma_rho/2.0_r_def - gamma_exner = gamma_exner/2.0_r_def - gamma_moist_dyn = gamma_moist_dyn/2.0_r_def gamma_mr = gamma_mr/2.0_r_def call set_bundle_scalar( 0.0_r_def, n2_rhs, bundle_size ) @@ -212,7 +215,7 @@ module tl_test_transport_control_mod call invoke( a_times_X( p_state(igh_u), gamma_u, random(igh_u) ), & a_times_X( p_state(igh_t), gamma_theta, random(igh_t) ), & a_times_X( p_state(igh_d), gamma_rho, random(igh_d) ), & - a_times_X( p_state(igh_p), gamma_exner, random(igh_p) ), & + setval_c( p_state(igh_p), 0.0_r_def ), & a_times_X( p_advected_u, gamma_u, random(igh_u) ) ) ! state = ls_state + p_state @@ -221,13 +224,6 @@ module tl_test_transport_control_mod call invoke( setval_X( n_advected_u, ls_advected_u), & inc_X_plus_Y(n_advected_u, p_advected_u ) ) - do i = 1, num_moist_factors - call invoke( & - a_times_X( p_moist_dyn(i), gamma_moist_dyn, r_moist_dyn(i) ), & - setval_X( n_moist_dyn(i), ls_moist_dyn(i) ), & - inc_X_plus_Y( n_moist_dyn(i), p_moist_dyn(i) ) ) - end do - do i = 1, nummr call invoke( & a_times_X( p_mr_in(i), gamma_mr, r_mr(i) ), & @@ -263,63 +259,43 @@ module tl_test_transport_control_mod call invoke( & inc_X_minus_Y( diff(igh_u), p_rhs(igh_u) ), & inc_X_minus_Y( diff(igh_t), p_rhs(igh_t) ), & - inc_X_minus_Y( diff(igh_d), p_rhs(igh_d) ), & - inc_X_minus_Y( diff(igh_p), p_rhs(igh_p) ) ) + inc_X_minus_Y( diff(igh_d), p_rhs(igh_d) ) ) - call invoke( X_innerproduct_X( norm_u, diff(igh_u) ) , & + call invoke( & + X_innerproduct_X( norm_u, diff(igh_u) ) , & X_innerproduct_X( norm_rho, diff(igh_d) ) , & - X_innerproduct_X( norm_theta, diff(igh_t) ) , & - X_innerproduct_X( norm_exner, diff(igh_p) ) ) + X_innerproduct_X( norm_theta, diff(igh_t) ) ) - norm_moist_dyn=0.0_r_def - do i = 1, num_moist_factors + + do i = 1, nummr call invoke( & - X_minus_Y( diff_moist_dyn(i), n_moist_dyn(i), moist_dyn(i) ), & - inc_X_minus_Y( diff_moist_dyn(i), p_moist_dyn(i) ), & - X_innerproduct_X( norm_moist_dyn_tmp, diff_moist_dyn(i) ) ) - norm_moist_dyn = norm_moist_dyn + norm_moist_dyn_tmp + X_minus_Y( diff_mr(i), n_mr_out(i), ls_mr_out(i) ), & + inc_X_minus_Y( diff_mr(i), p_mr_out(i) ) ) end do norm_mr=0.0_r_def do i = 1, nummr - call invoke( & - X_minus_Y( diff_mr(i), n_mr_out(i), ls_mr_out(i) ), & - inc_X_minus_Y( diff_mr(i), p_mr_out(i) ), & + call invoke( & X_innerproduct_X( norm_mr_tmp, diff_mr(i) ) ) norm_mr = norm_mr + norm_mr_tmp end do - write( log_scratch_space, & - '(A, E32.12, A, E32.12 )' ) & - ' norm_u = ' , norm_u, & - ' norm_rho = ' , norm_rho - - call log_event( log_scratch_space, LOG_LEVEL_INFO ) - - write( log_scratch_space, & - '(A, E32.12, A, E32.12 )' ) & - ' norm_theta = ' , norm_theta, & - ' norm_moist_dyn = ' , norm_moist_dyn - - call log_event( log_scratch_space, LOG_LEVEL_INFO ) - - write( log_scratch_space, & - '(A, E32.12, A, E32.12 )' ) & - ' norm_exner = ' , norm_exner, & - ' norm_mr = ' , norm_mr - - call log_event( log_scratch_space, LOG_LEVEL_INFO ) - - norm_diff = norm_rho + norm_theta + norm_exner + norm_u & - + norm_moist_dyn + norm_mr - norm_diff = sqrt(norm_diff) + array_norm(1) = norm_u + array_norm(2) = norm_rho + array_norm(3) = norm_theta + array_norm(4) = norm_mr + array_names(1) = 'norm_u' + array_names(2) = 'norm_rho' + array_names(3) = 'norm_theta' + array_names(4) = 'norm_mr' if (n == 2) then - call convergence_rate_check( norm_diff, norm_diff_prev, & - label, tol=tol ) + call array_convergence_rate_check( array_norm, array_norm_prev, & + array_names, n_variables, label, tol=tol, indiv_tol=indiv_tol) end if - norm_diff_prev = norm_diff + array_norm_prev = array_norm + end do end subroutine test_transport_control diff --git a/science/linear/integration-test/tl_test/tl_test_vorticity_mod.x90 b/science/linear/integration-test/tl_test/tl_test_vorticity_mod.x90 index 47e5cdfc..c129dd3f 100644 --- a/science/linear/integration-test/tl_test/tl_test_vorticity_mod.x90 +++ b/science/linear/integration-test/tl_test/tl_test_vorticity_mod.x90 @@ -203,8 +203,6 @@ module tl_test_vorticity_mod inc_X_minus_Y( diff, p_rhs_u ), & X_innerproduct_X( norm_diff, diff ) ) - norm_diff = sqrt(norm_diff) - write( log_scratch_space, '(A, E32.12, A, E32.12)' ) & 'gamma = ', gamma, ' norm = ', norm_diff diff --git a/applications/linear_model/plot_convergence/plot_convergence.py b/science/linear/plot_convergence/plot_convergence.py similarity index 54% rename from applications/linear_model/plot_convergence/plot_convergence.py rename to science/linear/plot_convergence/plot_convergence.py index 811421c8..6d7c42ea 100644 --- a/applications/linear_model/plot_convergence/plot_convergence.py +++ b/science/linear/plot_convergence/plot_convergence.py @@ -14,6 +14,7 @@ import os import pandas as pd import matplotlib.pyplot as plt +import numpy as np def plot_data(filename, axes, variable, color, shape): @@ -37,12 +38,55 @@ def plot_data(filename, axes, variable, color, shape): norm_df = pd.DataFrame(norm_line, columns=['gamma', 'norm']) datafile.close() + + # Square root (as the read data is only the inner product and + # has not included the square root to give the norm) + norm_df['norm'] = np.sqrt(norm_df['norm']) + + # Normalise - to give a relative error + normalise = norm_df['norm'].iloc[4] + norm_df['norm'] = norm_df['norm'] / normalise + + # Plot the data + if (CONFIG == 'nwp_gal9'): + ymin = 10**-2 + ymax = 10**2 + xmin = 10**-3 + xmax = 10**1 + elif (CONFIG == 'semi_implicit'): + ymin = 10**-2 + ymax = 10**2 + xmin = 10**-1 + xmax = 10**4 + elif (CONFIG == 'runge_kutta'): + ymin = 10**-3 + ymax = 10**2 + xmin = 10**-1 + xmax = 10**4 + else: + print(CONFIG+' not listed') + + norm_df.plot.scatter(x='gamma', y='norm', loglog=True, + xlim=(xmin, xmax), + ylim=(ymin, ymax), + ax=axes, color=color, + marker=shape, label = variable) + + # Check extremes + norm_min = norm_df['norm'].min() + norm_max = norm_df['norm'].max() + if (norm_min < ymin) : + print('Warning: Min value is'+ str(norm_min)) + if (norm_max > ymax) : + print('Warning: Max value is'+ str(norm_max)) + + # Plot the expected line + centre = norm_df['gamma'].iloc[4] + expected_x = [10**-2 *centre, centre, 100* centre] + expected_y = [10**-2, 10**0, 100] + plt.plot(expected_x, expected_y) - norm_df.plot.scatter(x='gamma', y='norm', loglog=True, xlim=(10**0, 10**5), - ylim=(10**-5, 10**0), ax=axes, color=color, - marker=shape) - - + def make_plot(directory, filename): ''' Plot the data for the different prognostic variables, together with the @@ -55,27 +99,25 @@ def make_plot(directory, filename): plot_data(directory + filename, axs, 'gamma_u', 'r', 'o') plot_data(directory + filename, axs, 'gamma_exner', 'b', 's') plot_data(directory + filename, axs, 'gamma_theta', 'g', '^') - plot_data(directory + filename, axs, 'gamma_total', 'black', 'x') - plot_data(directory + filename, axs, 'gamma_mr', 'm', '*') + plot_data(directory + filename, axs, 'gamma_mr', 'black', 'x') - expected_x = [10**0, 10**5] - expected_y = [10**-5, 10**0] - plt.plot(expected_x, expected_y) - - plt.legend(['Expected gradient (linear)', 'density', 'momentum', - 'exner pressure', 'potential temperature', 'total', - 'moisture mixing ratios'], - loc='lower right') + plt.legend(loc='lower right') plt.xlabel('Gamma') plt.ylabel('Relative error') plt.title('Validity of the tangent linear model') - plt.show() + # To show the plot to the screen, uncommment plt.show() + #plt.show() + + # Save the plot to a file + plt.savefig(directory + filename + "convergence_plot.png") if __name__ == "__main__": DATA_DIRECTORY = os.getcwd()+'/' DATA_FILENAME = 'outfile' + CONFIG = os.getenv('CONFIG') + print(CONFIG) make_plot(DATA_DIRECTORY, DATA_FILENAME) diff --git a/science/linear/plot_convergence/plot_convergence.sh b/science/linear/plot_convergence/plot_convergence.sh new file mode 100755 index 00000000..8b60e533 --- /dev/null +++ b/science/linear/plot_convergence/plot_convergence.sh @@ -0,0 +1,118 @@ +############################################################################## +# (c) Crown copyright 2022 Met Office. All rights reserved. +# The file LICENCE, distributed with this code, contains details of the terms +# under which the code may be used. +############################################################################## + +#---------------------------------------------------------------------- +# Plots the convergence rate of the tangent linear model. +#---------------------------------------------------------------------- + +# INSTRUCTIONS TO RUN +# 1. Specify the config_list +# 2. Run using . plot_convergence.sh, from the plot_convergence directory + +# SCIENCE DETAILS +# The relative linearisation error is +# E = || N(x+ gamma x') - N(x) - L(x) gamma x' || / || L(x) gamma x' || +# where N=nonlinear model, L=linear model, x=linearisation state +# x'=perturbation, gamma=scalar. +# From the Taylor series expansion, E(gamma) = O(gamma) i.e. of the order gamma +# So the relative error should be a linear function of gamma + +# SCRIPT STEPS +# 1. Build the model +# 2. Produce the data: The integration test tl_test_timesteps is extended by +# running over 10 values of gamma, rather than 2 values of gamma. +# 3. Plot the data: The data is plotted for each prognostic variable. + +# EXTENSION +# The plot_configuration.nml can also be extended to other configurations e.g +# * increase the number of timesteps (timesteps_end) +# * increase the number of timesteps between updating the linearisation state +# (update_ls_frequency) + +#-------------------------------------------------------------------------- + +######################## Functions ########################################## +build(){ + echo $CONFIG " Building the executable" + + # Integration tests executable name + exe=$Linear_dir/test/$CONFIG/$CONFIG + + # Build the integration tests, unless that has already been completed + if [ -f $exe ] ; then + echo "Do not need to build the executable as $exe exists" + else + echo "$exe does not exist, so now building the executable" + cd $Root_dir/build + python3 local_build.py linear -t integration-tests + + if [$? -ne 0 ]; then + echo "Error building the executable" + return + fi + fi +} + +integration_test(){ + # Setup the configuration - to test with 10 values of gamma + echo $CONFIG " Setting up the configuration" + cd $Linear_dir/integration-test/$CONFIG/resources/ + cp ${CONFIG}_configuration.nml plot_configuration.nml + sed -i 's/number_gamma_values=2/number_gamma_values=10/g' plot_configuration.nml + if [ $? -ne 0 ]; then + echo "Error in creating plot_configuration.nml" + return + else + echo "plot_configuration.nml created in " $PWD + fi + + # Run the tl_test_timesteps integration test + echo $CONFIG " Running the integration test" + cd $Linear_dir/integration-test/$CONFIG + echo $PWD + $exe resources/plot_configuration.nml test_timesteps > outfile + if [ $? -ne 0 ]; then + echo "Error in creating outfile data" + return + else + echo "Data created successfully" + fi +} + +plot(){ + # Plot the data, together with the expected gradient + echo $CONFIG " Plotting the data" + cd $Linear_dir/integration-test/$CONFIG + python $Linear_dir/plot_convergence/plot_convergence.py +} + +#################### MAIN PROGRAM ################################################# + +# Modules +module purge +module use /home/users/lfricadmin/lmod +module load lfric/vn3.0 +module load scitools + +config_list=(nwp_gal9 semi_implicit runge_kutta) + +# Directory of this script +SCRIPT_DIR="$(dirname "$(realpath "$BASH_SOURCE")")" +# And parent directories +export Linear_dir="$(dirname $SCRIPT_DIR)" +export Parent_dir="$(dirname $Linear_dir)" +export Root_dir="$(dirname $Parent_dir)" +echo "Linear_dir" $Linear_dir +echo "Root_dir" $Root_dir + +for configuration in "${config_list[@]}"; do + echo $configuration + export CONFIG="$configuration" + + build + integration_test + plot +done