petsc_solver_mod Module Reference

PETSc solver component to call out to PETSc for solving the Poisson equation for pressure. More...

Functions/Subroutines
type(component_descriptor_type) function, public	petsc_solver_get_descriptor ()
	Provides the descriptor back to the caller and is used in component registration. More...
subroutine	finalisation_callback (current_state)
subroutine	init_callback (current_state)
	Called upon model initialisation. Will basically read from the options database and set options in the database that are appropriate. More...
subroutine	ksp_monitor (ksp, n, rnorm, dummy, ierr)
	Monitor procedure called on each iteration, this is provided only at debugging level. More...
subroutine	timestep_callback (current_state)
	Timestep hook which is called at each timestep to solve the pressure equation via PETSc. This will call to set and precondition the matrix on the first call (first timestep) only, and will set the RHS and X initial guess on each call as one would expect. At the end a halo swap is performed on p, which is needed for pstep. More...
character(len=string_length) function	determine_convegence_reason (r_code)
subroutine	compute_initial_guess (ksp, x, dummy, ierr)
	Determines the initial guess, this is called from within PETSc and sets it to be the last value of P or zero if it is the first timestep. More...
subroutine	compute_rhs (ksp, b, dummy, ierr)
	Callback issued from the PETSc library to compute the RHS, this is called every timestep (as we have a different RHS) More...
subroutine	copy_data_to_petsc_pointer (x, zs, ys, xs, zm, ym, xm, src_values)
	Copies data into a data pointer which is provided by PETSc, doing it this ways allows us to give a shape to the pointer data and hence the data copy is simpler code wise. More...
subroutine	copy_petsc_pointer_to_data (x, z_start_idx, z_end_idx, y_start_idx, y_end_idx, x_start_idx, x_end_idx, target_values)
	Copies the data in a pointer that was provided by PETSc into some target data, doing it this way means we can give a shape to the pointer. More...
subroutine	compute_matrix (ksp, A, B, dummy, ierr)
	Call back issued from within PETSc to create the matrix, this is only called once (the first PETSc run) More...
subroutine	apply_dimension_bounds (dim_size, dim_processes, length_array)
	Applies dimension bounds to determine the number of elements held locally for each process in that dimension. More...
subroutine	complete_psrce_calculation (current_state, y_halo_size, x_halo_size)
	Completes the psrce calculation by waiting on all outstanding psrce communications to complete and then combine the received values with the P field for U and V. More...
subroutine	copy_p_to_halo_buffer (current_state, neighbour_description, dim, source_index, pid_location, current_page, source_data)
	Copies the p field data to halo buffers for a specific process in a dimension and halo cell. More...
subroutine	copy_halo_buffer_to_p (current_state, neighbour_description, dim, target_index, neighbour_location, current_page, source_data)
subroutine	perform_local_data_copy_for_p (current_state, halo_depth, involve_corners, source_data)
	Does local data copying for P variable halo swap. More...

Variables
integer	z_start
integer	z_end
integer	y_start
integer	y_end
integer	x_start
integer	x_end
type(halo_communication_type), save	halo_swap_state
real(kind=default_precision), dimension(:,:,:), allocatable	p_source
real(kind=default_precision), dimension(:,:,:), allocatable	prev_p
real(kind=default_precision)	cx
real(kind=default_precision)	cy
real(kind=default_precision), dimension(:), allocatable	rdzn
real(kind=default_precision), dimension(:), allocatable	rdz
real(kind=default_precision), dimension(:), allocatable	rho
real(kind=default_precision), dimension(:), allocatable	rhon
real(kind=default_precision), dimension(:), allocatable	dz
real(kind=default_precision), dimension(:), allocatable	dzn

Detailed Description

PETSc solver component to call out to PETSc for solving the Poisson equation for pressure.

Dummy stub when not compiling with PETSc iterative solver.

Function/Subroutine Documentation

apply_dimension_bounds()

subroutine petsc_solver_mod::apply_dimension_bounds ( integer, intent(in)  dim_size, integer, intent(in)  dim_processes, integer, dimension(:), intent(out)  length_array  )

private

Applies dimension bounds to determine the number of elements held locally for each process in that dimension.

Parameters

dim_size	The global size of the dimension
dim_processes	The number of processes in that dimension
length_array	Output array of size dim_processes which contains the number of elements held locally for each process

Definition at line 418 of file petsc_solver.F90.

    integer, intent(in) :: dim_size, dim_processes
    integer, intent(out) :: length_array(:)
 
    integer :: i, dimension_division, dimension_extra
 
    dimension_division = dim_size / dim_processes
    length_array= dimension_division    
    
    dimension_extra = dim_size - (dimension_division * dim_processes)
    if (dimension_extra .gt. 0) then
      do i=1, dimension_extra
        length_array(i)=length_array(i)+1
      end do
    end if

Here is the caller graph for this function:

complete_psrce_calculation()

subroutine petsc_solver_mod::complete_psrce_calculation ( type(model_state_type), intent(inout), target  current_state, integer, intent(in)  y_halo_size, integer, intent(in)  x_halo_size  )

private

Completes the psrce calculation by waiting on all outstanding psrce communications to complete and then combine the received values with the P field for U and V.

Parameters

current_state	The current model state
y_halo_size	The halo size in the Y dimension
x_halo_size	The halo size in the X dimension

Definition at line 440 of file petsc_solver.F90.

    type(model_state_type), target, intent(inout) :: current_state
    integer, intent(in) :: y_halo_size, x_halo_size
 
    integer :: ierr, combined_handles(2), i, j, k
 
    combined_handles(1)=current_state%psrce_x_hs_recv_request
    combined_handles(2)=current_state%psrce_y_hs_recv_request
    call mpi_waitall(2, combined_handles, mpi_statuses_ignore, ierr)
 
    do j=current_state%local_grid%local_domain_start_index(y_index), current_state%local_grid%local_domain_end_index(y_index)
      do k=2,current_state%local_grid%size(z_index)
#ifdef U_ACTIVE
        current_state%p%data(k,j,x_halo_size+1)=current_state%p%data(k,j,x_halo_size+1)-&
               current_state%psrce_recv_buffer_x(k-1,j-x_halo_size)
#endif
#ifdef V_ACTIVE
        if (j .gt. y_halo_size+1) current_state%p%data(k, j, x_halo_size+1)=current_state%p%data(k, j, x_halo_size+1)-&
             current_state%global_grid%configuration%horizontal%cy * current_state%sv%data(k, j-1, x_halo_size+1)
#endif
      end do
    end do
      
#ifdef V_ACTIVE
    do i=current_state%local_grid%local_domain_start_index(x_index), current_state%local_grid%local_domain_end_index(x_index)
      do k=2,current_state%local_grid%size(z_index)
        current_state%p%data(k,y_halo_size+1,i)=current_state%p%data(k,y_halo_size+1,i)-&
             current_state%psrce_recv_buffer_y(k-1,i-y_halo_size)
      end do
    end do
#endif
 
    combined_handles(1)=current_state%psrce_x_hs_send_request
    combined_handles(2)=current_state%psrce_y_hs_send_request
    call mpi_waitall(2, combined_handles, mpi_statuses_ignore, ierr)

Here is the caller graph for this function:

compute_initial_guess()

subroutine petsc_solver_mod::compute_initial_guess (   ksp,   x, integer, dimension(*)  dummy,   ierr  )

private

Determines the initial guess, this is called from within PETSc and sets it to be the last value of P or zero if it is the first timestep.

Parameters

ksp	The PETSc KSP data structure which contains the context of the solver
X	The X vector to create an initial guess for
dummy	Dummy argument to maintain consistency with the PETSc C interface
ierr	PETSc error code

Definition at line 234 of file petsc_solver.F90.

    petscerrorcode :: ierr
    ksp :: ksp
    vec :: x
    integer :: dummy(*)
 
    dm dm
    petscscalar, pointer :: xx(:)
    petscint :: zs, ys, xs, zm, ym, xm
 
    call kspgetdm(ksp, dm, ierr)
    call dmdagetcorners(dm, zs, ys, xs, zm, ym, xm, ierr)
    call vecgetarrayf90(x, xx, ierr)
    call copy_data_to_petsc_pointer(xx, zs, ys, xs, zm, ym, xm, prev_p)
    call vecrestorearrayf90(x, xx, ierr)

Here is the call graph for this function:

Here is the caller graph for this function:

compute_matrix()

subroutine petsc_solver_mod::compute_matrix (   ksp,   A,   B, integer, dimension(*)  dummy,   ierr  )

private

Call back issued from within PETSc to create the matrix, this is only called once (the first PETSc run)

Parameters

ksp	The solver data structure which provides context for this run
A	Matrix A this is the linear operator
B	Matrix B this is the preconditioning matrix
dummy	Dummy argument to ensure consistency with C PETSc interface
ierr	The PETSc error code

Definition at line 323 of file petsc_solver.F90.

    petscerrorcode :: ierr
    ksp :: ksp
    mat :: a, b
    integer :: dummy(*)
 
    dm dm
    petscint :: zs, ys,xs, zm, ym, xm, mz, my, mx
    real(kind=default_precision) ::  v(7), hx, hy
    real(kind=default_precision), dimension(:), allocatable ::  hzu, hzd, d_sc  
    matstencil :: row(4),col(4,7)
    integer :: i, j, k, num
 
    call kspgetdm(ksp, dm, ierr)
    call dmdagetinfo(dm,petsc_null_integer, mz, my, mx, petsc_null_integer, petsc_null_integer, petsc_null_integer, &
     petsc_null_integer, petsc_null_integer, petsc_null_integer, petsc_null_integer, petsc_null_integer, petsc_null_integer, ierr)
    call dmdagetcorners(dm, zs, ys, xs, zm, ym, xm, ierr)
 
    allocate(hzu(0:mz-1), hzd(0:mz-1), d_sc(0:mz-1))
 
    do k=0, mz-1
      d_sc(k)=rdz(k+1) / rhon(k+1)
      if (k .gt. 0) then
        hzd(k)=rho(k) * rdzn(k+1)
      else
        hzd(k)=0.0_default_precision
      end if
      if (k .lt. mz-1) then
        hzu(k)=rho(k+1) * rdzn(k+2)
      else
        hzu(k)=0.0_default_precision
      end if
    end do   
 
    hx = cx * cx
    hy = cy * cy  
    do i=xs, xs+xm-1
      do j=ys,ys+ym-1
        do k=zs, zs+zm-1
          row(matstencil_i) = k
          row(matstencil_j) = j
          row(matstencil_k) = i
 
          v(1)=hy
          col(matstencil_i, 1)=k
          col(matstencil_j, 1)=j-1
          col(matstencil_k, 1)=i
          v(2)=hx
          col(matstencil_i, 2)=k
          col(matstencil_j, 2)=j
          col(matstencil_k, 2)=i-1
          v(3)=d_sc(k) * (-(hzu(k) + hzd(k))) - (2.0*(hx + hy))
          col(matstencil_i, 3)=k
          col(matstencil_j, 3)=j
          col(matstencil_k, 3)=i
          v(4)=hx
          col(matstencil_i, 4)=k
          col(matstencil_j, 4)=j 
          col(matstencil_k, 4)=i+1
          v(5)=hy
          col(matstencil_i, 5)=k
          col(matstencil_j, 5)=j+1
          col(matstencil_k, 5)=i
          num=6
          if (k .gt. zs) then
            v(num)=(d_sc(k)*hzd(k))
            col(matstencil_i, num)=k-1
            col(matstencil_j, num)=j 
            col(matstencil_k, num)=i
            num=num+1
          end if
          if (k .lt. mz-1) then
            v(num)=(d_sc(k)*hzu(k))
            col(matstencil_i, num)=k+1
            col(matstencil_j, num)=j
            col(matstencil_k, num)=i
            num=num+1
          end if
          call matsetvaluesstencil(b, 1, row, num-1, col, v, insert_values, ierr)
        end do
      end do
    end do
    call matassemblybegin(b, mat_final_assembly, ierr)
    call matassemblyend(b, mat_final_assembly, ierr)
    if ( a .ne. b) then
      call matassemblybegin(a, mat_final_assembly, ierr)
      call matassemblyend(a, mat_final_assembly, ierr)
    endif
    deallocate(hzu, hzd, d_sc)

Here is the caller graph for this function:

compute_rhs()

subroutine petsc_solver_mod::compute_rhs (   ksp,   b, integer, dimension(*)  dummy,   ierr  )

private

Callback issued from the PETSc library to compute the RHS, this is called every timestep (as we have a different RHS)

Parameters

ksp	The KSP data structure which contains the solver context
b	The RHS
dummy	Dummy argument for compatability with the C PETSc interface
ierr	Error code

Definition at line 256 of file petsc_solver.F90.

    petscerrorcode :: ierr
    ksp :: ksp
    vec :: b
    integer :: dummy(*)
    
    dm dm
    petscscalar, pointer :: xx(:)
    petscint :: zs, ys, xs, zm, ym, xm
 
    call kspgetdm(ksp, dm, ierr)
    call dmdagetcorners(dm, zs, ys, xs, zm, ym, xm, ierr)
    call vecgetarrayf90(b, xx, ierr)
    call copy_data_to_petsc_pointer(xx, zs, ys, xs, zm, ym, xm, p_source)
    call vecrestorearrayf90(b, xx, ierr)

Here is the call graph for this function:

Here is the caller graph for this function:

copy_data_to_petsc_pointer()

subroutine petsc_solver_mod::copy_data_to_petsc_pointer ( dimension(zs:zs+zm-1, ys:ys+ym-1, xs:xs+xm-1), intent(inout)  x, integer, intent(in)  zs, integer, intent(in)  ys, integer, intent(in)  xs, integer, intent(in)  zm, integer, intent(in)  ym, integer, intent(in)  xm, real(kind=default_precision), dimension(:,:,:), intent(in)  src_values  )

private

Copies data into a data pointer which is provided by PETSc, doing it this ways allows us to give a shape to the pointer data and hence the data copy is simpler code wise.

Parameters

x	The target pointer to copy into
zs	Local starting point in Z
ys	Local starting point in Y
xs	Local starting point in X
zm	Number of points held locally in Z
ym	Number of points held locally in Y
xm	Number of points held locally in X
src_values	The data values to copy from into this target data array

Definition at line 283 of file petsc_solver.F90.

    integer, intent(in) :: zs, ys, xs, zm, ym, xm
    petscscalar, intent(inout) :: x(zs:zs+zm-1, ys:ys+ym-1, xs:xs+xm-1)
    real(kind=default_precision), dimension(:,:,:), intent(in) :: src_values
 
    integer :: i, j, k
 
    do i=xs, xs+xm-1
      do j=ys, ys+ym-1
        do k=zs, zs+zm-1
          x(k, j, i)=src_values((k-zs)+1, (j-ys)+1, (i-xs)+1)
        end do
      end do
    end do

Here is the caller graph for this function:

copy_halo_buffer_to_p()

subroutine petsc_solver_mod::copy_halo_buffer_to_p ( type(model_state_type), intent(inout)  current_state, type(neighbour_description_type), intent(inout)  neighbour_description, integer, intent(in)  dim, integer, intent(in)  target_index, integer, intent(in)  neighbour_location, integer, dimension(:), intent(inout)  current_page, type(field_data_wrapper_type), dimension(:), intent(in), optional  source_data  )

private

Definition at line 504 of file petsc_solver.F90.

    type(model_state_type), intent(inout) :: current_state
    integer, intent(in) :: dim, target_index, neighbour_location
    integer, intent(inout) :: current_page(:)
    type(neighbour_description_type), intent(inout) :: neighbour_description
    type(field_data_wrapper_type), dimension(:), intent(in), optional :: source_data
 
    call copy_buffer_to_field(current_state%local_grid, neighbour_description%recv_halo_buffer, current_state%p%data, &
         dim, target_index, current_page(neighbour_location))
 
    current_page(neighbour_location)=current_page(neighbour_location)+1

Here is the caller graph for this function:

copy_p_to_halo_buffer()

subroutine petsc_solver_mod::copy_p_to_halo_buffer ( type(model_state_type), intent(inout)  current_state, type(neighbour_description_type), intent(inout)  neighbour_description, integer, intent(in)  dim, integer, intent(in)  source_index, integer, intent(in)  pid_location, integer, dimension(:), intent(inout)  current_page, type(field_data_wrapper_type), dimension(:), intent(in), optional  source_data  )

private

Copies the p field data to halo buffers for a specific process in a dimension and halo cell.

Parameters

current_state	The current model state
neighbour_descriptions	Description of the neighbour halo swapping status
dim	Dimension to copy from
source_index	The source index of the dimension we are reading from in the prognostic field
pid_location	Location of the neighbouring process in the local stored data structures
current_page	The current (next) buffer page to copy into
source_data	Optional source data which is read from

Definition at line 485 of file petsc_solver.F90.

    type(model_state_type), intent(inout) :: current_state
    integer, intent(in) :: dim, pid_location, source_index
    integer, intent(inout) :: current_page(:)
    type(neighbour_description_type), intent(inout) :: neighbour_description
    type(field_data_wrapper_type), dimension(:), intent(in), optional :: source_data
 
    call copy_field_to_buffer(current_state%local_grid, neighbour_description%send_halo_buffer, current_state%p%data, &
         dim, source_index, current_page(pid_location))
 
    current_page(pid_location)=current_page(pid_location)+1

Here is the caller graph for this function:

copy_petsc_pointer_to_data()

subroutine petsc_solver_mod::copy_petsc_pointer_to_data (   x, integer, intent(in)  z_start_idx, integer, intent(in)  z_end_idx, integer, intent(in)  y_start_idx, integer, intent(in)  y_end_idx, integer, intent(in)  x_start_idx, integer, intent(in)  x_end_idx, real(kind=default_precision), dimension(:,:,:), intent(inout)  target_values  )

private

Copies the data in a pointer that was provided by PETSc into some target data, doing it this way means we can give a shape to the pointer.

Parameters

x	The pointer provided by PETSc that we copy from
z_start	Starting point for the target array in Z
z_end	Ending point for the target array in Z
y_start	Starting point for the target array in Y
y_end	Ending point for the target array in Y
x_start	Starting point for the target array in X
x_end	Ending point for the target array in X
target_values	Target array to copy into

Definition at line 309 of file petsc_solver.F90.

     integer, intent(in) :: z_start_idx, z_end_idx, y_start_idx, y_end_idx, x_start_idx, x_end_idx
     petscscalar :: x((z_end_idx-z_start_idx)+1, (y_end_idx-y_start_idx)+1, (x_end_idx-x_start_idx)+1)
     real(kind=default_precision), dimension(:,:,:), intent(inout) :: target_values
 
     target_values(z_start_idx:z_end_idx, y_start_idx:y_end_idx, x_start_idx:x_end_idx)=x

Here is the caller graph for this function:

determine_convegence_reason()

character(len=string_length) function petsc_solver_mod::determine_convegence_reason ( integer, intent(in)  r_code )

private

Definition at line 208 of file petsc_solver.F90.

    integer, intent(in) :: r_code
 
    if (r_code == 1 .or. r_code == 2) then
      determine_convegence_reason="relative tolerance of residual norm"
    else if (r_code == 9 .or. r_code == 3) then
      determine_convegence_reason="residual norm less than absolute tolerance"
    else if (r_code == 4) then
      determine_convegence_reason="number of iterations exceeded maximum"
    else if (r_code == -3) then
      determine_convegence_reason="diverged as number of iterations exceeded maximum before any convergence criteria"
    else if (r_code == -4) then
      determine_convegence_reason="diverged as divergence tolerance exceeded"
    else if (r_code .gt. 0) then
      determine_convegence_reason="convergence code of "//conv_to_string(r_code)
    else if (r_code .lt. 0) then
      determine_convegence_reason="divergence code of "//conv_to_string(r_code)
    end if    

Here is the caller graph for this function:

finalisation_callback()

subroutine petsc_solver_mod::finalisation_callback ( type(model_state_type), intent(inout), target  current_state )

private

Definition at line 51 of file petsc_solver.F90.

    type(model_state_type), target, intent(inout) :: current_state
 
    petscerrorcode :: ierr
 
    call petscfinalize(ierr)

Here is the caller graph for this function:

init_callback()

subroutine petsc_solver_mod::init_callback ( type(model_state_type), intent(inout), target  current_state )

private

Called upon model initialisation. Will basically read from the options database and set options in the database that are appropriate.

Parameters

current_state

The current model state_mod

Definition at line 62 of file petsc_solver.F90.

    type(model_state_type), target, intent(inout) :: current_state
 
    petscerrorcode :: ierr
    ksptype :: ksp_type
    pc :: pc_instance
    pctype :: pc_type
    integer :: y_lengths(current_state%parallel%dim_sizes(Y_INDEX)), x_lengths(current_state%parallel%dim_sizes(X_INDEX))
    integer new_communicator, my_transposed_rank, x_rank_val, y_rank_val
 
    call apply_dimension_bounds(current_state%global_grid%size(y_index), current_state%parallel%dim_sizes(y_index), y_lengths)
    call apply_dimension_bounds(current_state%global_grid%size(x_index), current_state%parallel%dim_sizes(x_index), x_lengths)
 
    allocate(p_source(current_state%local_grid%size(z_index)-1, current_state%local_grid%size(y_index), &
         current_state%local_grid%size(x_index)), prev_p(current_state%local_grid%size(z_index)-1, &
         current_state%local_grid%size(y_index), current_state%local_grid%size(x_index)))
 
    cx=current_state%global_grid%configuration%horizontal%cx
    cy=current_state%global_grid%configuration%horizontal%cy
    allocate(rdzn(current_state%local_grid%size(z_index)-1), rdz(current_state%local_grid%size(z_index)-1), &
         rho(current_state%local_grid%size(z_index)-1), rhon(current_state%local_grid%size(z_index)-1), &
         dz(current_state%local_grid%size(z_index)-1), dzn(current_state%local_grid%size(z_index)-1))
    rdzn=current_state%global_grid%configuration%vertical%rdzn(2:)
    rdz=current_state%global_grid%configuration%vertical%rdz(2:)
    rho=current_state%global_grid%configuration%vertical%rho(2:)
    rhon=current_state%global_grid%configuration%vertical%rhon(2:)
    dz=current_state%global_grid%configuration%vertical%dz(2:)
    dzn=current_state%global_grid%configuration%vertical%dzn(2:)
 
    z_start=current_state%local_grid%halo_size(z_index)+2
    z_end=current_state%local_grid%halo_size(z_index)+current_state%local_grid%size(z_index)
    y_start=current_state%local_grid%halo_size(y_index)+1
    y_end=current_state%local_grid%halo_size(y_index)+current_state%local_grid%size(y_index)
    x_start=current_state%local_grid%halo_size(x_index)+1
    x_end=current_state%local_grid%halo_size(x_index)+current_state%local_grid%size(x_index)
 
    y_rank_val = current_state%parallel%my_rank / current_state%parallel%dim_sizes(x_index)
    x_rank_val = mod(current_state%parallel%my_rank, current_state%parallel%dim_sizes(x_index))    
    my_transposed_rank = x_rank_val*current_state%parallel%dim_sizes(y_index) + y_rank_val
 
    call mpi_comm_split(current_state%parallel%monc_communicator, 1, my_transposed_rank, new_communicator, ierr)
    petsc_comm_world = new_communicator
 
    call petscinitialize(petsc_null_character, ierr)    
    call kspcreate(petsc_comm_world, ksp, ierr)
    call dmdacreate3d(petsc_comm_world, dm_boundary_none, dm_boundary_periodic, dm_boundary_periodic, dmda_stencil_star, &
         current_state%global_grid%size(z_index)-1, current_state%global_grid%size(y_index), &
         current_state%global_grid%size(x_index), 1, current_state%parallel%dim_sizes(y_index), &
         current_state%parallel%dim_sizes(x_index), 1, 1, (/ current_state%global_grid%size(z_index)-1 /), &
         y_lengths, x_lengths, da, ierr)
    call kspsetdm(ksp, da, ierr)
    call kspsetcomputerhs(ksp, compute_rhs, petsc_null_object, ierr)
    call kspsetcomputeoperators(ksp, compute_matrix, petsc_null_object, ierr)
    call kspsetcomputeinitialguess(ksp, compute_initial_guess, petsc_null_object, ierr)
    call petscoptionssetvalue("-options_left", "no", ierr)
    if (trim(options_get_string(current_state%options_database, "solver_type")) .ne. "auto") then
      call petscoptionssetvalue("-ksp_type", trim(options_get_string(current_state%options_database, "solver_type")), ierr)
    end if
    if (trim(options_get_string(current_state%options_database, "preconditioner_type")) .ne. "auto") then
      call petscoptionssetvalue("-pc_type", trim(options_get_string(current_state%options_database, "preconditioner_type")), ierr)
    end if
    if (trim(options_get_string(current_state%options_database, "norm_type")) .ne. "auto") then
      if (trim(options_get_string(current_state%options_database, "norm_type")) .eq. "preconditioned" .or. &
           trim(options_get_string(current_state%options_database, "norm_type")) .eq. "unpreconditioned" .or. &
           trim(options_get_string(current_state%options_database, "norm_type")) .eq. "natural" .or. &
           trim(options_get_string(current_state%options_database, "norm_type")) .eq. "none") then
        call petscoptionssetvalue("-ksp_norm_type", trim(options_get_string(current_state%options_database, "norm_type")), ierr)
      else
        call log_master_log(log_error, "Configured PETSc norm type of '"//&
             trim(options_get_string(current_state%options_database, "norm_type"))//"' not recognised")
      end if      
    end if    
    call petscoptionssetvalue("-ksp_rtol", conv_to_string(options_get_real(current_state%options_database, "tolerance")), ierr)
    if (options_get_integer(current_state%options_database, "max_iterations") .gt. 0) then
      call petscoptionssetvalue("-ksp_max_it", &
           conv_to_string(options_get_integer(current_state%options_database, "max_iterations")), ierr)
    end if
    call kspsetfromoptions(ksp, ierr)
    call kspsetinitialguessnonzero(ksp, petsc_true, ierr)
 
    if (log_is_master() .and. log_get_logging_level() == log_debug) then
      call kspmonitorset(ksp,ksp_monitor,petsc_null_object, petsc_null_function,ierr)
    end if    
    prev_p=0.0_default_precision
    call init_halo_communication(current_state, get_single_field_per_halo_cell, halo_swap_state, 1, .false.)
    call kspgettype(ksp, ksp_type, ierr)
    call kspgetpc(ksp, pc_instance, ierr)
    call pcgettype(pc_instance, pc_type, ierr)
    call log_master_log(log_info, "PETSc iterative solver initialised, using "//trim(pc_type)//" preconditioner with "//&
         trim(ksp_type)//" solver")

Here is the call graph for this function:

Here is the caller graph for this function:

ksp_monitor()

subroutine petsc_solver_mod::ksp_monitor (   ksp,   n,   rnorm,   dummy,   ierr  )

private

Monitor procedure called on each iteration, this is provided only at debugging level.

Parameters

ksp	The PETSc KSP data structure which contains the context of the solver
n	The iteration
rnorm	The relative norm
dummy	Dummy argument to maintain consistency with the PETSc C interface
ierr	PETSc error code

Definition at line 160 of file petsc_solver.F90.

    ksp              ksp
    petscerrorcode ierr
    petscint n,dummy
    petscreal rnorm
 
    call log_log(log_debug, "Iteration="//trim(conv_to_string(n))//" Rnorm="//trim(conv_to_string(rnorm, &
         exponent_small_numbers=.true.)))

Here is the caller graph for this function:

perform_local_data_copy_for_p()

subroutine petsc_solver_mod::perform_local_data_copy_for_p ( type(model_state_type), intent(inout)  current_state, integer, intent(in)  halo_depth, logical, intent(in)  involve_corners, type(field_data_wrapper_type), dimension(:), intent(in), optional  source_data  )

private

Does local data copying for P variable halo swap.

Parameters

current_state	The current model state_mod
source_data	Optional source data which is written into

Definition at line 521 of file petsc_solver.F90.

    type(model_state_type), intent(inout) :: current_state
    integer, intent(in) :: halo_depth
    logical, intent(in) :: involve_corners
    type(field_data_wrapper_type), dimension(:), intent(in), optional :: source_data
 
    call perform_local_data_copy_for_field(current_state%p%data, current_state%local_grid, &
         current_state%parallel%my_rank, halo_depth, involve_corners)

Here is the caller graph for this function:

petsc_solver_get_descriptor()

type(component_descriptor_type) function, public petsc_solver_mod::petsc_solver_get_descriptor

Provides the descriptor back to the caller and is used in component registration.

Returns

The PETSc solvre component descriptor

The termination check component descriptor

Definition at line 43 of file petsc_solver.F90.

    petsc_solver_get_descriptor%name="petsc_solver"
    petsc_solver_get_descriptor%version=0.1
    petsc_solver_get_descriptor%initialisation=>init_callback
    petsc_solver_get_descriptor%timestep=>timestep_callback
    petsc_solver_get_descriptor%finalisation=>finalisation_callback

Here is the call graph for this function:

timestep_callback()

subroutine petsc_solver_mod::timestep_callback ( type(model_state_type), intent(inout), target  current_state )

private

Timestep hook which is called at each timestep to solve the pressure equation via PETSc. This will call to set and precondition the matrix on the first call (first timestep) only, and will set the RHS and X initial guess on each call as one would expect. At the end a halo swap is performed on p, which is needed for pstep.

Parameters

current_state

The current model state_mod

Definition at line 174 of file petsc_solver.F90.

    type(model_state_type), target, intent(inout) :: current_state
    
    petscerrorcode :: ierr
    vec x
    petscscalar, pointer :: xx(:)
    integer :: its, i, j, k
    petscreal :: rnorm
    kspconvergedreason :: reason
 
    call complete_psrce_calculation(current_state, current_state%local_grid%halo_size(y_index), &
         current_state%local_grid%halo_size(x_index))    
 
    p_source=current_state%p%data(z_start:z_end, y_start:y_end, x_start:x_end)
    call kspsolve(ksp, petsc_null_object, petsc_null_object, ierr)
    call kspgetsolution(ksp, x, ierr)
    call vecgetarrayreadf90(x, xx, ierr)
    
    call copy_petsc_pointer_to_data(xx, z_start, z_end, y_start, y_end, x_start, x_end, current_state%p%data)
    call vecrestorearrayreadf90(x, xx, ierr)
    
    prev_p=current_state%p%data(z_start:z_end, y_start:y_end, x_start:x_end)
 
    if (log_is_master() .and. log_get_logging_level() == log_debug) then
      call kspgetiterationnumber(ksp, its, ierr)
      call kspgetresidualnorm(ksp, rnorm, ierr)
      call kspgetconvergedreason(ksp, reason, ierr)
      call log_log(log_debug, "Converged in "//trim(conv_to_string(its))//" rnorm="//trim(conv_to_string(&
           rnorm, exponent_small_numbers=.true.))// " ("//trim(determine_convegence_reason(reason))//")")
    end if
    call blocking_halo_swap(current_state, halo_swap_state, copy_p_to_halo_buffer, &
         perform_local_data_copy_for_p, copy_halo_buffer_to_p)

Here is the call graph for this function:

Here is the caller graph for this function:

Variable Documentation

cx

real(kind=default_precision) petsc_solver_mod::cx

private

Definition at line 34 of file petsc_solver.F90.

  real(kind=default_precision) :: cx, cy

cy

real(kind=default_precision) petsc_solver_mod::cy

private

Definition at line 34 of file petsc_solver.F90.

dz

real(kind=default_precision), dimension(:), allocatable petsc_solver_mod::dz

private

Definition at line 35 of file petsc_solver.F90.

dzn

real(kind=default_precision), dimension(:), allocatable petsc_solver_mod::dzn

private

Definition at line 35 of file petsc_solver.F90.

halo_swap_state

type(halo_communication_type), save petsc_solver_mod::halo_swap_state

private

Definition at line 31 of file petsc_solver.F90.

  type(halo_communication_type), save :: halo_swap_state

p_source

real(kind=default_precision), dimension(:,:,:), allocatable petsc_solver_mod::p_source

private

Definition at line 33 of file petsc_solver.F90.

  real(kind=default_precision), dimension(:,:,:), allocatable :: p_source, prev_p

prev_p

real(kind=default_precision), dimension(:,:,:), allocatable petsc_solver_mod::prev_p

private

Definition at line 33 of file petsc_solver.F90.

rdz

real(kind=default_precision), dimension(:), allocatable petsc_solver_mod::rdz

private

Definition at line 35 of file petsc_solver.F90.

rdzn

real(kind=default_precision), dimension(:), allocatable petsc_solver_mod::rdzn

private

Definition at line 35 of file petsc_solver.F90.

  real(kind=default_precision), dimension(:), allocatable :: rdzn, rdz, rho, rhon, dz, dzn

rho

real(kind=default_precision), dimension(:), allocatable petsc_solver_mod::rho

private

Definition at line 35 of file petsc_solver.F90.

rhon

real(kind=default_precision), dimension(:), allocatable petsc_solver_mod::rhon

private

Definition at line 35 of file petsc_solver.F90.

x_end

integer petsc_solver_mod::x_end

private

Definition at line 30 of file petsc_solver.F90.

x_start

integer petsc_solver_mod::x_start

private

Definition at line 30 of file petsc_solver.F90.

y_end

integer petsc_solver_mod::y_end

private

Definition at line 30 of file petsc_solver.F90.

y_start

integer petsc_solver_mod::y_start

private

Definition at line 30 of file petsc_solver.F90.

z_end

integer petsc_solver_mod::z_end

private

Definition at line 30 of file petsc_solver.F90.

z_start

integer petsc_solver_mod::z_start

private

Definition at line 30 of file petsc_solver.F90.

  integer :: z_start, z_end, y_start, y_end, x_start, x_end