diff --git a/src/formal/tensors_3D_m.F90 b/src/formal/tensors_3D_m.F90 index 695ceed..41af54f 100644 --- a/src/formal/tensors_3D_m.F90 +++ b/src/formal/tensors_3D_m.F90 @@ -151,7 +151,7 @@ pure module function construct_3D_scalar_from_components(tensor_3D, gradient_ope generic :: to_centers_extended => vector_3D_to_centers_extended generic :: operator(.div.) => vector_3D_divergence generic :: operator(.dot.) => vector_3D_dot_vector - generic :: operator(*) => vector_3D_postmultiply_scalar, vector_3D_premultiply_scalar + generic :: operator(*) => vector_3D_postmultiply_scalar_3D, vector_3D_premultiply_scalar_3D generic :: to_file => vector_3D_to_file procedure, non_overridable, private :: vector_3D_to_file procedure, non_overridable, private :: vector_3D_grid @@ -160,8 +160,8 @@ pure module function construct_3D_scalar_from_components(tensor_3D, gradient_ope procedure, non_overridable, private :: vector_3D_consistent procedure, non_overridable, private :: vector_3D_to_centers_extended procedure, non_overridable, private :: vector_3D_dot_vector - procedure, non_overridable, private :: vector_3D_postmultiply_scalar - procedure, non_overridable, private, pass(vector_3D) :: vector_3D_premultiply_scalar + procedure, non_overridable, private :: vector_3D_postmultiply_scalar_3D + procedure, non_overridable, private, pass(vector_3D) :: vector_3D_premultiply_scalar_3D end type interface vector_3D_t @@ -418,7 +418,7 @@ pure module function vector_3D_divergence(self) result(divergence_3D) type(divergence_3D_t) divergence_3D end function - pure module function vector_3D_postmultiply_scalar(vector_3D, scalar_3D) result(vector_x_scalar) + pure module function vector_3D_postmultiply_scalar_3D(vector_3D, scalar_3D) result(vector_x_scalar) !! Result is product of the 3D vector and scalar arguments implicit none class(vector_3D_t), intent(in) :: vector_3D @@ -426,7 +426,7 @@ pure module function vector_3D_postmultiply_scalar(vector_3D, scalar_3D) result( type(vector_3D_t) vector_x_scalar end function - pure module function vector_3D_premultiply_scalar(scalar_3D, vector_3D) result(scalar_x_vector) + pure module function vector_3D_premultiply_scalar_3D(scalar_3D, vector_3D) result(scalar_x_vector) !! Result is product of the 3D vector and scalar arguments implicit none class(vector_3D_t), intent(in) :: vector_3D diff --git a/src/formal/vector_3D_s.F90 b/src/formal/vector_3D_s.F90 index c0c4dc0..056f97d 100644 --- a/src/formal/vector_3D_s.F90 +++ b/src/formal/vector_3D_s.F90 @@ -186,7 +186,7 @@ pure function description(coordinate, component) result(point_cloud) end procedure - module procedure vector_3D_postmultiply_scalar + module procedure vector_3D_postmultiply_scalar_3D call_julienne_assert(vector_3D%conformable(scalar_3D)) @@ -221,7 +221,7 @@ pure function description(coordinate, component) result(point_cloud) end procedure - module procedure vector_3D_premultiply_scalar + module procedure vector_3D_premultiply_scalar_3D scalar_x_vector = vector_3D * scalar_3D end procedure diff --git a/test/driver.f90 b/test/driver.f90 index 8cc6f35..35d7202 100644 --- a/test/driver.f90 +++ b/test/driver.f90 @@ -13,6 +13,7 @@ program test_suite_driver use scalar_2D_test_m, only : scalar_2D_test_t use scalar_3D_test_m, only : scalar_3D_test_t use vector_2D_test_m, only : vector_2D_test_t + use vector_3D_test_m, only : vector_3D_test_t implicit none associate(test_harness => test_harness_t([ & @@ -25,6 +26,7 @@ program test_suite_driver ,test_fixture_t(scalar_2D_test_t()) & ,test_fixture_t(scalar_3D_test_t()) & ,test_fixture_t(vector_2D_test_t()) & + ,test_fixture_t(vector_3D_test_t()) & ])) call test_harness%report_results end associate diff --git a/test/vector_2D_test_m.F90 b/test/vector_2D_test_m.F90 index d4dd283..0ee6a60 100644 --- a/test/vector_2D_test_m.F90 +++ b/test/vector_2D_test_m.F90 @@ -52,7 +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)) & + ,test_description_t('computing the product of a vector field and a scalar field', usher(check_vector_scalar_product)) & ]) end function diff --git a/test/vector_3D_test_m.F90 b/test/vector_3D_test_m.F90 new file mode 100644 index 0000000..5e39137 --- /dev/null +++ b/test/vector_3D_test_m.F90 @@ -0,0 +1,173 @@ +! Copyright (c) 2026, The Regents of the University of California +! Terms of use are as specified in LICENSE.txt + +#include "julienne-assert-macros.h" + +module vector_3D_test_m + use julienne_m, only : & + call_julienne_assert_ & + ,operator(//) & + ,operator(.all.) & + ,operator(.also.) & + ,operator(.approximates.) & + ,operator(.equalsExpected.) & + ,operator(.within.) & + ,passing_test & + ,string_t & + ,test_description_t & + ,test_diagnosis_t & + ,test_result_t & + ,test_t & + ,usher + use formal_m, only : & + scalar_3D_t & + ,scalar_3D_initializer_i & + ,vector_3D_t & + ,vector_3D_initializer_i & + ,divergence_3D_t & + ,divergence_3D_initializer_i & + ,x_dir & + ,y_dir & + ,z_dir + + implicit none + + type, extends(test_t) :: vector_3D_test_t + contains + procedure, nopass :: subject + procedure, nopass :: results + end type + + double precision, parameter :: tolerance = 1D-12 + +contains + + pure function subject() result(test_subject) + character(len=:), allocatable :: test_subject + test_subject = 'The vector_3D_t derived type' + end function + + function results() result(test_results) + type(vector_3D_test_t) vector_3D_test + type(test_result_t), allocatable :: test_results(:) + + test_results = vector_3D_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 field', usher(check_vector_scalar_product)) & + ]) + end function + + pure function rotated_stagnation_point_potential(x,y,z) result(s) + !! Define a stagnation-point scalar potential in a plane tilted 45 deg from the x-y plane + double precision, intent(in), dimension(:) :: x, y, z + double precision s(size(x),size(y),size(z)) + double precision, parameter :: pi = acos(-1D0) + do concurrent(integer :: j=1:size(y), k=1:size(z)) default(none) shared(x,y,z,s) + associate(eta => y(j)*cos(pi/4) + z(k)*sin(pi/4)) ! x-eta plane rotated around x axis pi/4 radians from x-y plane + s(:,j,k) = (x**2 - eta**2)/2 + end associate + end do + end function + + pure function rotated_stagnation_point_velocity(x,y,z) result(gradient) + !! Define a stagnation-point velocity field as the gradient of the stagantion-point scalar potential: + !! gradient(x,eta) = [ds/dx, ds/dy, ds/dz] + !! = [ x, (ds/deta)(deta/dy), (ds/deta)(deta/dz)] + !! = [ x, -eta * cos(theta), -eta * sin(theta)] + double precision, intent(in), dimension(:) :: x, y, z + double precision, parameter :: pi = acos(-1D0), theta = pi/4 + integer, parameter :: dimensionality = 3 + double precision gradient(size(x),size(y),size(z),dimensionality) + do concurrent(integer :: i=1:size(x), j=1:size(y), k=1:size(z)) default(none) shared(gradient,x,y,z) + associate(eta => y(j)*cos(theta) + z(k)*sin(theta)) ! x-eta plane rotated around x axis theta radians from x-y plane + gradient(i,j,k,:) = [x(i), -eta * cos(theta), -eta * sin(theta)] + end associate + end do + end function + + pure function velocity_squared(x,y,z) result(v_sq) + double precision, intent(in), dimension(:) :: x, y, z + double precision, parameter :: pi = acos(-1D0), theta = pi/4 + integer, parameter :: dimensionality = 3 + double precision v_sq(size(x),size(y),size(z)) + do concurrent(integer :: i=1:size(x), j=1:size(y), k=1:size(z)) default(none) shared(v_sq, x, y, z) + associate(eta => y(j)*cos(theta) + z(k)*sin(theta)) ! x-eta plane rotated around x axis theta radians from x-y plane + associate(v => [x(i), -eta * cos(theta), -eta * sin(theta)]) + v_sq(i,j,k) = dot_product(v,v) + end associate + end associate + end do + end function + + function check_divergence() result(test_diagnosis) + type(test_diagnosis_t) test_diagnosis + procedure(vector_3D_initializer_i), pointer :: vector_3D_initializer + procedure(divergence_3D_initializer_i), pointer :: expected_divergence_initializer + integer order + + test_diagnosis = passing_test() + vector_3D_initializer => rotated_stagnation_point_velocity + !expected_divergence_initializer => cubic_divergence + + do order = 2, 4, 2 + associate(vector_3D => vector_3D_t(vector_3D_initializer, order, cells=[20,20,20], x_min=[-100D0,-10D0,-10D0], x_max=[10D0,10D0,10D0])) + associate(div_vector => .div. vector_3D) + !associate(expected_divergence => divergence_3D_t(expected_divergence_initializer, mold=vector_3D)) + test_diagnosis = test_diagnosis .also. & + (.all. (div_vector%values() .approximates. 0D0 .within. tolerance)) & + // string_t(" for order ") // string_t(order) + !end associate + end associate + end associate + end do + end function + + function check_dot_product() result(test_diagnosis) + type(test_diagnosis_t) test_diagnosis + procedure(vector_3D_initializer_i), pointer :: v_init + procedure(scalar_3D_initializer_i), pointer :: v_sq_init + integer order + + test_diagnosis = passing_test() + + v_init => rotated_stagnation_point_velocity + v_sq_init => velocity_squared + + do order = 2, 4, 2 + associate( & + v => vector_3D_t(v_init, order, cells=[20,20,20], x_min=[-10D0, -10D0, -10D0], x_max=[10D0, 10D0, 10D0]) & + ,v_squared => scalar_3D_t(v_sq_init, order, cells=[20,20,20], x_min=[-10D0, -10D0, -10D0], x_max=[10D0, 10D0, 10D0]) & + ) + associate(v_dot_v => v .dot. v) + test_diagnosis = test_diagnosis .also. (.all. (v_dot_v%values() .approximates. v_squared%values() .within. tolerance)) + end associate + end associate + end do + end function + + function check_vector_scalar_product() result(test_diagnosis) + type(test_diagnosis_t) test_diagnosis + procedure(scalar_3D_initializer_i), pointer :: s_init + procedure(vector_3D_initializer_i), pointer :: v_init + integer order + + test_diagnosis = passing_test() + + v_init => rotated_stagnation_point_velocity + s_init => rotated_stagnation_point_potential + + do order = 2, 4, 2 + associate(s => scalar_3D_t(s_init, order, cells=[10,10,10], x_min=[-10D0, -10D0, -10D0], x_max=[10D0, 10D0, 10D0])) + associate(v => vector_3D_t(v_init, mold = s)) + associate(vs => v * s) + test_diagnosis = test_diagnosis .also. (.all. (vs%values(x_dir) .approximates. v%values(x_dir)*s%to_faces(x_dir) .within. tolerance)) + test_diagnosis = test_diagnosis .also. (.all. (vs%values(y_dir) .approximates. v%values(y_dir)*s%to_faces(y_dir) .within. tolerance)) + test_diagnosis = test_diagnosis .also. (.all. (vs%values(z_dir) .approximates. v%values(z_dir)*s%to_faces(z_dir) .within. tolerance)) + end associate + end associate + end associate + end do + end function + +end module vector_3D_test_m \ No newline at end of file