diff --git a/applications/gungho_model/example/configuration.nml b/applications/gungho_model/example/configuration.nml index ea60b8aa..8b500d9b 100644 --- a/applications/gungho_model/example/configuration.nml +++ b/applications/gungho_model/example/configuration.nml @@ -77,11 +77,12 @@ use_wavedynamics=.true., vector_invariant=.false., / &helmholtz_solver +fail_on_non_converged=.false., gcrk=8, method='prec_only', -monitor_convergence=.false., +monitor_convergence=.true., normalise=.true., -preconditioner='multigrid', +preconditioner='tridiagonal', si_pressure_a_tol=1.0e-8, si_pressure_maximum_iterations=400, si_pressure_tolerance=1.0e-4, @@ -187,8 +188,8 @@ monitor_convergence=.true., normalise=.true., reference_reset_time=3600.0, si_maximum_iterations=10, -si_method='block_gcr', -si_preconditioner='pressure', +si_method='prec_only', +si_preconditioner='multigrid', si_tolerance=1.0e-1, split_w=.true., / @@ -204,7 +205,7 @@ multigrid_chain_nitems=4, n_coarsesmooth=4, n_postsmooth=2, n_presmooth=2, -smooth_relaxation=0.8, +smooth_relaxation=0.5, / &esm_couple l_esm_couple_test=.false., diff --git a/build/local_build.py b/build/local_build.py index 2c9d445e..c79b9e09 100755 --- a/build/local_build.py +++ b/build/local_build.py @@ -82,14 +82,11 @@ def determine_project_path(project, root_dir): ) -def clone_dependency(values, temp_dep): +def clone_dependency(source, ref, temp_dep): """ Clone the physics dependencies into a temporary directory """ - source = values["source"] - ref = values["ref"] - commands = ( f"git -C {temp_dep} init", f"git -C {temp_dep} remote add origin {source}", diff --git a/dependencies.yaml b/dependencies.yaml index caab3648..9ad110cf 100644 --- a/dependencies.yaml +++ b/dependencies.yaml @@ -30,8 +30,8 @@ lfric_apps: ref: lfric_core: - source: git@github.com:MetOffice/lfric_core.git - ref: 2025.12.1 + source: git@github.com:thomasmelvin/lfric_core.git + ref: mixed_multigrid moci: source: git@github.com:MetOffice/moci.git diff --git a/rose-stem/app/gungho_model/opt/rose-app-mixed-mg.conf b/rose-stem/app/gungho_model/opt/rose-app-mixed-mg.conf new file mode 100644 index 00000000..4c74bb63 --- /dev/null +++ b/rose-stem/app/gungho_model/opt/rose-app-mixed-mg.conf @@ -0,0 +1,12 @@ +[namelist:helmholtz_solver] +method='bicgstab' +preconditioner='tridiagonal', +si_pressure_tolerance=1.0e-2 + +[namelist:mixed_solver] +si_method='prec_only' +si_preconditioner='multigrid', + +[namelist:multigrid] +n_coarsesmooth=4 +smooth_relaxation=0.5 diff --git a/rose-stem/site/common/gungho_model/tasks_gungho_model.cylc b/rose-stem/site/common/gungho_model/tasks_gungho_model.cylc index 17c6843c..08405835 100644 --- a/rose-stem/site/common/gungho_model/tasks_gungho_model.cylc +++ b/rose-stem/site/common/gungho_model/tasks_gungho_model.cylc @@ -33,6 +33,39 @@ "tsteps": 240, "plot_str": "baroclinic.py $NODAL_DATA_DIR/lfric_diag.nc $PLOT_DIR dcmip dry surface_pressure_temperature", }) %} +{% elif task_ns.conf_name == "baroclinic-mixed-mg-C24_MG" %} + + {% do task_dict.update({ + "opt_confs": ["baroclinic", "mixed-mg"], + "resolution": "C24_MG", + "DT": 3600, + "threads": 4, + "mpi_parts": 6, + "tsteps": 240, + "plot_str": "baroclinic.py $NODAL_DATA_DIR/lfric_diag.nc $PLOT_DIR dcmip dry surface_pressure_temperature", + }) %} +{% elif task_ns.conf_name == "baroclinic-C48_MG" %} + + {% do task_dict.update({ + "opt_confs": ["baroclinic"], + "resolution": "C48_MG", + "DT": 3600, + "threads": 4, + "mpi_parts": 6, + "tsteps": 240, + "plot_str": "baroclinic.py $NODAL_DATA_DIR/lfric_diag.nc $PLOT_DIR dcmip dry surface_pressure_temperature", + }) %} +{% elif task_ns.conf_name == "baroclinic-mixed-mg-C48_MG" %} + + {% do task_dict.update({ + "opt_confs": ["baroclinic", "mixed-mg"], + "resolution": "C48_MG", + "DT": 3600, + "threads": 4, + "mpi_parts": 6, + "tsteps": 240, + "plot_str": "baroclinic.py $NODAL_DATA_DIR/lfric_diag.nc $PLOT_DIR dcmip dry surface_pressure_temperature", + }) %} {% elif task_ns.conf_name == "baroclinic-pert-C24_MG" %} diff --git a/rose-stem/site/meto/groups/groups_gungho_model.cylc b/rose-stem/site/meto/groups/groups_gungho_model.cylc index 0605786b..fb36e7f0 100644 --- a/rose-stem/site/meto/groups/groups_gungho_model.cylc +++ b/rose-stem/site/meto/groups/groups_gungho_model.cylc @@ -67,6 +67,14 @@ {# EX1A groups #} {% do site_groups.update({ + "gungho_model_ex1a_test": [ + "gungho_model_baroclinic-C24_MG_ex1a_gnu_fast-debug-64bit", + "gungho_model_baroclinic-mixed-mg-C24_MG_ex1a_gnu_fast-debug-64bit", + "gungho_model_baroclinic-C24_MG_ex1a_cce_fast-debug-64bit", + "gungho_model_baroclinic-mixed-mg-C24_MG_ex1a_cce_fast-debug-64bit", + "gungho_model_baroclinic-C48_MG_ex1a_gnu_fast-debug-64bit", + "gungho_model_baroclinic-mixed-mg-C48_MG_ex1a_gnu_fast-debug-64bit", + ], "gungho_model_ex1a_developer": [ "gungho_model_agnesi_hyd_cart-BiP120x8-2000x2000_ex1a_gnu_fast-debug-64bit", "gungho_model_baroclinic-C24_MG_ex1a_gnu_fast-debug-64bit", diff --git a/science/gungho/rose-meta/lfric-gungho/HEAD/rose-meta.conf b/science/gungho/rose-meta/lfric-gungho/HEAD/rose-meta.conf index 15f3d968..25f6972c 100644 --- a/science/gungho/rose-meta/lfric-gungho/HEAD/rose-meta.conf +++ b/science/gungho/rose-meta/lfric-gungho/HEAD/rose-meta.conf @@ -3493,8 +3493,8 @@ help=?????? =?????? !kind=default sort-key=Panel-A03 -value-titles=None, Diagonal, Pressure -values='none', 'diagonal', 'pressure' +value-titles=None, Diagonal, Pressure, Multigrid +values='none', 'diagonal', 'pressure', 'multigrid' [namelist:mixed_solver=si_tolerance] compulsory=true diff --git a/science/gungho/source/algorithm/core_dynamics/intermesh_mappings_alg_mod.x90 b/science/gungho/source/algorithm/core_dynamics/intermesh_mappings_alg_mod.x90 index 19f953a4..868d2fce 100644 --- a/science/gungho/source/algorithm/core_dynamics/intermesh_mappings_alg_mod.x90 +++ b/science/gungho/source/algorithm/core_dynamics/intermesh_mappings_alg_mod.x90 @@ -55,6 +55,8 @@ module intermesh_mappings_alg_mod private :: prolong_rho_reversible ! Conservative and reversible prolong for W3 fields public :: map_rho_intermesh ! Public routine for conservative W3 fields public :: map_w2_intermesh ! Routine for mapping any W2 field + public :: map_w2h_intermesh ! Routine for mapping any W2h field + public :: map_w2h_intermesh_rsolver ! Routine for mapping any W2h rsolver field public :: map_w0_intermesh ! Routine for mapping any W0 field private :: restrict_mr ! Conservative restriction for mixing ratio private :: inject_mr ! Conservative injection for mixing ratio @@ -566,6 +568,105 @@ contains end subroutine map_w2_intermesh + !> @details An algorithm for mapping W2h vector fields + !> @param[in,out] target_field Target W2h field + !> @param[in] source_field Source W2h field + subroutine map_w2h_intermesh( target_field, source_field ) + + use sci_prolong_w2h_kernel_mod, only: prolong_w2h_kernel_type + use sci_restrict_w2h_kernel_mod, only: restrict_w2h_kernel_type + use sci_mapping_constants_mod, only: get_intermesh_weights_w2h + + implicit none + + type(field_type), intent(in) :: source_field + type(field_type), intent(inout) :: target_field + + type(field_type), pointer :: weights + type(mesh_type), pointer :: source_mesh + type(mesh_type), pointer :: target_mesh + type(integer_field_type), pointer :: face_selector_ew + type(integer_field_type), pointer :: face_selector_ns + integer(kind=i_def) :: source_ncells + integer(kind=i_def) :: target_ncells + + ! Get source and target meshes + source_mesh => source_field%get_mesh() + target_mesh => target_field%get_mesh() + + ! Figure out which of the meshes has fewer cells + source_ncells = source_mesh%get_ncells() + target_ncells = target_mesh%get_ncells() + + ! Map W2 vector field from source to target + if ( source_ncells > target_ncells ) then + face_selector_ew => get_face_selector_ew(target_mesh%get_id()) + face_selector_ns => get_face_selector_ns(target_mesh%get_id()) + call invoke( restrict_w2h_kernel_type(target_field, source_field, & + face_selector_ew, face_selector_ns) ) + + else if ( source_ncells < target_ncells ) then + weights => get_intermesh_weights_w2h(target_mesh, source_mesh) + ! Prolong kernel increments target_field, so need to set it to be 0 first + call invoke( setval_c(target_field, 0.0_r_def), & + prolong_w2h_kernel_type(target_field, source_field, weights) ) + + else + call invoke( setval_X(target_field, source_field) ) + end if + + end subroutine map_w2h_intermesh + + !> @details An algorithm for mapping W2h rsolver fields + !> @param[in,out] target_field Target W2h rsolver field + !> @param[in] source_field Source W2h rsolver field + subroutine map_w2h_intermesh_rsolver( target_field_rsolver, source_field_rsolver) + + use sci_prolong_w2h_kernel_mod, only: prolong_w2h_kernel_type + use sci_restrict_w2h_kernel_mod, only: restrict_w2h_kernel_type + use sci_mapping_constants_mod, only: get_intermesh_weights_w2h + use r_solver_field_mod, only: r_solver_field_type + + implicit none + + type(r_solver_field_type), intent(in) :: source_field_rsolver + type(r_solver_field_type), intent(inout) :: target_field_rsolver + + type(field_type), pointer :: weights + type(mesh_type), pointer :: source_mesh + type(mesh_type), pointer :: target_mesh + type(integer_field_type), pointer :: face_selector_ew + type(integer_field_type), pointer :: face_selector_ns + integer(kind=i_def) :: source_ncells + integer(kind=i_def) :: target_ncells + + ! Get source and target meshes + source_mesh => source_field_rsolver%get_mesh() + target_mesh => target_field_rsolver%get_mesh() + + ! Figure out which of the meshes has fewer cells + source_ncells = source_mesh%get_ncells() + target_ncells = target_mesh%get_ncells() + + ! Map W2 vector field from source to target + if ( source_ncells > target_ncells ) then + face_selector_ew => get_face_selector_ew(target_mesh%get_id()) + face_selector_ns => get_face_selector_ns(target_mesh%get_id()) + call invoke( restrict_w2h_kernel_type(target_field_rsolver, source_field_rsolver, & + face_selector_ew, face_selector_ns) ) + + else if ( source_ncells < target_ncells ) then + weights => get_intermesh_weights_w2h(target_mesh, source_mesh) + ! Prolong kernel increments target_field, so need to set it to be 0 first + call invoke( setval_c(target_field_rsolver, 0.0_r_def), & + prolong_w2h_kernel_type(target_field_rsolver, source_field_rsolver, weights) ) + + else + call invoke( setval_X(target_field_rsolver, source_field_rsolver) ) + end if + + end subroutine map_w2h_intermesh_rsolver + !> @details An algorithm for mapping W0 scalar fields !> @param[in,out] target_field Target field diff --git a/science/gungho/source/algorithm/solver/create_pressure_preconditioner_mod.F90 b/science/gungho/source/algorithm/solver/create_pressure_preconditioner_mod.F90 new file mode 100644 index 00000000..434afef6 --- /dev/null +++ b/science/gungho/source/algorithm/solver/create_pressure_preconditioner_mod.F90 @@ -0,0 +1,95 @@ +module create_pressure_preconditioner_mod + + !=============================================================================! + !> @file create_pressure_preconditioner_mod.F90 + !> @brief Create iterative preconditioner for (Helmholtz) pressure problem + !> @details This module provides a subroutine to create an iterative preconditioner for the + !! (Helmholtz) pressure problem, which is used in the semi-implicit preconditioner. + !=============================================================================! + + use log_mod, only: log_event, LOG_LEVEL_INFO, LOG_LEVEL_ERROR + use constants_mod, only: i_def + use pressure_operator_alg_mod, only: pressure_operator_type + use sci_preconditioner_mod, only: abstract_preconditioner_type + use pressure_precon_alg_mod, only: pressure_preconditioner_type + use pressure_diag_precon_alg_mod, only: pressure_diag_preconditioner_type + use multigrid_preconditioner_alg_mod, only: multigrid_preconditioner_type + use sci_null_preconditioner_alg_mod, only: null_preconditioner_type + + implicit none + + private + public :: create_pressure_preconditioner + + contains + +!=============================================================================! + !> @brief Create operator and preconditioner for (Helmholtz) pressure problem + !> @details Called by init method of this module, but also by adj_semi_implicit_solver_alg_mod, + !! adjt_mixed_schur_preconditioner_alg_mod and adjt_mixed_solver_alg_mod + !> @param[in] fs_id ID of the pressure function space + !> @param[out] pressure_operator_out Output (Helmholtz) pressure operator + !> @param[out] pressure_preconditioner_out Output (Helmholtz) pressure preconditioner + !> @param[in] level Multigrid level to create preconditioner on + subroutine create_pressure_preconditioner( fs_id, pressure_operator_out, pressure_preconditioner_out, level ) + + use helmholtz_solver_config_mod, only: helmholtz_preconditioner => preconditioner, & + preconditioner_none, & + preconditioner_diagonal, & + preconditioner_tridiagonal, & + preconditioner_multigrid + use function_space_mod, only: function_space_type + + implicit none + + integer(kind=i_def), intent(in) :: fs_id + integer(kind=i_def), intent(in) :: level + + ! Output operator and preconditioner for (Helmholtz) pressure problem + type(pressure_operator_type), intent(out) :: pressure_operator_out + class(abstract_preconditioner_type), allocatable, intent(out) :: pressure_preconditioner_out + + ! Vertical pressure preconditioner + type(pressure_preconditioner_type) :: Hz_preconditioner + + integer(kind=i_def) :: n_fields + + pressure_operator_out = pressure_operator_type(level) + + call log_event( "create_pressure_preconditioner: starting", LOG_LEVEL_INFO ) + + ! Allocate pressure preconditioner of correct type + select case(helmholtz_preconditioner) + case(PRECONDITIONER_NONE) + call log_event( "No pressure preconditioner specified", LOG_LEVEL_INFO ) + allocate( pressure_preconditioner_out, & + source = null_preconditioner_type() ) + case(PRECONDITIONER_DIAGONAL) + call log_event( "Creating diagonal pressure preconditioner", LOG_LEVEL_INFO ) + allocate( pressure_preconditioner_out, & + source = pressure_diag_preconditioner_type() ) + case(PRECONDITIONER_TRIDIAGONAL) + call log_event( "Creating tridiagonal pressure preconditioner", LOG_LEVEL_INFO ) + allocate( pressure_preconditioner_out, & + source = pressure_preconditioner_type(level) ) + case(PRECONDITIONER_MULTIGRID) + call log_event( "Creating multigrid pressure preconditioner", LOG_LEVEL_INFO ) + if (level /= 1_i_def) & + call log_event( "Pressure Multigrid has to be constructed on level 1", LOG_LEVEL_ERROR) + + Hz_preconditioner = pressure_preconditioner_type(level=1_i_def) + n_fields = 1_i_def + allocate( pressure_preconditioner_out, & + source = multigrid_preconditioner_type( (/fs_id/), & + n_fields, & + pressure_operator_out, & + Hz_preconditioner ) ) + case default + call log_event( "Unknown pressure preconditioner specified", LOG_LEVEL_ERROR) + end select + + call log_event( "create_pressure_preconditioner: done", LOG_LEVEL_INFO ) + + end subroutine create_pressure_preconditioner + + end module create_pressure_preconditioner_mod diff --git a/science/gungho/source/algorithm/solver/create_pressure_solver_mod.F90 b/science/gungho/source/algorithm/solver/create_pressure_solver_mod.F90 new file mode 100644 index 00000000..d70b6d68 --- /dev/null +++ b/science/gungho/source/algorithm/solver/create_pressure_solver_mod.F90 @@ -0,0 +1,138 @@ +module create_pressure_solver_mod + + !=============================================================================! + !> @file create_pressure_solver_mod.F90 + !> @brief Create iterative solver for (Helmholtz) pressure problem + !> @details This module provides a subroutine to create an iterative solver for the + !! (Helmholtz) pressure problem, which is used in the semi-implicit solver. + !=============================================================================! + + use sci_iterative_solver_mod, only: abstract_iterative_solver_type, & + conjugate_gradient_type, bicgstab_type, & + gmres_type, fgmres_type, gcr_type, & + precondition_only_type, jacobi_type + + use sci_preconditioner_mod, only: abstract_preconditioner_type + use pressure_operator_alg_mod, only: pressure_operator_type + use log_mod, only: log_event, LOG_LEVEL_INFO, LOG_LEVEL_ERROR + + implicit none + + private + public :: create_pressure_solver + + contains + +!=============================================================================! + !> @brief Create iterative solver for (Helmholtz) pressure problem + !> @details Called by init method of this module, but also by + !! adjt_mixed_schur_preconditioner_alg_mod and adjt_mixed_solver_alg_mod + !> @param[in] pressure_operator_in Input (Helmholtz) pressure operator + !> @param[in] pressure_preconditioner_in Input (Helmholtz) pressure preconditioner + !> @param[out] pressure_solver_out Output (Helmholtz) pressure solver + subroutine create_pressure_solver( pressure_operator_in, pressure_preconditioner_in, pressure_solver_out ) + + use helmholtz_solver_config_mod, only: si_pressure_maximum_iterations, & + helmholtz_gcrk => gcrk, & + si_pressure_tolerance, & + si_pressure_a_tol, & + helmholtz_method => method, & + method_cg, & + method_bicgstab, & + method_gmres, & + method_fgmres, & + method_gcr, & + method_prec_only, & + method_jacobi, & + si_pressure_monitor_convergence => & + monitor_convergence, & + si_pressure_fail_on_non_converged => & + fail_on_non_converged, & + si_pressure_jacobi_relaxation => & + jacobi_relaxation + + implicit none + + ! Input operator and preconditioner for (Helmholtz) pressure problem + type(pressure_operator_type), intent(in) :: pressure_operator_in + class(abstract_preconditioner_type), allocatable, intent(in) :: pressure_preconditioner_in + + ! Output iterative solver for (Helmholtz) pressure problem + class(abstract_iterative_solver_type), allocatable, intent(out) :: pressure_solver_out + + call log_event( "create_pressure_solver: starting", LOG_LEVEL_INFO ) + + ! Allocate pressure solver of correct type + select case( helmholtz_method ) + case (METHOD_BICGSTAB) + allocate( pressure_solver_out, & + source = bicgstab_type( pressure_operator_in, & + pressure_preconditioner_in, & + si_pressure_tolerance, & + si_pressure_a_tol, & + si_pressure_maximum_iterations, & + si_pressure_monitor_convergence, & + si_pressure_fail_on_non_converged ) ) + case(METHOD_CG) + allocate( pressure_solver_out, & + source = conjugate_gradient_type( pressure_operator_in, & + pressure_preconditioner_in, & + si_pressure_tolerance, & + si_pressure_a_tol, & + si_pressure_maximum_iterations, & + si_pressure_monitor_convergence, & + si_pressure_fail_on_non_converged ) ) + case(METHOD_GMRES) + allocate( pressure_solver_out, & + source = gmres_type( pressure_operator_in, & + pressure_preconditioner_in, & + helmholtz_gcrk, & + si_pressure_tolerance, & + si_pressure_a_tol, & + si_pressure_maximum_iterations, & + si_pressure_monitor_convergence, & + si_pressure_fail_on_non_converged ) ) + case(METHOD_FGMRES) + allocate( pressure_solver_out, & + source = fgmres_type( pressure_operator_in, & + pressure_preconditioner_in, & + helmholtz_gcrk, & + si_pressure_tolerance, & + si_pressure_a_tol, & + si_pressure_maximum_iterations, & + si_pressure_monitor_convergence, & + si_pressure_fail_on_non_converged ) ) + case(METHOD_GCR) + allocate( pressure_solver_out, & + source = gcr_type( pressure_operator_in, & + pressure_preconditioner_in, & + helmholtz_gcrk, & + si_pressure_tolerance, & + si_pressure_a_tol, & + si_pressure_maximum_iterations, & + si_pressure_monitor_convergence, & + si_pressure_fail_on_non_converged ) ) + case(METHOD_PREC_ONLY) + allocate( pressure_solver_out, & + source = precondition_only_type( pressure_operator_in, & + pressure_preconditioner_in, & + si_pressure_monitor_convergence) ) + case(METHOD_JACOBI) + allocate( pressure_solver_out, & + source = jacobi_type( pressure_operator_in, & + pressure_preconditioner_in, & + si_pressure_tolerance, & + si_pressure_a_tol, & + si_pressure_maximum_iterations, & + si_pressure_monitor_convergence, & + si_pressure_fail_on_non_converged, & + si_pressure_jacobi_relaxation ) ) + case default + call log_event("Unknown pressure solver specified",LOG_LEVEL_ERROR) + end select + + call log_event( "create_pressure_solver: done", LOG_LEVEL_INFO ) + + end subroutine create_pressure_solver + +end module create_pressure_solver_mod diff --git a/science/gungho/source/algorithm/solver/mixed_operator_alg_mod.x90 b/science/gungho/source/algorithm/solver/mixed_operator_alg_mod.x90 index 72454b17..691e8942 100644 --- a/science/gungho/source/algorithm/solver/mixed_operator_alg_mod.x90 +++ b/science/gungho/source/algorithm/solver/mixed_operator_alg_mod.x90 @@ -57,7 +57,8 @@ module mixed_operator_alg_mod eliminate_variables, & eliminate_variables_analytic, & eliminate_variables_discrete - use sci_linear_operator_mod, only: abstract_linear_operator_type + use sci_hierarchical_linear_operator_mod, & + only: abstract_hierarchical_linear_operator_type use sci_r_solver_field_vector_mod, & only: r_solver_field_vector_type use split_w2_field_kernel_mod, only: split_w2_field_kernel_type @@ -69,11 +70,14 @@ module mixed_operator_alg_mod private - type, public, extends(abstract_linear_operator_type) :: & + type, public, extends(abstract_hierarchical_linear_operator_type) :: & mixed_operator_type private + !> Level of the operator (used for multigrid) + integer(kind=i_def) :: level + contains !> Over-ride the abstract interface !> param[in,out] self A linear operator @@ -82,12 +86,34 @@ module mixed_operator_alg_mod procedure, public :: apply => apply_mixed_operator !> Applies the operator as LMA matrix matrix-vector kernels in all blocks procedure, private :: apply_mixed_operator + !> Returns a version of the operator on the next-coarser space + procedure, public :: coarsen => coarsen_mixed_operator + procedure, private :: coarsen_mixed_operator !> Destroys the object final :: destroy_mixed_operator end type mixed_operator_type + interface mixed_operator_type + module procedure mixed_operator_constructor + end interface + contains + !> @brief Construct new instance of the mixed operator on a specified mesh level. + !> + !> @param[in] level The mesh level the space is on + !> @return Instance of the Helmholtz operator + function mixed_operator_constructor(level) result(self) + implicit none + + type(mixed_operator_type) :: self + integer(kind=i_def), intent(in) :: level + + ! Temporaries required in operator application + self%level = level + + end function mixed_operator_constructor + !> @brief Applies the mixed operator to the vector, \f$ y = M x \f$. !> !> @param[in,out] self Instance of the mixed operator @@ -203,11 +229,11 @@ contains mm_vel => get_w2_mass_matrix_r_solver(w2_si_matrix, mesh%get_id()) ! Obtain operators from SI operators - p2theta => get_p2theta() - div_star => get_div_star() - ptheta2 => get_ptheta2() - m3_exner_star => get_m3_exner_star() - p3theta => get_p3theta() + p2theta => get_p2theta(self%level) + div_star => get_div_star(self%level) + ptheta2 => get_ptheta2(self%level) + m3_exner_star => get_m3_exner_star(self%level) + p3theta => get_p3theta(self%level) select type (x) type is (r_solver_field_vector_type) @@ -254,7 +280,7 @@ contains select case ( eliminate_variables ) case ( eliminate_variables_discrete ) - q32_op => get_eliminated_q32() + q32_op => get_eliminated_q32(self%level) if ( optimised_operator ) then call invoke( name="apply_mixed_operator_new", & setval_c( yvec_uv, 0.0_r_solver ), & @@ -296,8 +322,8 @@ contains ! For analytic elimination: ! Q32 = tau*dt*M3^rho * D * rho^ref + ! and x_t = 0 - q32_op => get_eliminated_q32() - q22_op => get_eliminated_q22() + q32_op => get_eliminated_q32(self%level) + q22_op => get_eliminated_q22(self%level) call q22u%initialise( vector_space = u_fs ) call invoke( name="analytic_elim_mixed_lhs", & setval_c( x_t, 0.0_r_solver ), & @@ -370,4 +396,19 @@ contains type(mixed_operator_type), intent(inout) :: self end subroutine destroy_mixed_operator + !> @brief Return a coarsened version of the operator on the + !! next level of the multigrid hierarchy. + !> + !>@param[in] self Instance of pressure_operator_type + !>@param[in,out] coarse_operator Coarsened version of the operator on the next multigrid level + subroutine coarsen_mixed_operator(self, coarse_operator) + implicit none + class(mixed_operator_type), intent(in) :: self + class(abstract_hierarchical_linear_operator_type), allocatable, intent(inout) :: coarse_operator + + allocate(coarse_operator, & + source=mixed_operator_type(self%level+1)) + + end subroutine coarsen_mixed_operator + end module mixed_operator_alg_mod diff --git a/science/gungho/source/algorithm/solver/mixed_schur_preconditioner_alg_mod.x90 b/science/gungho/source/algorithm/solver/mixed_schur_preconditioner_alg_mod.x90 index e9f613b8..4bc08706 100644 --- a/science/gungho/source/algorithm/solver/mixed_schur_preconditioner_alg_mod.x90 +++ b/science/gungho/source/algorithm/solver/mixed_schur_preconditioner_alg_mod.x90 @@ -96,7 +96,8 @@ module mixed_schur_preconditioner_alg_mod get_face_selector_ns use integer_field_mod, only: integer_field_type use sci_iterative_solver_mod, only: abstract_iterative_solver_type - use sci_preconditioner_mod, only: abstract_preconditioner_type + use sci_hierarchical_preconditioner_mod, & + only: abstract_hierarchical_preconditioner_type use sci_r_solver_field_vector_mod, only: r_solver_field_vector_type use field_indices_mod, only: isol_u, isol_p, & isol_uv, isol_w @@ -106,6 +107,7 @@ module mixed_schur_preconditioner_alg_mod use log_mod, only: log_event, & LOG_LEVEL_ERROR, & LOG_LEVEL_DEBUG, & + LOG_LEVEL_INFO, & log_scratch_space use mesh_mod, only: mesh_type use si_operators_alg_mod, only: get_m3_rho_star, & @@ -122,6 +124,11 @@ module mixed_schur_preconditioner_alg_mod use boundaries_config_mod, only: limited_area use mixed_solver_config_mod, only: eliminate_variables, & eliminate_variables_analytic + use create_pressure_preconditioner_mod, only: create_pressure_preconditioner + use create_pressure_solver_mod, only: create_pressure_solver + use pressure_operator_alg_mod, only: pressure_operator_type + use sci_preconditioner_mod, only: abstract_preconditioner_type + implicit none private @@ -130,7 +137,7 @@ module mixed_schur_preconditioner_alg_mod !> !> @details Implements mixed preconditioner which can be used in the !! iterative solver algorithms. - type, public, extends(abstract_preconditioner_type) :: & + type, public, extends(abstract_hierarchical_preconditioner_type) :: & mixed_schur_preconditioner_type private @@ -141,10 +148,18 @@ module mixed_schur_preconditioner_alg_mod type(r_solver_field_vector_type) :: pressure_b !> 1-component field vector for solution of pressure system type(r_solver_field_vector_type) :: pressure_x + !> Pressure (Helmholtz) operator + type( pressure_operator_type ) :: pressure_operator + !> Pressure (Helmholtz) preconditioner + class( abstract_preconditioner_type ), allocatable :: & + pressure_preconditioner !> Pressure (Helmholtz) solver object - class(abstract_iterative_solver_type), pointer :: & + class(abstract_iterative_solver_type), allocatable :: & pressure_solver + !> Level of the operator (used for multigrid) + integer(kind=i_def) :: level + contains ! Override the (abstract interface) for application of ! a preconditioner \f$y = P.x\f$ @@ -157,6 +172,14 @@ module mixed_schur_preconditioner_alg_mod !> Reconstruct the velocity and buoyancy from the solution of the !> Helmholtz equation procedure, private :: back_substitute + !> Returns a version of the preconditioner on the next-coarser level + procedure, public :: coarsen => coarsen_mixed_schur_preconditioner + procedure, private :: coarsen_mixed_schur_preconditioner + !> Copy a mixed_schur_preconditioner_type object + procedure, public :: mixed_schur_preconditioner_assign + + !> Override default assignment for pressure_operator_type pairs. + generic :: assignment(=) => mixed_schur_preconditioner_assign !> Destructor final :: destroy_mixed_schur_preconditioner @@ -182,11 +205,10 @@ contains !> @param [in] mesh Mesh to create function spaces on !> @param [in] element_order_h Horizontal element order for function spaces !> @param [in] element_order_v Vertical element order for function spaces - !> @param [in] pressure_solver Solver object for Helmholtz system !> @return self the constructed preconditioner object function mixed_schur_preconditioner_constructor(mesh, element_order_h, & element_order_v, & - pressure_solver) & + level) & result(self) use function_space_mod, only: function_space_type @@ -198,15 +220,16 @@ contains integer(i_def), intent(in) :: element_order_h integer(i_def), intent(in) :: element_order_v type(mesh_type), target, intent(in) :: mesh - class(abstract_iterative_solver_type), target, intent(in) :: pressure_solver + integer(i_def), intent(in) :: level type(mixed_schur_preconditioner_type) :: self type(mesh_type), pointer :: mesh_ptr type(function_space_type), pointer :: fs - call log_event( 'Constructing approximate Schur mixed preconditioner...', & - LOG_LEVEL_DEBUG ) + write(log_scratch_space, '(A,I3)') & + 'mixed_schur_preconditioner_alg_mod: constructing mixed Schur preconditioner on level ', level + call log_event( log_scratch_space, LOG_LEVEL_DEBUG ) mesh_ptr => mesh fs => function_space_collection%get_fs( mesh_ptr, element_order_h, & @@ -223,13 +246,46 @@ contains call self%pressure_b%initialise_field( 1, fs ) call self%pressure_x%initialise_field( 1, fs ) - ! Set pressure solver - self%pressure_solver => pressure_solver + self%level = level + ! Create solver for the pressure system on the next multigrid level + call create_pressure_preconditioner( fs%get_id(), self%pressure_operator, & + self%pressure_preconditioner, self%level ) + call create_pressure_solver( self%pressure_operator, self%pressure_preconditioner, & + self%pressure_solver ) call log_event( 'done', LOG_LEVEL_DEBUG ) end function mixed_schur_preconditioner_constructor + !> @brief Performs a deep copy between mixed_schur_preconditioner_type pairs + !! (for overriding the "=" operator). + !> + !> @param[out] dest mixed_schur_preconditioner_type lhs + !> @param[in] source mixed_schur_preconditioner_type rhs + subroutine mixed_schur_preconditioner_assign(dest, source) + + implicit none + class(mixed_schur_preconditioner_type), intent(in) :: source + class(mixed_schur_preconditioner_type), intent(out) :: dest + + ! Deep copy of the contents of the pressure_operator_type + dest%level = source%level + + call dest%rhs_u%initialise( vector_space = source%rhs_u%get_function_space() ) + + dest%pressure_b = r_solver_field_vector_type(1) + dest%pressure_x = r_solver_field_vector_type(1) + call dest%pressure_b%initialise_field( 1, source%pressure_b%vector(1)%get_function_space() ) + call dest%pressure_x%initialise_field( 1, source%pressure_x%vector(1)%get_function_space() ) + + ! Copy the pressure solver + dest%pressure_operator = source%pressure_operator + allocate(dest%pressure_preconditioner, source = source%pressure_preconditioner) + call create_pressure_solver( dest%pressure_operator, dest%pressure_preconditioner, & + dest%pressure_solver ) + + end subroutine mixed_schur_preconditioner_assign + !> @brief Apply the preconditioner to calculate \f$y = P.x\f$ for the mixed !! system in velocity, pressure, density and potential temperature. !> @@ -250,7 +306,6 @@ contains class(abstract_vector_type), intent(inout) :: y if ( subroutine_timers ) call timer('mixed_schur_preconditioner_alg') - select type(x) type is(r_solver_field_vector_type) select type(y) @@ -362,11 +417,11 @@ contains call invoke( setval_X( rhs_uvw, rhs_uv ) ) end if - ptheta2 => get_ptheta2v() - p3theta => get_p3theta() - m3_rho_star => get_m3_rho_star() - compound_div => get_compound_div() - Hb_lumped_inv => get_Hb_lumped_inv() + ptheta2 => get_ptheta2v(self%level) + p3theta => get_p3theta(self%level) + m3_rho_star => get_m3_rho_star(self%level) + compound_div => get_compound_div(self%level) + Hb_lumped_inv => get_Hb_lumped_inv(self%level) ! Compute rhs_u terms call invoke( name = "compute_elim_helmholtz_ru", & @@ -379,7 +434,7 @@ contains ! Compute H(rhs_u) if ( eliminate_variables == eliminate_variables_analytic ) then ! r_pi = rhs_pi - Q32*rhs_u - q32_op => get_eliminated_q32() + q32_op => get_eliminated_q32(self%level) call r_p%initialise( rhs%get_function_space() ) call invoke( dg_matrix_vector_kernel_type(r_p, self%rhs_u, q32_op), & X_minus_Y(rhs, rhs_p, r_p ) ) @@ -463,8 +518,8 @@ contains ! u increment u_normalisation => get_normalisation_r_solver(W2, mesh_id) - div_star => get_div_star() - Hb_lumped_inv => get_Hb_lumped_inv() + div_star => get_div_star(self%level) + Hb_lumped_inv => get_Hb_lumped_inv(self%level) ! u' = ru + BC[HB * unorm * D * p'] if ( split_w ) then state_uv => state%get_field_from_position( isol_uv ) @@ -515,4 +570,40 @@ contains end subroutine destroy_mixed_schur_preconditioner + !> @brief Construct a coarsened version of the preconditioner on the next + !! level of the multigrid hierarchy. + !> + !> @param[in,out] self Instance of type mixed_schur_preconditioner_type + !> @param[in,out] other Coarsed version on next multigrid level + subroutine coarsen_mixed_schur_preconditioner(self, other) + use function_space_mod, only: function_space_type + use function_space_chain_mod, only: multigrid_function_space_chain + + implicit none + class(mixed_schur_preconditioner_type), intent(inout) :: self + class(abstract_hierarchical_preconditioner_type), allocatable, intent(inout) :: other + + type(function_space_type), pointer :: fs + type(r_solver_field_type), pointer :: b_vec + + b_vec => self%pressure_b%get_field_from_position(1) + fs => b_vec%get_function_space() + + ! Get the next function space in the multigrid hierarchy + call multigrid_function_space_chain%set_current( fs%get_id() ) + fs => multigrid_function_space_chain%get_next( ) + + ! Create the mixed Schur preconditioner on the next level + write(log_scratch_space, '(A,I3,A,I3)') & + "Coarsening mixed Schur preconditioner from level ", & + self%level, " to level ", self%level+1 + call log_event( log_scratch_space, LOG_LEVEL_INFO) + + other = mixed_schur_preconditioner_type( & + fs%get_mesh(), & + fs%get_element_order_h(), & + fs%get_element_order_v(), & + self%level+1) + end subroutine coarsen_mixed_schur_preconditioner + end module mixed_schur_preconditioner_alg_mod diff --git a/science/gungho/source/algorithm/solver/multigrid_preconditioner_alg_mod.x90 b/science/gungho/source/algorithm/solver/multigrid_preconditioner_alg_mod.x90 index deff82fa..18621aa9 100644 --- a/science/gungho/source/algorithm/solver/multigrid_preconditioner_alg_mod.x90 +++ b/science/gungho/source/algorithm/solver/multigrid_preconditioner_alg_mod.x90 @@ -10,11 +10,13 @@ module multigrid_preconditioner_alg_mod use base_mesh_config_mod, only: prime_mesh_name use boundaries_config_mod, only: limited_area - use constants_mod, only: i_def, r_solver, l_def + use constants_mod, only: i_def, r_solver, l_def, r_def use function_space_mod, only: function_space_type use function_space_chain_mod, & - only: function_space_chain_type, & - multigrid_function_space_chain + only: function_space_chain_type, & + multigrid_function_space_chain, & + w2h_multigrid_function_space_chain, & + w2v_multigrid_function_space_chain use log_mod, only: log_event, & LOG_LEVEL_INFO, & LOG_LEVEL_ERROR, & @@ -30,6 +32,9 @@ module multigrid_preconditioner_alg_mod only: abstract_hierarchical_preconditioner_type use sci_preconditioner_mod, only: abstract_preconditioner_type use vector_mod, only: abstract_vector_type + use field_indices_mod, only: isol_p, isol_w, isol_uv + + use psykal_lite_mod, only: print_dofs implicit none @@ -40,18 +45,20 @@ module multigrid_preconditioner_alg_mod !> on each level of the multigrid hierarchy. The concrete types are deduced !> from the type of the operator/preconditioner that is passed to the constructor !> of the multigrid type. - class(abstract_hierarchical_linear_operator_type), allocatable :: H_operator - class(abstract_hierarchical_preconditioner_type), allocatable :: Hz_preconditioner + class(abstract_hierarchical_linear_operator_type), allocatable :: Linear_operator + class(abstract_hierarchical_preconditioner_type), allocatable :: Linear_preconditioner end type linear_operators_type type, public, extends(abstract_preconditioner_type) :: multigrid_preconditioner_type private !> Number of multigrid levels integer(kind=i_def) :: n_level + !> Number of fields in the multigrid system + integer(kind=i_def) :: n_fields !> relaxation parameter real(kind=r_solver) :: rho_smooth !> Function space ids - integer(kind=i_def), dimension(:), allocatable :: fs_id + integer(kind=i_def), dimension(:,:), allocatable :: fs_id !> Mesh ids integer(kind=i_def), dimension(:), allocatable :: mesh_id !> Solution on all levels @@ -69,7 +76,7 @@ module multigrid_preconditioner_alg_mod integer(kind=i_def) :: n_postsmooth !> Jacobi smoothing parameter, for coarsest mesh integer(kind=i_def) :: n_coarsesmooth - !> Helmholtz operator and precondition for tridiagonal solve on each multigrid level + !> Linear operator and precondition for solve on each multigrid level type(linear_operators_type), allocatable, dimension(:) :: linear_operators contains !> Construct fields on all multigrid levels @@ -96,22 +103,26 @@ contains !! and calling the subroutine which constructs the multigrid fields, !! operators and preconditioners on all levels of the multigrid hierarchy. !> - !> @param[in] p_fs Function space of the pressure field - !> @param[in] H_op Pressure operator on finest level - !> @param[in] Hz_prec Pressure preconditioner on finest level + !> @param[in] fine_fs_ids IDs of the function spaces on the fine mesh + !> @param[in] n_fields Number of fields in the multigrid system + !> @param[in] MG_op Operator on finest level + !> @param[in] MG_prec Preconditioner on finest level !> @return self The constructed multigrid_preconditioner type - function multigrid_preconditioner_constructor(p_fs, & - H_op, & - Hz_prec) result(self) + function multigrid_preconditioner_constructor(fine_fs_ids, & + n_fields, & + MG_op, & + MG_prec) result(self) use multigrid_config_mod, only : smooth_relaxation, & n_presmooth, & n_postsmooth, & n_coarsesmooth, & multigrid_chain_nitems implicit none - type(function_space_type), intent(in) :: p_fs - class(abstract_hierarchical_linear_operator_type) :: H_op - class(abstract_hierarchical_preconditioner_type) :: Hz_prec + + integer(kind=i_def), intent(in) :: n_fields + integer(kind=i_def), dimension(n_fields), intent(in) :: fine_fs_ids + class(abstract_hierarchical_linear_operator_type) :: MG_op + class(abstract_hierarchical_preconditioner_type) :: MG_prec type(multigrid_preconditioner_type) :: self self%rho_smooth = real(smooth_relaxation, r_solver) @@ -119,6 +130,10 @@ contains self%n_postsmooth = n_postsmooth self%n_coarsesmooth = n_coarsesmooth self%n_level = multigrid_chain_nitems + self%n_fields = n_fields + + allocate( self%fs_id(self%n_fields, self%n_level) ) + self%fs_id(:,1) = fine_fs_ids write(log_scratch_space,'(A,I0,A,":",F3.1,3(":",I0))') & "Multigrid_preconditioner_constructor[",self%n_level,"]: setting smoothing parameters", & @@ -126,7 +141,7 @@ contains call log_event(log_scratch_space,LOG_LEVEL_INFO) ! Construct fields and operators on all multigrid levels - call self%construct_multigrid_fields(p_fs, H_op, Hz_prec) + call self%construct_multigrid_fields(MG_op, MG_prec) end function multigrid_preconditioner_constructor @@ -137,41 +152,45 @@ contains !! guarantees that they are allocated with the correct concrete type. !> !> @param[in,out] self Instance of type multigrid_preconditioner_type - !> @param[in] p_fs Function space of pressure field - !> @param[in] H_operator_fine Pressure operator on finest level - !> @param[in] Hz_preconditioner_fine Pressure preconditioner on finest level + !> @param[in] Linear_operator_fine Pressure operator on finest level + !> @param[in] Linear_preconditioner_fine Pressure preconditioner on finest level subroutine construct_multigrid_fields(self, & - p_fs, & - H_operator_fine, & - Hz_preconditioner_fine ) + Linear_operator_fine, & + Linear_preconditioner_fine ) use limited_area_constants_mod, only: get_mask_r_solver - use fs_continuity_mod, only: W3 + use fs_continuity_mod, only: W3, W2h, W2v implicit none class(multigrid_preconditioner_type), intent(inout) :: self - type(function_space_type), target, intent(in) :: p_fs - class(abstract_hierarchical_linear_operator_type), intent(in) :: H_operator_fine - class(abstract_hierarchical_preconditioner_type), intent(in) :: Hz_preconditioner_fine - type(function_space_type), pointer :: fs => null() - integer(kind=i_def) :: fs_fine_id + class(abstract_hierarchical_linear_operator_type), intent(in) :: Linear_operator_fine + class(abstract_hierarchical_preconditioner_type), intent(in) :: Linear_preconditioner_fine + + type(function_space_type), pointer :: p_fs, w_fs, uv_fs integer(kind=i_def) :: level type(r_solver_field_type), pointer :: mask => null(), & mask_mg_ptr => null() write(log_scratch_space,'(A,I1,A)') & - "construct_multigrid_preconditioner:make fields and operator on ", & + "construct_multigrid_preconditioner: make fields and operator on ", & self%n_level," levels" call log_event(log_scratch_space,LOG_LEVEL_INFO) - ! Extract ids - fs => p_fs - fs_fine_id = fs%get_id() - ! Point the function space chain to the current, fine level function space - call multigrid_function_space_chain%set_current(fs_fine_id) + !write(6,*) 'pressure function space id:', self%fs_id(isol_p, 1) + call multigrid_function_space_chain%set_current(self%fs_id(isol_p, 1)) + p_fs => multigrid_function_space_chain%get_current() + + if ( self%n_fields == 3) then + ! Multiple fields so this is a mixed system multigrid + !write(6,*) 'w2v function space id:', self%fs_id(isol_w, 1) + !write(6,*) 'w2h function space id:', self%fs_id(isol_uv, 1) + call w2v_multigrid_function_space_chain%set_current(self%fs_id(isol_w, 1)) + call w2h_multigrid_function_space_chain%set_current(self%fs_id(isol_uv, 1)) + w_fs => w2v_multigrid_function_space_chain%get_current() + uv_fs => w2h_multigrid_function_space_chain%get_current() + end if ! Allocate memory on all multigrid levels - allocate( self%fs_id(self%n_level) ) allocate( self%mesh_id(self%n_level) ) allocate( self%u_mg( self%n_level ) ) allocate( self%b_mg( self%n_level ) ) @@ -179,43 +198,66 @@ contains allocate( self%linear_operators(self%n_level) ) allocate( self%mask_mg( self%n_level ) ) - ! Set function space- and mesh ids as well as fields on all levels - self%fs_id(1) = fs_fine_id ! Create operators on the finest level; the concrete types will be deduced from ! the passed operator and preconditioner - allocate(self%linear_operators(1)%H_operator, & - source=H_operator_fine) - allocate(self%linear_operators(1)%Hz_preconditioner, & - source=Hz_preconditioner_fine) + allocate(self%linear_operators(1)%Linear_operator, & + source=Linear_operator_fine) + allocate(self%linear_operators(1)%Linear_preconditioner, & + source=Linear_preconditioner_fine) do level = 1, self%n_level - self%fs_id(level) = fs%get_id() - self%mesh_id(level) = fs%get_mesh_id() + !write(6,*) 'Creating MG fields on level ',level,' with nfields = ',self%n_fields + self%fs_id(isol_p, level) = p_fs%get_id() + self%mesh_id(level) = p_fs%get_mesh_id() ! Construct fields on a particular level - self%u_mg(level) = r_solver_field_vector_type( 1 ) - self%b_mg(level) = r_solver_field_vector_type( 1 ) - self%r_mg(level) = r_solver_field_vector_type( 1 ) - call self%u_mg(level)%initialise_field( 1, fs ) - call self%b_mg(level)%initialise_field( 1, fs ) - call self%r_mg(level)%initialise_field( 1, fs ) + self%u_mg(level) = r_solver_field_vector_type( self%n_fields ) + self%b_mg(level) = r_solver_field_vector_type( self%n_fields ) + self%r_mg(level) = r_solver_field_vector_type( self%n_fields ) + call self%u_mg(level)%initialise_field( isol_p, p_fs ) + call self%b_mg(level)%initialise_field( isol_p, p_fs ) + call self%r_mg(level)%initialise_field( isol_p, p_fs ) + if ( self%n_fields == 3 ) then + ! Multiple fields so this is a mixed system multigrid + self%fs_id(isol_w, level) = w_fs%get_id() + self%fs_id(isol_uv, level) = uv_fs%get_id() + call self%u_mg(level)%initialise_field( isol_w, w_fs ) + call self%u_mg(level)%initialise_field( isol_uv, uv_fs ) + call self%b_mg(level)%initialise_field( isol_w, w_fs ) + call self%b_mg(level)%initialise_field( isol_uv, uv_fs ) + call self%r_mg(level)%initialise_field( isol_w, w_fs ) + call self%r_mg(level)%initialise_field( isol_uv, uv_fs ) + end if if (level > 1) then ! Coarsen linear operator - call self%linear_operators(level-1)%H_operator%coarsen( & - self%linear_operators(level)%H_operator ) + call self%linear_operators(level-1)%Linear_operator%coarsen( & + self%linear_operators(level)%Linear_operator ) ! Coarsen preconditioner - call self%linear_operators(level-1)%Hz_preconditioner%coarsen( & - self%linear_operators(level)%Hz_preconditioner) + call self%linear_operators(level-1)%Linear_preconditioner%coarsen( & + self%linear_operators(level)%Linear_preconditioner) end if if (limited_area) then - self%mask_mg(level) = r_solver_field_vector_type( 1 ) - call self%mask_mg(level)%initialise_field( 1,fs ) + self%mask_mg(level) = r_solver_field_vector_type( self%n_fields ) + call self%mask_mg(level)%initialise_field( self%n_fields, p_fs ) mask => get_mask_r_solver( W3, self%mesh_id(level), prime_mesh_name ) - mask_mg_ptr => self%mask_mg(level)%get_field_from_position(1) + mask_mg_ptr => self%mask_mg(level)%get_field_from_position(isol_p) call invoke( setval_X( mask_mg_ptr, mask )) + if ( self%n_fields == 3 ) then + ! Multiple fields so this is a mixed system multigrid + mask => get_mask_r_solver( W2v, self%mesh_id(level) ) + mask_mg_ptr => self%mask_mg(level)%get_field_from_position(isol_w) + call invoke( setval_X( mask_mg_ptr, mask )) + mask => get_mask_r_solver( W2h, self%mesh_id(level) ) + mask_mg_ptr => self%mask_mg(level)%get_field_from_position(isol_uv) + call invoke( setval_X( mask_mg_ptr, mask )) + end if end if ! If we are not yet on the coarsest level point the function space ! chain to the function space on the next coarser level if ( level < self%n_level ) & - fs => multigrid_function_space_chain%get_next() + p_fs => multigrid_function_space_chain%get_next() + if ( self%n_fields == 3 ) then + w_fs => w2v_multigrid_function_space_chain%get_next() + uv_fs => w2h_multigrid_function_space_chain%get_next() + end if end do end subroutine construct_multigrid_fields @@ -228,16 +270,18 @@ contains !> @param[in] level Current multigrid level !> @param[in] initial_solution_is_zero Can we implicitly assume that the !! solution is zero? - recursive subroutine vcycle(self,level,initial_solution_is_zero) + recursive subroutine vcycle(self, level, initial_solution_is_zero) use sci_restrict_scalar_unweighted_kernel_mod, & - only: restrict_scalar_unweighted_kernel_type + only: restrict_scalar_unweighted_kernel_type + use sci_restrict_w2h_kernel_mod, only: restrict_w2h_kernel_type use sci_prolong_scalar_unweighted_kernel_mod, & - only: prolong_scalar_unweighted_kernel_type + only: prolong_scalar_unweighted_kernel_type use sci_restrict_scalar_masked_kernel_mod, & - only: restrict_scalar_masked_kernel_type + only: restrict_scalar_masked_kernel_type use sci_prolong_scalar_masked_kernel_mod, & - only: prolong_scalar_masked_kernel_type + only: prolong_scalar_masked_kernel_type + use intermesh_mappings_alg_mod, only: map_w2h_intermesh_rsolver implicit none @@ -251,53 +295,77 @@ contains bvec_mg_coarse => null(), & uvec_mg_coarse => null(), & mask => null() - + integer(kind=i_def) :: n + real(kind=r_def) :: error ! Call smoother on coarsest level if (level == self%n_level) then - call self%smooth(self%n_coarsesmooth,level,.true.) + call self%smooth(self%n_coarsesmooth, level, .true.) else ! Presmooth solution u^{(ell)} call self%smooth(self%n_presmooth, level, initial_solution_is_zero) ! Calculate residual r^{(ell)} = b^{(ell)} - H.u^{(ell)} - call self%linear_operators(level)%H_operator%apply(self%u_mg(level), & - self%r_mg(level) ) - - rvec_mg => self%r_mg(level)%get_field_from_position(1) - bvec_mg => self%b_mg(level)%get_field_from_position(1) - bvec_mg_coarse => self%b_mg(level+1)%get_field_from_position(1) - - call invoke( inc_aX_plus_y (-1.0_r_solver, rvec_mg, bvec_mg) ) - ! Restrict b^{(ell+1)} = R r^{(ell)} - if (limited_area) then - mask => self%mask_mg(level)%get_field_from_position(1) - call invoke( name='masked_restriction', & - restrict_scalar_masked_kernel_type(bvec_mg_coarse, & - rvec_mg, & - mask )) - else - call invoke( name='unweighted_restriction', & - restrict_scalar_unweighted_kernel_type(bvec_mg_coarse, rvec_mg) ) - - end if + call self%linear_operators(level)%Linear_operator%apply(self%u_mg(level), & + self%r_mg(level) ) + + do n = 1, self%n_fields + rvec_mg => self%r_mg(level)%get_field_from_position(n) + bvec_mg => self%b_mg(level)%get_field_from_position(n) + bvec_mg_coarse => self%b_mg(level+1)%get_field_from_position(n) + + call invoke( inc_aX_plus_y (-1.0_r_solver, rvec_mg, bvec_mg) ) + + ! Restrict b^{(ell+1)} = R r^{(ell)} + if (limited_area) then + mask => self%mask_mg(level)%get_field_from_position(n) + call invoke( name='masked_restriction', & + restrict_scalar_masked_kernel_type(bvec_mg_coarse, & + rvec_mg, & + mask )) + else + if ( n == isol_uv ) then + call map_w2h_intermesh_rsolver(bvec_mg_coarse, rvec_mg) + else + call invoke( name='unweighted_restriction', & + restrict_scalar_unweighted_kernel_type(bvec_mg_coarse, rvec_mg) ) + + end if + end if + call invoke( X_innerproduct_X(error, rvec_mg) ) + write(6,*) 'Pre Restriction: ',level, n, error + call invoke( X_innerproduct_X(error, bvec_mg_coarse) ) + write(6,*) 'Post Restriction: ',level, n, error + end do ! Recursive call to V-cycle call self%vcycle(level+1,.true.) ! Prolong and add solution from coarser level ! u^{(ell)} = u^{(ell)} + P u^{(ell+1)} - uvec_mg => self%u_mg(level)%get_field_from_position(1) - uvec_mg_coarse => self%u_mg(level+1)%get_field_from_position(1) - if (limited_area) then - mask => self%mask_mg(level)%get_field_from_position(1) - call invoke( name='masked_prolong', & - prolong_scalar_masked_kernel_type(uvec_mg, & - uvec_mg_coarse, & - mask ) ) - else - call invoke( name='unweighted_prolong', & - prolong_scalar_unweighted_kernel_type(uvec_mg, uvec_mg_coarse) ) - endif + do n = 1, self%n_fields + uvec_mg => self%u_mg(level)%get_field_from_position(n) + uvec_mg_coarse => self%u_mg(level+1)%get_field_from_position(n) + + call print_dofs(uvec_mg_coarse, 'Coarse field', 12) + if (limited_area) then + mask => self%mask_mg(level)%get_field_from_position(n) + call invoke( name='masked_prolong', & + prolong_scalar_masked_kernel_type(uvec_mg, & + uvec_mg_coarse, & + mask ) ) + else + if ( n == isol_uv ) then + call map_w2h_intermesh_rsolver(uvec_mg, uvec_mg_coarse) + else + call invoke( name='unweighted_prolong', & + prolong_scalar_unweighted_kernel_type(uvec_mg, uvec_mg_coarse) ) + end if + endif + call invoke( X_innerproduct_X(error, uvec_mg_coarse) ) + write(6,*) 'Pre Prolongation: ',level, n, error + call invoke( X_innerproduct_X(error, uvec_mg) ) + write(6,*) 'Post Prolongation: ',level, n, error + end do ! Postsmooth u^{(ell)} call self%smooth(self%n_postsmooth, level, .false.) @@ -323,20 +391,35 @@ contains y_vec => null(), & bvec_mg => null(), & uvec_mg => null() + + integer(kind=i_def) :: n + + call log_event('Applying multigrid preconditioner', log_level_info) + select type(x) type is(r_solver_field_vector_type) select type(y) type is(r_solver_field_vector_type) ! Copy in RHS on finest level - x_vec => x%get_field_from_position(1) - bvec_mg => self%b_mg(1)%get_field_from_position(1) - call invoke(setval_X(bvec_mg, x_vec )) + do n = 1, self%n_fields + !write(6,*) 'copy x_vec: ', n + x_vec => x%get_field_from_position(n) + bvec_mg => self%b_mg(1)%get_field_from_position(n) + !fs => x_vec%get_function_space() + !write(6,*) 'x_fs undf = ', fs%get_undf() + !fs => bvec_mg%get_function_space() + !write(6,*) 'b_fs undf = ', fs%get_undf() + !flush(6) + call invoke(setval_X(bvec_mg, x_vec )) + end do ! Solve using a multigrid V-cycle call self%vcycle(1, .true.) ! Copy out solution on finest level - y_vec => y%get_field_from_position(1) - uvec_mg => self%u_mg(1)%get_field_from_position(1) - call invoke(setval_X(y_vec, uvec_mg )) + do n = 1, self%n_fields + y_vec => y%get_field_from_position(n) + uvec_mg => self%u_mg(1)%get_field_from_position(n) + call invoke(setval_X(y_vec, uvec_mg )) + end do class default write(log_scratch_space, '(A)') & "multigrid_preconditioner_mod: incorrect vector_type argument y" @@ -367,34 +450,75 @@ contains !! \f$ u^{(\ell)} \f$ is zero before !! the first iteration? subroutine smooth(self,nsmooth,level,initial_solution_is_zero) + use mixed_solver_config_mod, only: mixed_solver_a_tol, & + si_tolerance, & + mixed_gcrk => gcrk, & + si_monitor_convergence => & + monitor_convergence, & + si_fail_on_non_converged => & + fail_on_non_converged, & + si_maximum_iterations + use sci_iterative_solver_mod, only: block_gcr_type, & + abstract_iterative_solver_type + implicit none class(multigrid_preconditioner_type), intent(inout) :: self integer(kind=i_def), intent(in) :: nsmooth integer(kind=i_def), intent(in) :: level logical(kind=l_def), intent(in) :: initial_solution_is_zero - type(function_space_type), pointer :: fs => null() + type(function_space_type), pointer :: fs type(r_solver_field_vector_type) :: res - type(r_solver_field_type), pointer :: uvec_mg => null(), & - rvec_mg => null(), & - bvec_mg => null() + type(r_solver_field_type), pointer :: uvec_mg, & + rvec_mg, & + bvec_mg real(kind=r_solver) :: const - integer(kind=i_def) :: i, i_first + real(kind=r_def), allocatable :: error(:), error0(:) + integer(kind=i_def) :: i, i_first, n + + class( abstract_iterative_solver_type ), allocatable :: coarse_level_solver + + + if ( level == self%n_level .and. .false. ) then + ! Iterative solve on coarsest level + allocate( coarse_level_solver, & + source = block_gcr_type( self%linear_operators(level)%Linear_operator, & + self%linear_operators(level)%Linear_preconditioner, & + mixed_gcrk, & + si_tolerance, & + mixed_solver_a_tol, & + si_maximum_iterations, & + si_monitor_convergence, & + si_fail_on_non_converged ) ) + call self%u_mg(level)%set_scalar(0.0_r_def) + call coarse_level_solver%apply(self%u_mg(level), self%b_mg(level)) + else + + allocate( error(self%n_fields), error0(self%n_fields) ) ! Construct temporary fields (we probably do not want to construct them ! in every smoother application in the future) - fs => self%u_mg(level)%vector(1)%get_function_space() - res = r_solver_field_vector_type( 1 ) - call res%initialise_field( 1, fs ) + res = r_solver_field_vector_type( self%n_fields ) + do i = 1, self%n_fields + uvec_mg => self%u_mg(level)%get_field_from_position(i) + fs => uvec_mg%get_function_space() + call res%initialise_field( i, fs ) + end do + do n = 1, self%n_fields + rvec_mg => self%b_mg(level)%get_field_from_position(n) + call invoke( X_innerproduct_X(error0(n), rvec_mg) ) + end do if (initial_solution_is_zero) then ! First iteration - call self%linear_operators(level)%Hz_preconditioner%apply(self%b_mg(level), & - self%r_mg(level)) + call self%linear_operators(level)%Linear_preconditioner%apply(self%b_mg(level), & + self%r_mg(level)) const = self%rho_smooth - uvec_mg => self%u_mg(level)%get_field_from_position(1) - rvec_mg => self%r_mg(level)%get_field_from_position(1) - call invoke(a_times_X(uvec_mg, const, rvec_mg)) + do n = 1, self%n_fields + uvec_mg => self%u_mg(level)%get_field_from_position(n) + rvec_mg => self%r_mg(level)%get_field_from_position(n) + call invoke(a_times_X(uvec_mg, const, rvec_mg)) + end do i_first = 2 else i_first = 1 @@ -405,18 +529,35 @@ contains const = -self%rho_smooth do i = i_first, nsmooth ! Calculate res = H.u^{(ell)} - b^{(ell)} - call self%linear_operators(level)%H_operator%apply(self%u_mg(level), res) - bvec_mg => self%b_mg(level)%get_field_from_position(1) - rvec_mg => res%get_field_from_position(1) - call invoke( inc_X_minus_Y (rvec_mg, bvec_mg)) + call self%linear_operators(level)%Linear_operator%apply(self%u_mg(level), res) + do n = 1, self%n_fields + bvec_mg => self%b_mg(level)%get_field_from_position(n) + rvec_mg => res%get_field_from_position(n) + call invoke( inc_X_minus_Y (rvec_mg, bvec_mg)) + end do ! Calculate r^{(ell)} = H_z^{-1}.(H.u^{(ell)} - b^{(ell)}) - call self%linear_operators(level)%Hz_preconditioner%apply(res, self%r_mg(level)) + call self%linear_operators(level)%Linear_preconditioner%apply(res, self%r_mg(level)) ! Calculate u^{(ell)} = u^{(ell)} - rho*H_z^{-1}.(H.u^{(ell)} - b^{(ell)}) - uvec_mg => self%u_mg(level)%get_field_from_position(1) - rvec_mg => self%r_mg(level)%get_field_from_position(1) - call invoke(inc_X_plus_bY (uvec_mg, const, rvec_mg) ) + do n = 1, self%n_fields + uvec_mg => self%u_mg(level)%get_field_from_position(n) + rvec_mg => self%r_mg(level)%get_field_from_position(n) + call invoke(inc_X_plus_bY (uvec_mg, const, rvec_mg) ) + rvec_mg => res%get_field_from_position(n) + call invoke( X_innerproduct_X(error(n), rvec_mg) ) + call print_dofs(rvec_mg, 'Residual field', 12) + call print_dofs(uvec_mg, 'Smoothed field', 12) + end do + if ( self%n_fields == 1 ) then + write(log_scratch_space,'(A,2I4,2E16.8)') 'Smoothing Error: ',level,i,error0, error!/error0 + else + write(log_scratch_space,'(A,2I4,6E16.8)') 'Smoothing Error: ',level,i,error0, error!/error0 + end if + call log_event(log_scratch_space, LOG_LEVEL_INFO) end do + + end if + end subroutine smooth !> @brief Finalizer for the multigrid preconditioner. @@ -429,8 +570,8 @@ contains if (allocated(self%linear_operators)) then do level = 1, self%n_level - deallocate(self%linear_operators(level)%H_operator) - deallocate(self%linear_operators(level)%Hz_preconditioner) + deallocate(self%linear_operators(level)%Linear_operator) + deallocate(self%linear_operators(level)%Linear_preconditioner) end do deallocate(self%linear_operators) end if diff --git a/science/gungho/source/algorithm/solver/pressure_precon_alg_mod.x90 b/science/gungho/source/algorithm/solver/pressure_precon_alg_mod.x90 index f5d1260e..740f523b 100644 --- a/science/gungho/source/algorithm/solver/pressure_precon_alg_mod.x90 +++ b/science/gungho/source/algorithm/solver/pressure_precon_alg_mod.x90 @@ -112,6 +112,7 @@ contains type(r_solver_field_type), pointer :: tri(:) => null() type(r_solver_field_type), pointer :: x_vec => null(), & y_vec => null() + select type(x) type is(r_solver_field_vector_type) select type(y) diff --git a/science/gungho/source/algorithm/solver/semi_implicit_solver_alg_mod.x90 b/science/gungho/source/algorithm/solver/semi_implicit_solver_alg_mod.x90 index 893deb5b..40f9b081 100644 --- a/science/gungho/source/algorithm/solver/semi_implicit_solver_alg_mod.x90 +++ b/science/gungho/source/algorithm/solver/semi_implicit_solver_alg_mod.x90 @@ -39,6 +39,9 @@ module semi_implicit_solver_alg_mod igh_p, igh_t, igh_d, & igh_u, igh_w, igh_uv use copy_field_alg_mod, only: copy_field + use function_space_mod, only: function_space_type + use function_space_collection_mod, only: function_space_collection + use fs_continuity_mod, only: W2v, W2h ! Pointers use sci_geometric_constants_mod, only: get_face_selector_ew, & @@ -47,9 +50,6 @@ module semi_implicit_solver_alg_mod ! Algorithms use mixed_operator_alg_mod, only: mixed_operator_type use mixed_schur_preconditioner_alg_mod, only: mixed_schur_preconditioner_type - use pressure_operator_alg_mod, only: pressure_operator_type - use pressure_precon_alg_mod, only: pressure_preconditioner_type - use pressure_diag_precon_alg_mod, only: pressure_diag_preconditioner_type use multigrid_preconditioner_alg_mod, only: multigrid_preconditioner_type use sci_null_preconditioner_alg_mod, only: null_preconditioner_type @@ -82,14 +82,6 @@ module semi_implicit_solver_alg_mod class( abstract_preconditioner_type ), allocatable :: mixed_preconditioner class( abstract_iterative_solver_type ), allocatable :: mixed_solver - !> Operator, preconditioner and iterative solver for - !> Helmholtz (pressure) problem - type( pressure_operator_type ) :: pressure_operator - class( abstract_preconditioner_type ), allocatable :: pressure_preconditioner - class( abstract_iterative_solver_type ), allocatable :: pressure_solver - - public :: create_pressure_preconditioner - public :: create_pressure_solver public :: create_mixed_preconditioner public :: create_mixed_solver public :: semi_implicit_solver_alg_init @@ -97,185 +89,18 @@ module semi_implicit_solver_alg_mod public :: semi_implicit_solver_alg_step private :: construct_solver_state contains -!=============================================================================! - !> @brief Create operator and preconditioner for (Helmholtz) pressure problem - !> @details Called by init method of this module, but also by adj_semi_implicit_solver_alg_mod, - !! adjt_mixed_schur_preconditioner_alg_mod and adjt_mixed_solver_alg_mod - !> @param[in] state Prognostic state for the solver - !> @param[out] pressure_operator_out Output (Helmholtz) pressure operator - !> @param[out] pressure_preconditioner_out Output (Helmholtz) pressure preconditioner - subroutine create_pressure_preconditioner( state, pressure_operator_out, pressure_preconditioner_out ) - - use helmholtz_solver_config_mod, only: helmholtz_preconditioner => preconditioner, & - preconditioner_none, & - preconditioner_diagonal, & - preconditioner_tridiagonal, & - preconditioner_multigrid - - implicit none - - ! Prognostic fields - type(field_type), dimension(bundle_size), intent(in) :: state - - ! Output operator and preconditioner for (Helmholtz) pressure problem - type(pressure_operator_type), intent(out) :: pressure_operator_out - class(abstract_preconditioner_type), allocatable, intent(out) :: pressure_preconditioner_out - - ! Vertical pressure preconditioner - type(pressure_preconditioner_type) :: Hz_preconditioner - - pressure_operator_out = pressure_operator_type(level=1_i_def) - - call log_event( "create_pressure_preconditioner: starting", LOG_LEVEL_INFO ) - - ! Allocate pressure preconditioner of correct type - select case(helmholtz_preconditioner) - case(PRECONDITIONER_NONE) - allocate( pressure_preconditioner_out, & - source = null_preconditioner_type() ) - case(PRECONDITIONER_DIAGONAL) - allocate( pressure_preconditioner_out, & - source = pressure_diag_preconditioner_type() ) - case(PRECONDITIONER_TRIDIAGONAL) - allocate( pressure_preconditioner_out, & - source = pressure_preconditioner_type(level=1_i_def) ) - case(PRECONDITIONER_MULTIGRID) - Hz_preconditioner = pressure_preconditioner_type(level=1_i_def) - allocate( pressure_preconditioner_out, & - source = multigrid_preconditioner_type( state(igh_p)%get_function_space(), & - pressure_operator_out, & - Hz_preconditioner ) ) - case default - call log_event( "Unknown pressure preconditioner specified", LOG_LEVEL_ERROR) - end select - - call log_event( "create_pressure_preconditioner: done", LOG_LEVEL_INFO ) - - end subroutine create_pressure_preconditioner - -!=============================================================================! - !> @brief Create iterative solver for (Helmholtz) pressure problem - !> @details Called by init method of this module, but also by - !! adjt_mixed_schur_preconditioner_alg_mod and adjt_mixed_solver_alg_mod - !> @param[in] pressure_operator_in Input (Helmholtz) pressure operator - !> @param[in] pressure_preconditioner_in Input (Helmholtz) pressure preconditioner - !> @param[out] pressure_solver_out Output (Helmholtz) pressure solver - subroutine create_pressure_solver( pressure_operator_in, pressure_preconditioner_in, pressure_solver_out ) - - use helmholtz_solver_config_mod, only: si_pressure_maximum_iterations, & - helmholtz_gcrk => gcrk, & - si_pressure_tolerance, & - si_pressure_a_tol, & - helmholtz_method => method, & - method_cg, & - method_bicgstab, & - method_gmres, & - method_fgmres, & - method_gcr, & - method_prec_only, & - method_jacobi, & - si_pressure_monitor_convergence => & - monitor_convergence, & - si_pressure_fail_on_non_converged => & - fail_on_non_converged, & - si_pressure_jacobi_relaxation => & - jacobi_relaxation - - implicit none - - ! Input operator and preconditioner for (Helmholtz) pressure problem - type(pressure_operator_type), intent(in) :: pressure_operator_in - class(abstract_preconditioner_type), allocatable, intent(in) :: pressure_preconditioner_in - - ! Output iterative solver for (Helmholtz) pressure problem - class(abstract_iterative_solver_type), allocatable, intent(out) :: pressure_solver_out - - call log_event( "create_pressure_solver: starting", LOG_LEVEL_INFO ) - - ! Allocate pressure solver of correct type - select case( helmholtz_method ) - case (METHOD_BICGSTAB) - allocate( pressure_solver_out, & - source = bicgstab_type( pressure_operator_in, & - pressure_preconditioner_in, & - si_pressure_tolerance, & - si_pressure_a_tol, & - si_pressure_maximum_iterations, & - si_pressure_monitor_convergence, & - si_pressure_fail_on_non_converged ) ) - case(METHOD_CG) - allocate( pressure_solver_out, & - source = conjugate_gradient_type( pressure_operator_in, & - pressure_preconditioner_in, & - si_pressure_tolerance, & - si_pressure_a_tol, & - si_pressure_maximum_iterations, & - si_pressure_monitor_convergence, & - si_pressure_fail_on_non_converged ) ) - case(METHOD_GMRES) - allocate( pressure_solver_out, & - source = gmres_type( pressure_operator_in, & - pressure_preconditioner_in, & - helmholtz_gcrk, & - si_pressure_tolerance, & - si_pressure_a_tol, & - si_pressure_maximum_iterations, & - si_pressure_monitor_convergence, & - si_pressure_fail_on_non_converged ) ) - case(METHOD_FGMRES) - allocate( pressure_solver_out, & - source = fgmres_type( pressure_operator_in, & - pressure_preconditioner_in, & - helmholtz_gcrk, & - si_pressure_tolerance, & - si_pressure_a_tol, & - si_pressure_maximum_iterations, & - si_pressure_monitor_convergence, & - si_pressure_fail_on_non_converged ) ) - case(METHOD_GCR) - allocate( pressure_solver_out, & - source = gcr_type( pressure_operator_in, & - pressure_preconditioner_in, & - helmholtz_gcrk, & - si_pressure_tolerance, & - si_pressure_a_tol, & - si_pressure_maximum_iterations, & - si_pressure_monitor_convergence, & - si_pressure_fail_on_non_converged ) ) - case(METHOD_PREC_ONLY) - allocate( pressure_solver_out, & - source = precondition_only_type( pressure_operator_in, & - pressure_preconditioner_in, & - si_pressure_monitor_convergence) ) - case(METHOD_JACOBI) - allocate( pressure_solver_out, & - source = jacobi_type( pressure_operator_in, & - pressure_preconditioner_in, & - si_pressure_tolerance, & - si_pressure_a_tol, & - si_pressure_maximum_iterations, & - si_pressure_monitor_convergence, & - si_pressure_fail_on_non_converged, & - si_pressure_jacobi_relaxation ) ) - case default - call log_event("Unknown pressure solver specified",LOG_LEVEL_ERROR) - end select - - call log_event( "create_pressure_solver: done", LOG_LEVEL_INFO ) - - end subroutine create_pressure_solver !=============================================================================! !> @brief Create preconditioner for mixed problem !> @details Called by init method of this module, but also by !! adjt_mixed_schur_preconditioner_alg_mod and adjt_mixed_solver_alg_mod !> @param[in] state Prognostic state for the solver - !> @param[in] pressure_solver_in Input (Helmholtz) pressure solver !> @param[out] mixed_preconditioner_out Output mixed preconditioner - subroutine create_mixed_preconditioner( state, pressure_solver_in, mixed_preconditioner_out ) + subroutine create_mixed_preconditioner( state, mixed_preconditioner_out ) - use mixed_solver_config_mod, only: si_preconditioner, & - si_preconditioner_pressure, & + use mixed_solver_config_mod, only: si_preconditioner, & + si_preconditioner_pressure, & + si_preconditioner_multigrid, & si_preconditioner_none implicit none @@ -283,12 +108,17 @@ contains ! Prognostic fields type(field_type), dimension(bundle_size), intent(in) :: state - ! Input solver for (Helmholtz) pressure problem - class(abstract_iterative_solver_type), allocatable, intent(in) :: pressure_solver_in - ! Output preconditioner for mixed problem class(abstract_preconditioner_type), allocatable, intent(out) :: mixed_preconditioner_out + type(mixed_schur_preconditioner_type) :: mixed_multigrid_preconditioner + type(mixed_operator_type) :: mixed_multigrid_operator + integer(kind=i_def) :: nfields + type(function_space_type), pointer :: fs + integer(kind=i_def) :: fs_ids(3) + integer(kind=i_def) :: element_order_h, element_order_v + type(mesh_type), pointer :: mesh + call log_event( "create_mixed_preconditioner: starting", LOG_LEVEL_INFO ) ! Allocate mixed preconditioner of correct type @@ -299,8 +129,35 @@ contains state(igh_u)%get_mesh(), & state(igh_u)%get_element_order_h(), & state(igh_u)%get_element_order_v(), & - pressure_solver_in ) ) + level = 1_i_def ) ) + + case(SI_PRECONDITIONER_MULTIGRID) + mixed_multigrid_preconditioner = mixed_schur_preconditioner_type( & + state(igh_u)%get_mesh(), & + state(igh_u)%get_element_order_h(), & + state(igh_u)%get_element_order_v(), & + level = 1_i_def ) + mixed_multigrid_operator = mixed_operator_type(level=1_i_def) + nfields = 3_i_def + ! Function space IDs for the mixed operator (W3, W2H, W2V) + fs => state(igh_p)%get_function_space() + fs_ids(isol_p) = fs%get_id() + element_order_h = state(igh_u)%get_element_order_h() + element_order_v = state(igh_u)%get_element_order_v() + mesh => state(igh_u)%get_mesh() + fs => function_space_collection%get_fs( mesh, element_order_h, & + element_order_v, W2h ) + fs_ids(isol_uv) = fs%get_id() + fs => function_space_collection%get_fs( mesh, element_order_h, & + element_order_v, W2v ) + fs_ids(isol_w) = fs%get_id() + allocate( mixed_preconditioner_out, & + source = multigrid_preconditioner_type( fs_ids, & + nfields, & + mixed_multigrid_operator, & + mixed_multigrid_preconditioner ) ) + nullify(fs, mesh) case(SI_PRECONDITIONER_NONE) allocate( mixed_preconditioner_out, & source = null_preconditioner_type() ) @@ -350,6 +207,8 @@ contains call log_event( "create_mixed_solver: starting", LOG_LEVEL_INFO ) + mixed_operator_out = mixed_operator_type(level=1_i_def) + ! Allocate mixed solver of correct type select case(si_method) case(SI_METHOD_BICGSTAB) @@ -443,14 +302,15 @@ contains ! Prognostic fields type( field_type ), dimension(bundle_size), intent( in ) :: state + type( function_space_type ), pointer :: p_fs + write(log_scratch_space,'(A)') & 'SI solver built with '//trim(PRECISION_R_SOLVER)// & '-bit real numbers' call log_event( log_scratch_space, LOG_LEVEL_ALWAYS ) - call create_pressure_preconditioner( state, pressure_operator, pressure_preconditioner ) - call create_pressure_solver( pressure_operator, pressure_preconditioner, pressure_solver ) - call create_mixed_preconditioner( state, pressure_solver, mixed_preconditioner ) + p_fs => state(igh_p)%get_function_space() + call create_mixed_preconditioner( state, mixed_preconditioner ) call create_mixed_solver( mixed_preconditioner, mixed_operator, mixed_solver ) call log_event( "semi_implicit_solver_alg_init: Initialised semi-implicit solver", LOG_LEVEL_INFO ) @@ -472,14 +332,6 @@ contains if (allocated(mixed_solver)) then deallocate(mixed_solver) end if - ! Deallocate pressure preconditioner object - if (allocated(pressure_preconditioner)) then - deallocate(pressure_preconditioner) - end if - ! Deallocate pressure solver object - if (allocated(pressure_solver)) then - deallocate(pressure_solver) - end if end subroutine semi_implicit_solver_alg_final @@ -752,8 +604,6 @@ contains !>@param[in] bundle_state Field bundle to copy !>@param[in] import fields Import (copy) data from field_bundle subroutine construct_solver_state(vector_state, bundle_state, import_fields) - use function_space_collection_mod, only: function_space_collection - use fs_continuity_mod, only: W2v, W2h implicit none diff --git a/science/gungho/source/algorithm/solver/si_operators_alg_mod.x90 b/science/gungho/source/algorithm/solver/si_operators_alg_mod.x90 index 0b116dad..42b62943 100644 --- a/science/gungho/source/algorithm/solver/si_operators_alg_mod.x90 +++ b/science/gungho/source/algorithm/solver/si_operators_alg_mod.x90 @@ -67,7 +67,7 @@ module si_operators_alg_mod type(r_solver_operator_type), allocatable, dimension(:), target :: compound_div type(r_solver_operator_type), allocatable, dimension(:), target :: eliminated_q22 type(r_solver_operator_type), allocatable, dimension(:), target :: eliminated_q32 - type(r_solver_operator_type), target :: eliminated_q2t + type(r_solver_operator_type), allocatable, dimension(:), target :: eliminated_q2t type(r_solver_field_type), allocatable, dimension(:), target :: rho_at_u type(r_solver_field_type), allocatable, dimension(:,:), target :: tri_precon @@ -219,6 +219,7 @@ contains allocate(compound_div(multigrid_levels)) allocate(eliminated_q22(multigrid_levels)) allocate(eliminated_q32(multigrid_levels)) + allocate(eliminated_q2t(multigrid_levels)) !fields allocate(rho_at_u(multigrid_levels)) @@ -242,8 +243,6 @@ contains call wtheta_multigrid_function_space_chain%set_current(wt_fs%get_id()) end if - call eliminated_q2t%initialise( w2_fs, wt_fs ) - do i = 1, multigrid_levels write(log_scratch_space,'(A,I0,A)') "si_ops[",i,"]:creating ops" @@ -259,6 +258,7 @@ contains call compound_div(i)%initialise( w3_fs, w2_fs ) call eliminated_q22(i)%initialise( w2_fs, w2_fs ) call eliminated_q32(i)%initialise( w3_fs, w2_fs ) + call eliminated_q2t(i)%initialise( w2_fs, wt_fs ) call rho_at_u(i)%initialise(vector_space = w2_fs) if ( preconditioner == preconditioner_tridiagonal .or. & @@ -352,6 +352,10 @@ contains deallocate(eliminated_q32) end if + if(allocated(eliminated_q2t)) then + deallocate(eliminated_q2t) + end if + end subroutine final_si_operators !> @brief Subroutine to compute the si operators @@ -638,8 +642,8 @@ contains div, rho_at_u(level), const3), & operator_setval_x_kernel_type( eliminated_q32(level), compound_div(level) ) ) - if ( level == 1 .and. eliminate_variables == eliminate_variables_analytic ) then - call invoke( eliminated_theta_q2t_kernel_type(eliminated_q2t, u_normalisation, & + if ( eliminate_variables == eliminate_variables_analytic ) then + call invoke( eliminated_theta_q2t_kernel_type(eliminated_q2t(level), u_normalisation, & exner_in_wth, const1, qr) ) end if @@ -1130,10 +1134,16 @@ contains !> @brief Function to return a pointer to the q2t operator. !> @return The operator - function get_eliminated_q2t() result(op) + function get_eliminated_q2t(level) result(op) implicit none + integer(kind=i_def), optional, intent(in) :: level type(r_solver_operator_type), pointer :: op - op => eliminated_q2t + + if ( present(level) ) then + op => eliminated_q2t(level) + else + op => eliminated_q2t(1) + end if end function get_eliminated_q2t !> @brief Function to return a pointer to the q3t operator. diff --git a/science/gungho/source/field/field_indices_mod.F90 b/science/gungho/source/field/field_indices_mod.F90 index 96b00987..47714d79 100644 --- a/science/gungho/source/field/field_indices_mod.F90 +++ b/science/gungho/source/field/field_indices_mod.F90 @@ -33,9 +33,9 @@ module field_indices_mod integer, parameter :: isw_q = 4 ! vorticity ! For Semi-Implicit Solver - integer, parameter :: isol_u = 1 ! wind + integer, parameter :: isol_u = 2 ! wind integer, parameter :: isol_uv = isol_u ! uv wind if split - must be the same as isol_u - integer, parameter :: isol_p = 2 ! Exner pressure + integer, parameter :: isol_p = 1 ! Exner pressure integer, parameter :: isol_w = 3 ! w wind if split contains diff --git a/science/gungho/source/psy/psykal_lite_mod.F90 b/science/gungho/source/psy/psykal_lite_mod.F90 index 80509dc1..93be8dd3 100644 --- a/science/gungho/source/psy/psykal_lite_mod.F90 +++ b/science/gungho/source/psy/psykal_lite_mod.F90 @@ -39,6 +39,26 @@ module psykal_lite_mod !> Non pointwise Kernels + subroutine print_dofs(field, msg, ndofs) + + type(r_solver_field_type), intent(in) :: field + type(r_solver_field_proxy_type) :: field_proxy + character(len=*), intent(in) :: msg + + integer(kind=i_def), intent(in) :: ndofs + + integer(kind=i_def) :: df, ndf + + field_proxy = field%get_proxy() + + ndf = field_proxy%vspace%get_last_dof_owned() + + do df = 1, ndofs + write(6,*) trim(msg),df,ndf,field_proxy%data(df) + end do + + end subroutine print_dofs + !------------------------------------------------------------------------------- subroutine invoke_compute_dof_level_kernel(level)