Skip to content

Mixed multigrid #197

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
thomasmelvin wants to merge 2 commits into
base: main
Choose a base branch
from
Open

Mixed multigrid #197

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
c72b43a
merge of code from old branch
thomasmelvin Jan 23, 2026
1fa5870
add in core branch
thomasmelvin Jan 26, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 6 additions & 5 deletions applications/gungho_model/example/configuration.nml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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.,
/
Expand All @@ -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.,
Expand Down
5 changes: 1 addition & 4 deletions build/local_build.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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}",
Expand Down
4 changes: 2 additions & 2 deletions dependencies.yaml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
12 changes: 12 additions & 0 deletions rose-stem/app/gungho_model/opt/rose-app-mixed-mg.conf
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
33 changes: 33 additions & 0 deletions rose-stem/site/common/gungho_model/tasks_gungho_model.cylc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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" %}

Expand Down
8 changes: 8 additions & 0 deletions rose-stem/site/meto/groups/groups_gungho_model.cylc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand Down
4 changes: 2 additions & 2 deletions science/gungho/rose-meta/lfric-gungho/HEAD/rose-meta.conf
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
101 changes: 101 additions & 0 deletions science/gungho/source/algorithm/core_dynamics/intermesh_mappings_alg_mod.x90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
95 changes: 95 additions & 0 deletions science/gungho/source/algorithm/solver/create_pressure_preconditioner_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
Loading
Loading