Neutral Diffusion: Add option for discontinuous reconstruction of T/S

src/ALE/MOM_remapping.F90

@@ -37,8 +37,8 @@ module MOM_remapping
 
 ! The following routines are visible to the outside world
 public remapping_core_h, remapping_core_w
-public initialize_remapping, end_remapping, remapping_set_param
-public remapping_unit_tests
+public initialize_remapping, end_remapping, remapping_set_param, extract_member_remapping_CS
+public remapping_unit_tests, build_reconstructions_1d, average_value_ppoly
 public dzFromH1H2
 
 ! The following are private parameter constants
@@ -109,6 +109,26 @@ subroutine remapping_set_param(CS, remapping_scheme, boundary_extrapolation,  &
   endif
 end subroutine remapping_set_param
 
+subroutine extract_member_remapping_CS(CS, remapping_scheme, degree, boundary_extrapolation, check_reconstruction, &
+                                       check_remapping, force_bounds_in_subcell)
+  type(remapping_CS), intent(in) :: CS !< Control structure for remapping module
+  integer, optional, intent(out) :: remapping_scheme        !< Determines which reconstruction scheme to use
+  integer, optional, intent(out) :: degree                  !< Degree of polynomial reconstruction
+  logical, optional, intent(out) :: boundary_extrapolation  !< If true, extrapolate boundaries
+  logical, optional, intent(out) :: check_reconstruction    !< If true, reconstructions are checked for consistency.
+  logical, optional, intent(out) :: check_remapping         !< If true, the result of remapping are checked
+                                                            !!  for conservation and bounds.
+  logical, optional, intent(out) :: force_bounds_in_subcell !< If true, the intermediate values used in
+                                                            !! remapping are forced to be bounded.
+
+  if (present(remapping_scheme)) remapping_scheme = CS%remapping_scheme
+  if (present(degree)) degree = CS%degree
+  if (present(boundary_extrapolation)) boundary_extrapolation = CS%boundary_extrapolation
+  if (present(check_reconstruction)) check_reconstruction = CS%check_reconstruction
+  if (present(check_remapping)) check_remapping = CS%check_remapping
+  if (present(force_bounds_in_subcell)) force_bounds_in_subcell = CS%force_bounds_in_subcell
+
+end subroutine extract_member_remapping_CS
 !> Calculate edge coordinate x from cell width h
 subroutine buildGridFromH(nz, h, x)
   integer,               intent(in)    :: nz !< Number of cells
@@ -175,8 +195,7 @@ subroutine remapping_core_h(CS,  n0, h0, u0, n1, h1, u1)
   integer :: k
   real :: eps, h0tot, h0err, h1tot, h1err, u0tot, u0err, u0min, u0max, u1tot, u1err, u1min, u1max, uh_err
 
-  call build_reconstructions_1d( n0, h0, u0, CS%remapping_scheme, CS%degree, CS%boundary_extrapolation, &
-                                   ppoly_r_coefficients, ppoly_r_E, ppoly_r_S, iMethod )
+  call build_reconstructions_1d( CS, n0, h0, u0, ppoly_r_coefficients, ppoly_r_E, ppoly_r_S, iMethod )
 
   if (CS%check_reconstruction) call check_reconstructions_1d(n0, h0, u0, CS%degree, &
                                    CS%boundary_extrapolation, ppoly_r_coefficients, ppoly_r_E, ppoly_r_S)
@@ -243,8 +262,7 @@ subroutine remapping_core_w( CS, n0, h0, u0, n1, dx, u1 )
   real :: u0tot, u0err, u0min, u0max, u1tot, u1err, u1min, u1max, uh_err
   real, dimension(n1) :: h1 !< Cell widths on target grid
 
-  call build_reconstructions_1d( n0, h0, u0, CS%remapping_scheme, CS%degree, CS%boundary_extrapolation, &
-                                   ppoly_r_coefficients, ppoly_r_E, ppoly_r_S, iMethod )
+  call build_reconstructions_1d( CS, n0, h0, u0, ppoly_r_coefficients, ppoly_r_E, ppoly_r_S, iMethod )
 
   if (CS%check_reconstruction) call check_reconstructions_1d(n0, h0, u0, CS%degree, &
                                    CS%boundary_extrapolation, ppoly_r_coefficients, ppoly_r_E, ppoly_r_S)
@@ -301,28 +319,27 @@ subroutine remapping_core_w( CS, n0, h0, u0, n1, dx, u1 )
 end subroutine remapping_core_w
 
 !> Creates polynomial reconstructions of u0 on the source grid h0.
-subroutine build_reconstructions_1d( n0, h0, u0, remapping_scheme, deg, boundary_extrapolation, &
-                                     ppoly_r_coefficients, ppoly_r_E, ppoly_r_S, iMethod )
-  integer,                  intent(in)  :: n0 !< Number of cells on source grid
-  real, dimension(n0),      intent(in)  :: h0 !< Cell widths on source grid
-  real, dimension(n0),      intent(in)  :: u0 !< Cell averages on source grid
-  integer,                  intent(in)  :: remapping_scheme !< Remapping scheme
-  integer,                  intent(in)  :: deg !< Degree of polynomial reconstruction
-  logical,                  intent(in)  :: boundary_extrapolation !< Extrapolate at boundaries if true
-  real, dimension(n0,deg+1),intent(out) :: ppoly_r_coefficients !< Coefficients of polynomial
-  real, dimension(n0,2),    intent(out) :: ppoly_r_E !< Edge value of polynomial
-  real, dimension(n0,2),    intent(out) :: ppoly_r_S !< Edge slope of polynomial
-  integer,                  intent(out) :: iMethod !< Integration method
+subroutine build_reconstructions_1d( CS, n0, h0, u0, ppoly_r_coefficients, ppoly_r_E, ppoly_r_S, iMethod )
+  type(remapping_CS),              intent(in)  :: CS
+  integer,                         intent(in)  :: n0 !< Number of cells on source grid
+  real, dimension(n0),             intent(in)  :: h0 !< Cell widths on source grid
+  real, dimension(n0),             intent(in)  :: u0 !< Cell averages on source grid
+  real, dimension(n0,CS%degree+1), intent(out) :: ppoly_r_coefficients !< Coefficients of polynomial
+  real, dimension(n0,2),           intent(out) :: ppoly_r_E !< Edge value of polynomial
+  real, dimension(n0,2),           intent(out) :: ppoly_r_S !< Edge slope of polynomial
+  integer,                         intent(out) :: iMethod !< Integration method
   ! Local variables
   integer :: local_remapping_scheme
+  integer :: remapping_scheme !< Remapping scheme
+  logical :: boundary_extrapolation !< Extrapolate at boundaries if true
 
   ! Reset polynomial
   ppoly_r_E(:,:) = 0.0
   ppoly_r_S(:,:) = 0.0
   ppoly_r_coefficients(:,:) = 0.0
   iMethod = -999
 
-  local_remapping_scheme = remapping_scheme
+  local_remapping_scheme = CS%remapping_scheme
   if (n0<=1) then
     local_remapping_scheme = REMAPPING_PCM
   elseif (n0<=3) then
@@ -336,37 +353,37 @@ subroutine build_reconstructions_1d( n0, h0, u0, remapping_scheme, deg, boundary
       iMethod = INTEGRATION_PCM
     case ( REMAPPING_PLM )
       call PLM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_coefficients )
-      if ( boundary_extrapolation) then
+      if ( CS%boundary_extrapolation ) then
         call PLM_boundary_extrapolation( n0, h0, u0, ppoly_r_E, ppoly_r_coefficients)
       end if
       iMethod = INTEGRATION_PLM
     case ( REMAPPING_PPM_H4 )
       call edge_values_explicit_h4( n0, h0, u0, ppoly_r_E )
       call PPM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_coefficients )
-      if ( boundary_extrapolation) then
+      if ( CS%boundary_extrapolation ) then
         call PPM_boundary_extrapolation( n0, h0, u0, ppoly_r_E, ppoly_r_coefficients )
       end if
       iMethod = INTEGRATION_PPM
     case ( REMAPPING_PPM_IH4 )
       call edge_values_implicit_h4( n0, h0, u0, ppoly_r_E )
       call PPM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_coefficients )
-      if ( boundary_extrapolation) then
+      if ( CS%boundary_extrapolation ) then
         call PPM_boundary_extrapolation( n0, h0, u0, ppoly_r_E, ppoly_r_coefficients )
       end if
       iMethod = INTEGRATION_PPM
     case ( REMAPPING_PQM_IH4IH3 )
       call edge_values_implicit_h4( n0, h0, u0, ppoly_r_E )
       call edge_slopes_implicit_h3( n0, h0, u0, ppoly_r_S )
       call PQM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_S, ppoly_r_coefficients )
-      if ( boundary_extrapolation) then
+      if ( CS%boundary_extrapolation ) then
         call PQM_boundary_extrapolation_v1( n0, h0, u0, ppoly_r_E, ppoly_r_S, ppoly_r_coefficients )
       end if
       iMethod = INTEGRATION_PQM
     case ( REMAPPING_PQM_IH6IH5 )
       call edge_values_implicit_h6( n0, h0, u0, ppoly_r_E )
       call edge_slopes_implicit_h5( n0, h0, u0, ppoly_r_S )
       call PQM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_S, ppoly_r_coefficients )
-      if ( boundary_extrapolation) then
+      if ( CS%boundary_extrapolation ) then
         call PQM_boundary_extrapolation_v1( n0, h0, u0, ppoly_r_E, ppoly_r_S, ppoly_r_coefficients )
       end if
       iMethod = INTEGRATION_PQM

src/ALE/polynomial_functions.F90

@@ -12,7 +12,7 @@ module polynomial_functions
 
 implicit none ; private
 
-public :: evaluation_polynomial, integration_polynomial
+public :: evaluation_polynomial, integration_polynomial, first_derivative_polynomial
 
 ! -----------------------------------------------------------------------------
 ! This module contains the following routines
@@ -51,6 +51,37 @@ real function evaluation_polynomial( coefficients, nb_coefficients, x )
 
 end function evaluation_polynomial
 
+!> Calculates the first derivative of a polynomial with coefficients as above
+!! evaluated at a point x
+real function first_derivative_polynomial( coefficients, nb_coefficients, x )
+! -----------------------------------------------------------------------------
+! The polynomial is defined by the coefficients contained in the
+! array of the same name, as follows: C(1) + C(2)x + C(3)x^2 + C(4)x^3 + ...
+! where C refers to the array 'coefficients'.
+! The number of coefficients is given by nb_coefficients and x
+! is the coordinate where the polynomial's derivative is to be evaluated.
+!
+! The function returns the value of the polynomial at x.
+! -----------------------------------------------------------------------------
+
+  ! Arguments
+  real, dimension(:), intent(in)        :: coefficients
+  integer, intent(in)                   :: nb_coefficients
+  real, intent(in)                      :: x
+
+  ! Local variables
+  integer                               :: k
+  real                                  :: f    ! value of polynomial at x
+
+  f = 0.0
+  do k = 2,nb_coefficients
+    f = f + REAL(k-1)*coefficients(k) * ( x**(k-2) )
+  end do
+
+  first_derivative_polynomial = f
+
+end function first_derivative_polynomial
+
 ! -----------------------------------------------------------------------------
 ! Exact integration of polynomial of degree n
 ! -----------------------------------------------------------------------------