More MLP updates, not done yet

.vscode/launch.json

@@ -43,5 +43,19 @@
             "externalConsole": false,
             "MIMode": "gdb",
         },
+        {
+            "name": "(gdb) Sedov 2D",
+            "type": "cppdbg",
+            "request": "launch",
+            "program": "${workspaceFolder}/build/bin/cato.x",
+            "args": [
+                "input.ini"
+            ],
+            "stopAtEntry": false,
+            "cwd": "${workspaceFolder}/tests/integrated/sedov",
+            "environment": [],
+            "externalConsole": false,
+            "MIMode": "gdb",
+        },
     ]
 }

src/lib/fluid/mod_fluid.f90

@@ -1,5 +1,5 @@
 ! MIT License
-! Copyright (c) 2019 Sam Miller
+! Copyright (c) 2020 Sam Miller
 ! Permission is hereby granted, free of charge, to any person obtaining a copy
 ! of this software and associated documentation files (the "Software"), to deal
 ! in the Software without restriction, including without limitation the rights
@@ -541,7 +541,7 @@ subroutine get_continuity_sensor(self)
     !< or discontinuous (linear or non-linear)
     class(fluid_t), intent(inout) :: self
 
-    call distinguish(rho=self%rho, u=self%u, v=self%v, p=self%p, continuity_sensor=self%continuous_sensor)
+    call distinguish(lbounds=lbound(self%rho), rho=self%rho, u=self%u, v=self%v, p=self%p, continuity_sensor=self%continuous_sensor)
   end subroutine get_continuity_sensor
 
   subroutine residual_smoother(self)

src/lib/io/contour_writer.f90

@@ -284,8 +284,8 @@ subroutine write_hdf5(self, master, time, iteration)
                                     description='Cell j Index', units='dimensionless')
     deallocate(int_data_buffer)
 
-    call self%write_2d_integer_data(data=master%fluid%continuous_sensor, name='/continuity_sensor', &
-                                description='Continuity Sensor [0=continuous, 1=linear discontinuity, 2=non-linear discontinuity', &
+    call self%write_2d_integer_data(data=master%fluid%continuous_sensor(ilo:ihi, jlo:jhi), name='/continuity_sensor', &
+                               description='Continuity Sensor [0=continuous, 1=linear discontinuity, 2=non-linear discontinuity]', &
                                     units='dimensionless')
 
     ! Primitive Variables

src/lib/io/mod_input.f90

@@ -287,23 +287,23 @@ subroutine read_from_ini(self, filename)
       call cfg%get("initial_conditions", "init_pressure", self%init_pressure, 0.0_rk)
     end if
 
-    ! Grid
-    select case(trim(self%edge_interpolation_scheme))
-    case('TVD2', 'TVD3', 'TVD5', 'MLP3', 'MLP5')
-      required_n_ghost_layers = 2
-      call cfg%get("grid", "n_ghost_layers", self%n_ghost_layers, required_n_ghost_layers)
-    case default
-      call error_msg(module='mod_input', class="input_t", procedure='read_from_ini', &
-                     message="Unknown edge interpolation scheme, must be one of the following: "// &
-                     "'TVD2', 'TVD3', 'TVD5', 'MLP3', or 'MLP5'", &
-                     file_name=__FILE__, line_number=__LINE__)
-    end select
-
-    if(self%n_ghost_layers /= required_n_ghost_layers) then
-      call error_msg(module='mod_input', class="input_t", procedure='read_from_ini', &
-                     message="The number of required ghost cell layers doesn't match the edge interpolation order", &
-                     file_name=__FILE__, line_number=__LINE__)
-    end if
+    ! ! Grid
+    ! select case(trim(self%edge_interpolation_scheme))
+    ! case('TVD2', 'TVD3', 'TVD5', 'MLP3', 'MLP5')
+    !   required_n_ghost_layers = 2
+    !   call cfg%get("grid", "n_ghost_layers", self%n_ghost_layers, required_n_ghost_layers)
+    ! case default
+    !   call error_msg(module='mod_input', class="input_t", procedure='read_from_ini', &
+    !                  message="Unknown edge interpolation scheme, must be one of the following: "// &
+    !                  "'TVD2', 'TVD3', 'TVD5', 'MLP3', or 'MLP5'", &
+    !                  file_name=__FILE__, line_number=__LINE__)
+    ! end select
+
+    ! if(self%n_ghost_layers /= required_n_ghost_layers) then
+    !   call error_msg(module='mod_input', class="input_t", procedure='read_from_ini', &
+    !                  message="The number of required ghost cell layers doesn't match the edge interpolation order", &
+    !                  file_name=__FILE__, line_number=__LINE__)
+    ! end if
 
     call cfg%get("grid", "grid_type", char_buffer, '2d_regular')
     self%grid_type = trim(char_buffer)

src/lib/limiters/e_mlp_distinguisher.fypp

@@ -18,6 +18,12 @@
 ! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 ! SOFTWARE.
 
+#ifdef __SIMD_ALIGN_OMP__
+#define __CONT_ALIGN__ aligned(rho, u, v, p, d_bar_rho, d_bar_u, d_bar_v, d_bar_p:__ALIGNBYTES__)
+#else
+#define __CONT_ALIGN__
+#endif
+
 module mod_distinguisher
   !< Summary: Provide the procedures to conduct the "distinguisher" step in e-MLP. This process
   !<          scans the domain and finds the regions of linear/non-linear discontinuities as well
@@ -39,12 +45,14 @@ module mod_distinguisher
 
 contains
 
-  subroutine distinguish(rho, u, v, p, continuity_sensor)
+  subroutine distinguish(lbounds, rho, u, v, p, continuity_sensor)
     !< Scan the domain for continuous and discontinuous regions
-    real(rk), dimension(:, :), contiguous, intent(in) :: rho !< (i,j); cell-centered density
-    real(rk), dimension(:, :), contiguous, intent(in) :: u   !< (i,j); cell-centered x-velocity
-    real(rk), dimension(:, :), contiguous, intent(in) :: v   !< (i,j); cell-centered y-velocity
-    real(rk), dimension(:, :), contiguous, intent(in) :: p   !< (i,j); cell-centered pressure
+
+    integer(ik), dimension(2), intent(in) :: lbounds
+    real(rk), dimension(lbounds(1):, lbounds(2):), contiguous, intent(in) :: rho !< (i,j); cell-centered density
+    real(rk), dimension(lbounds(1):, lbounds(2):), contiguous, intent(in) :: u   !< (i,j); cell-centered x-velocity
+    real(rk), dimension(lbounds(1):, lbounds(2):), contiguous, intent(in) :: v   !< (i,j); cell-centered y-velocity
+    real(rk), dimension(lbounds(1):, lbounds(2):), contiguous, intent(in) :: p   !< (i,j); cell-centered pressure
 
     integer(ik), dimension(:, :), allocatable, intent(out) :: continuity_sensor
     !< (i,j); sensor value to tag wheter it is continuous, linear discontinuous, or non-linear discontinuous
@@ -70,12 +78,18 @@ contains
     jhi = ubound(rho, dim=2)
 
     allocate(continuity_sensor(ilo:ihi, jlo:jhi))
+
     continuity_sensor = CONTINUOUS_REGION
 
     allocate(d_bar_rho(ilo:ihi, jlo:jhi))
+    !dir$ assume_aligned d_bar_rho: __ALIGNBYTES__
     allocate(d_bar_u(ilo:ihi, jlo:jhi))
+    !dir$ assume_aligned d_bar_u: __ALIGNBYTES__
     allocate(d_bar_v(ilo:ihi, jlo:jhi))
+    !dir$ assume_aligned d_bar_v: __ALIGNBYTES__
     allocate(d_bar_p(ilo:ihi, jlo:jhi))
+    !dir$ assume_aligned d_bar_p: __ALIGNBYTES__
+
     d_bar_rho = 0.0_rk
     d_bar_u = 0.0_rk
     d_bar_v = 0.0_rk
@@ -89,35 +103,36 @@ contains
         if(abs(${F}$(i, j)) > 0.0_rk) then
           ! Eq 11a
           d_ij_i = abs(((-one_sixth * ${F}$(i - 2, j) + two_thirds * ${F}$(i - 1, j) + two_thirds * ${F}$(i + 1, j) - one_sixth * ${F}$(i + 2, j)) &
-                        / ${F}$(i,j)) &
-                        - 1.0_rk)
+                        / ${F}$(i,j)) - 1.0_rk)
 
           ! Eq 11b
           d_ij_j = abs(((-one_sixth * ${F}$(i, j - 2) + two_thirds * ${F}$(i, j - 1) + two_thirds * ${F}$(i, j + 1) - one_sixth * ${F}$(i, j + 2)) &
-                         / ${F}$(i,j)) &
-                         - 1.0_rk)
+                        / ${F}$(i,j)) - 1.0_rk)
         end if
-        d_bar_${F}$(i, j) = 0.5_rk * (d_ij_i + d_ij_j)
         ! Eq 11c
+        d_bar_${F}$(i, j) = 0.5_rk * (d_ij_i + d_ij_j)
       end do
     end do
+
     #:endfor
 
     ! Assign contiuous sensor based on the approximate values, e.g. d_bar_rho
     do j = jlo, jhi
       do i = ilo, ihi
         if(abs(d_bar_rho(i, j)) > EPS) continuity_sensor(i, j) = LINEAR_DISCONT_REGION
 
-        if(abs(u(i, j)) < abs(v(i, j))) then
-          if(abs(d_bar_v(i, j)) > EPS .and. abs(v(i, j)) > 1e-6_rk) continuity_sensor(i, j) = LINEAR_DISCONT_REGION
-        else ! abs(u(i,j)) >= abs(v(i,j))
+        if(abs(u(i, j)) > abs(v(i, j)) .or.  abs(u(i, j) - v(i, j)) < epsilon(1.0_rk)) then
           if(abs(d_bar_u(i, j)) > EPS .and. abs(u(i, j)) > 1e-6_rk) continuity_sensor(i, j) = LINEAR_DISCONT_REGION
+
+        else if(abs(u(i, j)) < abs(v(i, j))) then
+          if(abs(d_bar_v(i, j)) > EPS .and. abs(v(i, j)) > 1e-6_rk) continuity_sensor(i, j) = LINEAR_DISCONT_REGION
         end if
 
         if(abs(d_bar_p(i, j)) > EPS) continuity_sensor(i, j) = NONLINEAR_DISCONT_REGION
       end do
     end do
 
+
     deallocate(d_bar_rho)
     deallocate(d_bar_u)
     deallocate(d_bar_v)