diff --git a/src/formal/divergence_2D_s.F90 b/src/formal/divergence_2D_s.F90 index a8e756a..a910866 100644 --- a/src/formal/divergence_2D_s.F90 +++ b/src/formal/divergence_2D_s.F90 @@ -16,18 +16,6 @@ contains - module procedure divergence_2D_conformable_scalar - call_julienne_assert(scalar_2D_consistent(scalar_2D)) - call_julienne_assert(self%tensor_2D_conformable(scalar_2D)) - conformable = .true. - end procedure - - module procedure divergence_2D_conformable_vector - call_julienne_assert(vector_2D_consistent(vector_2D)) - call_julienne_assert(self%tensor_2D_conformable(vector_2D)) - conformable = .true. - end procedure - module procedure divergence_2D_values call_julienne_assert(self%consistent()) divergences = self%points_(1,1,1,1)%values_(:,:) diff --git a/src/formal/scalar_2D_s.F90 b/src/formal/scalar_2D_s.F90 index 96115bd..2eb5076 100644 --- a/src/formal/scalar_2D_s.F90 +++ b/src/formal/scalar_2D_s.F90 @@ -15,6 +15,7 @@ ,operator(.isAtLeast.) & ,string_t use tensors_1D_m, only : cell_centers_extended_1D, scalar_1D_t + use interpolator_1D_m, only : centers_to_faces_1D_t implicit none contains @@ -26,13 +27,6 @@ self_consistent = .true. end procedure - module procedure scalar_2D_conformable_scalar - call_julienne_assert(scalar_2D_consistent(self)) - call_julienne_assert(scalar_2D_consistent(scalar_2D)) - call_julienne_assert(self%tensor_2D_conformable(scalar_2D)) - conformable = .true. - end procedure - module procedure construct_2D_scalar_from_components call_julienne_assert(size(gradient_operator_1D) .equalsExpected. space_dimension) @@ -170,9 +164,11 @@ module procedure scalar_2D_assign_divergence - call_julienne_assert(lhs%conformable(rhs)) + call_julienne_assert(rhs%consistent()) + if (allocated(lhs%points_)) deallocate(lhs%points_) allocate(lhs%points_(1,1,1,1)) + if (allocated(lhs%points_(1,1,1,1)%values_)) deallocate(lhs%points_(1,1,1,1)%values_) allocate(lhs%points_(1,1,1,1)%values_(rhs%cells_(x_dir)+2, rhs%cells_(y_dir)+2)) associate( & @@ -192,6 +188,7 @@ lhs%order_ = rhs%order_ call_julienne_assert(lhs%consistent()) + call_julienne_assert(lhs%conformable(rhs)) end procedure @@ -235,4 +232,39 @@ file = file_t(lines) end procedure + module procedure scalar_2D_to_faces + + call_julienne_assert(self%consistent()) + + construct_interpolator_array: & + associate(interpolator => centers_to_faces_1D_t(order=self%order_, cells=self%cells_, dx=(self%x_max_ - self%x_min_)/self%cells_)) + + select case(direction) + + case(x_dir) + + allocate(scalars(self%cells_(x_dir)+1, self%cells_(y_dir)+2)) + + interpolate_centers_to_x_faces: & + do concurrent(integer :: j = 1:size(self%points_(1,1,1,1)%values_,y_dir)) + scalars(:,j) = interpolator(x_dir)%face_values(self%points_(1,1,1,1)%values_(:,j)) + end do interpolate_centers_to_x_faces + + case(y_dir) + + allocate(scalars(self%cells_(x_dir)+2, self%cells_(y_dir)+1)) + + interpolate_centers_to_y_faces: & + do concurrent(integer :: i = 1:size(self%points_(1,1,1,1)%values_,x_dir)) + scalars(i,:) = interpolator(y_dir)%face_values(self%points_(1,1,1,1)%values_(i,:)) + end do interpolate_centers_to_y_faces + + case default + error stop "scalar_2D_to_faces in scalar_2D_s: invalid direction" + end select + + end associate construct_interpolator_array + + end procedure scalar_2D_to_faces + end submodule scalar_2D_s \ No newline at end of file diff --git a/src/formal/tensors_2D_m.F90 b/src/formal/tensors_2D_m.F90 index efba343..9349d52 100644 --- a/src/formal/tensors_2D_m.F90 +++ b/src/formal/tensors_2D_m.F90 @@ -17,8 +17,10 @@ module tensors_2D_m public :: scalar_2D_initializer_i public :: vector_2D_initializer_i public :: divergence_2D_initializer_i + public :: x_dir + public :: y_dir - integer, parameter :: space_dimension = 2, max_tensor_rank = 4, x_dir = 1, y_dir = 2, z_dir = 3 + integer, parameter :: space_dimension = 2, max_tensor_rank = 4, x_dir = 1, y_dir = 2 abstract interface @@ -61,6 +63,7 @@ pure function vector_2D_initializer_i(x,y) result(v) integer cells_(space_dimension) !! number of grid cells spanning the domain integer order_ !! order of accuracy of mimetic discretization contains + generic :: conformable => tensor_2D_conformable procedure, non_overridable, private :: tensor_rank procedure, non_overridable, private :: tensor_2D_consistent procedure, non_overridable, private :: tensor_2D_conformable @@ -92,15 +95,15 @@ pure module function construct_2D_tensor_from_components(points, cells, x_min, x generic :: values => scalar_2D_values generic :: grid => scalar_2D_grid generic :: consistent => scalar_2D_consistent - generic :: conformable => scalar_2D_conformable_scalar + generic :: to_faces => scalar_2D_to_faces generic :: to_file => scalar_2D_to_file procedure, non_overridable, private :: scalar_2D_assign_divergence procedure, non_overridable, private :: scalar_2D_to_file + procedure, non_overridable, private :: scalar_2D_to_faces procedure, non_overridable, private :: scalar_2D_gradient procedure, non_overridable, private :: scalar_2D_values procedure, non_overridable, private :: scalar_2D_grid procedure, non_overridable, private :: scalar_2D_consistent - procedure, non_overridable, private :: scalar_2D_conformable_scalar procedure, non_overridable, private :: scalar_2D_postmultiply_double procedure, non_overridable, private :: scalar_2D_plus_scalar procedure, non_overridable, private, pass(rhs) :: scalar_2D_premultiply_double @@ -143,8 +146,8 @@ pure module function construct_2D_scalar_from_components(tensor_2D, gradient_ope type(divergence_operator_1D_t) divergence_operator_1D_(space_dimension) contains generic :: grid => vector_2D_grid + generic :: values => vector_2D_values generic :: consistent => vector_2D_consistent - generic :: conformable => vector_2D_conformable_vector, vector_2D_conformable_scalar generic :: to_centers_extended => vector_2D_to_centers_extended generic :: operator(.div.) => vector_2D_divergence generic :: operator(.dot.) => vector_2D_dot_vector @@ -152,10 +155,9 @@ pure module function construct_2D_scalar_from_components(tensor_2D, gradient_ope generic :: to_file => vector_2D_to_file procedure, non_overridable, private :: vector_2D_to_file procedure, non_overridable, private :: vector_2D_grid + procedure, non_overridable, private :: vector_2D_values procedure, non_overridable, private :: vector_2D_divergence procedure, non_overridable, private :: vector_2D_consistent - procedure, non_overridable, private :: vector_2D_conformable_vector - procedure, non_overridable, private :: vector_2D_conformable_scalar procedure, non_overridable, private :: vector_2D_to_centers_extended procedure, non_overridable, private :: vector_2D_dot_vector procedure, non_overridable, private :: vector_2D_postmultiply_scalar @@ -230,7 +232,6 @@ pure module function construct_2D_gradient_from_components(tensor_2D, divergence generic :: values => divergence_2D_values generic :: grid => divergence_2D_grid generic :: consistent => tensor_2D_consistent - generic :: conformable => divergence_2D_conformable_scalar, divergence_2D_conformable_vector generic :: operator(*) => divergence_2D_premultiply_constant, divergence_2D_postmultiply_constant generic :: operator(-) => divergence_2D_minus_scalar, divergence_2D_minus_divergence generic :: to_file => divergence_2D_to_file @@ -238,8 +239,6 @@ pure module function construct_2D_gradient_from_components(tensor_2D, divergence procedure, private, non_overridable :: divergence_2D_values procedure, private, non_overridable :: divergence_2D_grid procedure, private, non_overridable :: divergence_2D_minus_scalar - procedure, private, non_overridable :: divergence_2D_conformable_vector - procedure, private, non_overridable :: divergence_2D_conformable_scalar procedure, private, non_overridable :: divergence_2D_postmultiply_constant procedure, private, non_overridable :: divergence_2D_minus_divergence procedure, private, non_overridable, pass(rhs) :: divergence_2D_premultiply_constant @@ -298,49 +297,12 @@ pure module function scalar_2D_consistent(self) result(self_consistent) logical self_consistent end function - pure module function scalar_2D_conformable_scalar(self, scalar_2D) result(conformable) - !! Assert the arguments' components are conformable, self-consistent, and consistent with each other + pure module function scalar_2D_to_faces(self, direction) result(scalars) + !! Result is the scalar values interpolated form centers-extended to faces along the requested direction implicit none - class(scalar_2D_t), intent(in) :: self, scalar_2D - logical conformable - end function - - pure module function divergence_2D_conformable_scalar(self, scalar_2D) result(conformable) - !! Assert the arguments' components are conformable, self-consistent, and consistent with each other - implicit none - class(divergence_2D_t), intent(in) :: self - class(scalar_2D_t), intent(in) :: scalar_2D - logical conformable - end function - - pure module function divergence_2D_conformable_vector(self, vector_2D) result(conformable) - !! Assert the arguments' components are conformable, self-consistent, and consistent with each other - implicit none - class(divergence_2D_t), intent(in) :: self - class(vector_2D_t), intent(in) :: vector_2D - logical conformable - end function - - pure module function vector_2D_consistent(self) result(self_consistent) - !! Assert components allocated and self-consistent, including sufficient accuracy for divergence operator - implicit none - class(vector_2D_t), intent(in) :: self - logical self_consistent - end function - - pure module function vector_2D_conformable_vector(self, vector_2D) result(conformable) - !! Assert the arguments' components are conformable, self-consistent, and consistent with each other - implicit none - class(vector_2D_t), intent(in) :: self, vector_2D - logical conformable - end function - - pure module function vector_2D_conformable_scalar(self, scalar_2D) result(conformable) - !! Assert components allocated and self-consistent - implicit none - class(vector_2D_t), intent(in) :: self - class(scalar_2D_t), intent(in) :: scalar_2D - logical conformable + class(scalar_2D_t), intent(in) :: self + integer, intent(in) :: direction + double precision, allocatable :: scalars(:,:) end function pure module function scalar_2D_values(self) result(values) @@ -389,6 +351,14 @@ pure module function vector_2D_grid(self, component, coordinate) result(vector_g double precision, allocatable :: vector_grid_1D(:) !! grid points along the requested coordinate direction end function + pure module function vector_2D_values(self, direction) result(vector_values) + !! Result contains the vector values for the component designated by "direction" + implicit none + class(vector_2D_t), intent(in) :: self + integer, intent(in) :: direction + double precision, allocatable :: vector_values(:,:) + end function + pure module function divergence_2D_grid(self, direction) result(divergence_grid_1D) !! Result array contains divergence grid locations along the requested spatial direction implicit none @@ -404,6 +374,13 @@ pure module function vector_2D_to_centers_extended(self) result(vectors) double precision, allocatable :: vectors(:,:,:) end function + pure module function vector_2D_consistent(self) result(self_consistent) + !! Assert internal self consistency, including sufficient accuracy for divergence operator + implicit none + class(vector_2D_t), intent(in) :: self + logical self_consistent + end function + pure module function divergence_2D_values(self) result(divergences) !! Vector values getter implicit none diff --git a/src/formal/vector_2D_s.F90 b/src/formal/vector_2D_s.F90 index 240c360..b1e5101 100644 --- a/src/formal/vector_2D_s.F90 +++ b/src/formal/vector_2D_s.F90 @@ -41,18 +41,8 @@ self_consistent = .true. end procedure - module procedure vector_2D_conformable_vector - call_julienne_assert(vector_2D_consistent(self)) - call_julienne_assert(vector_2D_consistent(vector_2D)) - call_julienne_assert(self%tensor_2D_conformable(vector_2D)) - conformable = .true. - end procedure - - module procedure vector_2D_conformable_scalar - call_julienne_assert(vector_2D_consistent(self)) - call_julienne_assert(scalar_2D_consistent(scalar_2D)) - call_julienne_assert(self%tensor_2D_conformable(scalar_2D)) - conformable = .true. + module procedure vector_2D_values + vector_values = self%points_(direction,1,1,1)%values_ end procedure module procedure construct_2D_vector_from_function @@ -183,15 +173,18 @@ pure function description(coordinate, component) result(point_cloud) call_julienne_assert(vector_2D%conformable(scalar_2D)) - associate(vector => vector_2D%to_centers_extended()) - - call_julienne_assert(.all. ([size(vector,x_dir), size(vector,y_dir)] .equalsExpected. shape(scalar_2D%points(1,1,1,1)%values_))) + associate( & + scalar_x => scalar_2D%to_faces(x_dir) & + ,scalar_y => scalar_2D%to_faces(y_dir) & + ) + call_julienne_assert(.all. (shape(scalar_x) .equalsExpected. shape(vector_2D%points_(x_dir,1,1,1)%values_))) + call_julienne_assert(.all. (shape(scalar_y) .equalsExpected. shape(vector_2D%points_(y_dir,1,1,1)%values_))) - vector_x_scalar = construct_2D_vector_from_components( & + vector_x_scalar = vector_2D_t( & tensor_2D_t( & points = reshape( & - source = [ points_2D_t(vector(:,:,x_dir) * scalar_2D%points_(1,1,1,1)%values_) & - ,points_2D_t(vector(:,:,y_dir) * scalar_2D%points_(1,1,1,1)%values_)] & + source = [ points_2D_t(vector_2D%points_(x_dir,1,1,1)%values_ * scalar_x) & + ,points_2D_t(vector_2D%points_(y_dir,1,1,1)%values_ * scalar_y)] & ,shape = [space_dimension,1,1,1] & ) & ,cells = vector_2D%cells_ & diff --git a/src/formal_m.f90 b/src/formal_m.f90 index 355e05b..bc49cc6 100644 --- a/src/formal_m.f90 +++ b/src/formal_m.f90 @@ -25,7 +25,8 @@ module formal_m ,gradient_2D_t & ! result of `.grad. s` for a scalar_2D_t s ,scalar_2D_initializer_i & ! scalar_2D_t initializer abstract interface ,vector_2D_initializer_i & ! vector_2D_t initializar abstract interface - ,divergence_2D_initializer_i ! divergence_2D_t initializar abstract interface + ,divergence_2D_initializer_i & ! divergence_2D_t initializar abstract interface + ,x_dir, y_dir use differential_operators_1D_m, only : & gradient_operator_1D_t & ! matrix operator defining a 1D mimetic gradient diff --git a/test/vector_2D_test_m.F90 b/test/vector_2D_test_m.F90 index 0ccaa37..d4dd283 100644 --- a/test/vector_2D_test_m.F90 +++ b/test/vector_2D_test_m.F90 @@ -15,12 +15,16 @@ module vector_2D_test_m ,test_diagnosis_t & ,test_result_t & ,test_t & - ,usher + ,usher use formal_m, only : & - vector_2D_t & + scalar_2D_t & + ,scalar_2D_initializer_i & + ,vector_2D_t & ,vector_2D_initializer_i & ,divergence_2D_t & - ,divergence_2D_initializer_i + ,divergence_2D_initializer_i & + ,x_dir & + ,y_dir implicit none @@ -48,6 +52,7 @@ function results() result(test_results) test_results = vector_2D_test%run([ & test_description_t('computing the divergence of a vector field', usher(check_divergence)) & ,test_description_t('computing the dot product of two vector fields', usher(check_dot_product)) & + ,test_description_t('computing the product of a vector field and a scalar', usher(check_vector_scalar_product)) & ]) end function @@ -158,4 +163,54 @@ function check_dot_product() result(test_diagnosis) end do end function + pure function scalar_field(x,y) result(s) + double precision, intent(in) :: x(:), y(:) + double precision s(size(x),size(y)) + do concurrent(integer :: i=1:size(x), j=1:size(y)) + s(i,j) = x(i) + end do + end function + + pure function vector_field(x,y) result(v) + double precision, intent(in) :: x(:), y(:) + double precision v(size(x),size(y),space_dimension) + do concurrent(integer :: i=1:size(x), j=1:size(y)) + v(i,j,:) = y(j) + end do + end function + + pure function scalar_vector_product(x,y) result(u) + double precision, intent(in) :: x(:), y(:) + double precision u(size(x),size(y),space_dimension) + do concurrent(integer :: i=1:size(x), j=1:size(y)) + u(i,j,:) = x(i)*y(j) + end do + end function + + function check_vector_scalar_product() result(test_diagnosis) + type(test_diagnosis_t) test_diagnosis + procedure(scalar_2D_initializer_i), pointer :: s_init + procedure(vector_2D_initializer_i), pointer :: v_init, vs_init + integer order + + test_diagnosis = passing_test() + + s_init => scalar_field + v_init => vector_field + vs_init => scalar_vector_product + + do order = 2, 4, 2 + associate( & + s => scalar_2D_t(s_init, cells=[10,10], x_min=[0D0,0D0], x_max=[5D0,5D0], order=order) & + ,v => vector_2D_t(v_init, cells=[10,10], x_min=[0D0,0D0], x_max=[5D0,5D0], order=order) & + ,vs_expected => vector_2D_t(vs_init, cells=[10,10], x_min=[0D0,0D0], x_max=[5D0,5D0], order=order) & + ) + associate(vs => v * s) + test_diagnosis = test_diagnosis .also. (.all. (vs%values(x_dir) .approximates. vs_expected%values(x_dir) .within. 1D-6)) + test_diagnosis = test_diagnosis .also. (.all. (vs%values(y_dir) .approximates. vs_expected%values(y_dir) .within. 1D-6)) + end associate + end associate + end do + end function + end module vector_2D_test_m \ No newline at end of file