diff --git a/example/2D-advection-diffusion.F90 b/example/2D-advection-diffusion.F90 index 1916a51..d197640 100644 --- a/example/2D-advection-diffusion.F90 +++ b/example/2D-advection-diffusion.F90 @@ -53,15 +53,17 @@ program advection_diffusion_2D advance_time: & block double precision :: dt = 1D-6 - s = s + dt * d_dt(s, v) + associate(s_half => s + (dt/2) * d_dt(s, v)) + s = s + dt * d_dt(s_half, v) + end associate end block advance_time associate( & scalar_file => s%to_file("scalar") & ,velocity_file => v%to_file("vector") & ) - call scalar_file%write_lines("example/scripts/scalar.csv") - call velocity_file%write_lines("example/scripts/velocity.csv") + call scalar_file%write_lines("example/scripts/scalar-adv-dif.csv") + call velocity_file%write_lines("example/scripts/velocity-adv-dif.csv") end associate end associate @@ -73,7 +75,7 @@ pure function d_dt(s, v) result(ds_dt) type(vector_2D_t), intent(in) :: v type(scalar_2D_t) ds_dt double precision, parameter :: D = 1D0 - ds_dt = .div. (D * .grad. s) - (v .dot. .grad. s) + ds_dt = .div. (D * .grad. s) - .div. (v * s) end function end program diff --git a/src/formal/divergence_2D_s.F90 b/src/formal/divergence_2D_s.F90 index d78f5b0..a8e756a 100644 --- a/src/formal/divergence_2D_s.F90 +++ b/src/formal/divergence_2D_s.F90 @@ -68,6 +68,20 @@ end procedure + module procedure divergence_2D_minus_divergence + call_julienne_assert(lhs%conformable(rhs)) + + difference%tensor_2D_t = tensor_2D_t( & + points = reshape([points_2D_t(lhs%points_(1,1,1,1)%values_ - rhs%points_(1,1,1,1)%values_)], shape = [1,1,1,1]) & + ,cells = lhs%cells_ & + ,x_min = lhs%x_min_ & + ,x_max = lhs%x_max_ & + ,order = lhs%order_ & + ) + call_julienne_assert(difference%consistent()) + + end procedure + module procedure construct_2D_divergence_from_vector_mold call_julienne_assert(mold%consistent()) diff --git a/src/formal/scalar_2D_s.F90 b/src/formal/scalar_2D_s.F90 index 9e2c9ff..96115bd 100644 --- a/src/formal/scalar_2D_s.F90 +++ b/src/formal/scalar_2D_s.F90 @@ -168,6 +168,34 @@ end procedure + module procedure scalar_2D_assign_divergence + + call_julienne_assert(lhs%conformable(rhs)) + + allocate(lhs%points_(1,1,1,1)) + allocate(lhs%points_(1,1,1,1)%values_(rhs%cells_(x_dir)+2, rhs%cells_(y_dir)+2)) + + associate( & + x_last => size(rhs%points_(1,1,1,1)%values_,x_dir) - 1 & + ,y_last => size(rhs%points_(1,1,1,1)%values_,y_dir) - 1 & + ) + lhs%points_(1,1,1,1)%values_(2:x_last-1, 2:y_last-1) = rhs%points_(1,1,1,1)%values_(2:x_last-1, 2:y_last-1) ! internal points + lhs%points_(1,1,1,1)%values_(1 , : ) = 0D0 ! x_min boundary + lhs%points_(1,1,1,1)%values_( x_last , : ) = 0D0 ! x_max boundary + lhs%points_(1,1,1,1)%values_( : , 1 ) = 0D0 ! y_min boundary + lhs%points_(1,1,1,1)%values_( : , y_last ) = 0D0 ! y_max boundary + end associate + + lhs%cells_ = rhs%cells_ + lhs%x_min_ = rhs%x_min_ + lhs%x_max_ = rhs%x_max_ + lhs%order_ = rhs%order_ + + call_julienne_assert(lhs%consistent()) + + end procedure + + module procedure scalar_2D_to_file type(string_t), allocatable :: lines(:) integer i, j, l, m, n, p, q diff --git a/src/formal/tensors_2D_m.F90 b/src/formal/tensors_2D_m.F90 index ba48bf3..efba343 100644 --- a/src/formal/tensors_2D_m.F90 +++ b/src/formal/tensors_2D_m.F90 @@ -9,7 +9,7 @@ module tensors_2D_m use julienne_m, only : file_t implicit none - private + !private public :: scalar_2D_t public :: vector_2D_t public :: gradient_2D_t @@ -85,6 +85,7 @@ pure module function construct_2D_tensor_from_components(points, cells, x_min, x private type(gradient_operator_1D_t) gradient_operator_1D_(space_dimension) contains + generic :: assignment(=) => scalar_2D_assign_divergence generic :: operator(.grad.) => scalar_2D_gradient generic :: operator(*) => scalar_2D_postmultiply_double, scalar_2D_premultiply_double generic :: operator(+) => scalar_2D_plus_scalar @@ -93,6 +94,7 @@ pure module function construct_2D_tensor_from_components(points, cells, x_min, x generic :: consistent => scalar_2D_consistent generic :: conformable => scalar_2D_conformable_scalar 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_gradient procedure, non_overridable, private :: scalar_2D_values @@ -143,9 +145,10 @@ pure module function construct_2D_scalar_from_components(tensor_2D, gradient_ope generic :: grid => vector_2D_grid generic :: consistent => vector_2D_consistent generic :: conformable => vector_2D_conformable_vector, vector_2D_conformable_scalar - generic :: co_located_components => vector_2D_co_located_components + generic :: to_centers_extended => vector_2D_to_centers_extended generic :: operator(.div.) => vector_2D_divergence generic :: operator(.dot.) => vector_2D_dot_vector + generic :: operator(*) => vector_2D_postmultiply_scalar, vector_2D_premultiply_scalar generic :: to_file => vector_2D_to_file procedure, non_overridable, private :: vector_2D_to_file procedure, non_overridable, private :: vector_2D_grid @@ -153,8 +156,10 @@ pure module function construct_2D_scalar_from_components(tensor_2D, gradient_ope 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_co_located_components + 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 + procedure, non_overridable, private, pass(vector_2D) :: vector_2D_premultiply_scalar end type interface vector_2D_t @@ -227,7 +232,7 @@ pure module function construct_2D_gradient_from_components(tensor_2D, divergence 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 + generic :: operator(-) => divergence_2D_minus_scalar, divergence_2D_minus_divergence generic :: to_file => divergence_2D_to_file procedure, non_overridable, private :: divergence_2D_to_file procedure, private, non_overridable :: divergence_2D_values @@ -236,6 +241,7 @@ pure module function construct_2D_gradient_from_components(tensor_2D, divergence 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 end type @@ -278,6 +284,13 @@ pure module function tensor_2D_conformable(self, tensor_2D) result(conformable) logical conformable end function + pure module subroutine scalar_2D_assign_divergence(lhs, rhs) + !! Assign 2D divergence to 2D scalar at internal points + implicit none + class(scalar_2D_t), intent(inout) :: lhs + type(divergence_2D_t), intent(in) :: rhs + end subroutine + pure module function scalar_2D_consistent(self) result(self_consistent) !! Assert components allocated and self-consistent, including sufficient accuracy for gradient operator implicit none @@ -384,7 +397,7 @@ pure module function divergence_2D_grid(self, direction) result(divergence_grid_ double precision, allocatable :: divergence_grid_1D(:) !! grid points along the requested coordinate direction end function - pure module function vector_2D_co_located_components(self) result(vectors) + pure module function vector_2D_to_centers_extended(self) result(vectors) !! Vector values getter implicit none class(vector_2D_t), intent(in) :: self @@ -412,6 +425,22 @@ pure module function vector_2D_divergence(self) result(divergence_2D) type(divergence_2D_t) divergence_2D end function + pure module function vector_2D_postmultiply_scalar(vector_2D, scalar_2D) result(vector_x_scalar) + !! Result is product of the 2D vector and scalar arguments + implicit none + class(vector_2D_t), intent(in) :: vector_2D + type(scalar_2D_t), intent(in) :: scalar_2D + type(vector_2D_t) vector_x_scalar + end function + + pure module function vector_2D_premultiply_scalar(scalar_2D, vector_2D) result(scalar_x_vector) + !! Result is product of the 2D vector and scalar arguments + implicit none + class(vector_2D_t), intent(in) :: vector_2D + type(scalar_2D_t), intent(in) :: scalar_2D + type(vector_2D_t) scalar_x_vector + end function + pure module function vector_2D_dot_vector(lhs, rhs) result(scalar_2D) !! Result is scalar product of the 2D-vector arguments implicit none @@ -451,6 +480,13 @@ pure module function divergence_2D_premultiply_constant(lhs, rhs) result(lhs_x_r type(divergence_2D_t) lhs_x_rhs end function + pure module function divergence_2D_minus_divergence(lhs, rhs) result(difference) + !! Result is the pointwise difference between the lhs and rhs + implicit none + class(divergence_2D_t), intent(in) :: lhs, rhs + type(divergence_2D_t) difference + end function + pure module function divergence_2D_minus_scalar(lhs, rhs) result(difference) !! Result is the pointwise difference between the lhs and rhs implicit none diff --git a/src/formal/vector_2D_s.F90 b/src/formal/vector_2D_s.F90 index a225e0a..240c360 100644 --- a/src/formal/vector_2D_s.F90 +++ b/src/formal/vector_2D_s.F90 @@ -156,7 +156,7 @@ pure function description(coordinate, component) result(point_cloud) end procedure - module procedure vector_2D_co_located_components + module procedure vector_2D_to_centers_extended call_julienne_assert(self%consistent()) @@ -179,13 +179,46 @@ pure function description(coordinate, component) result(point_cloud) end procedure + module procedure vector_2D_postmultiply_scalar + + 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_))) + + vector_x_scalar = construct_2D_vector_from_components( & + 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_)] & + ,shape = [space_dimension,1,1,1] & + ) & + ,cells = vector_2D%cells_ & + ,x_min = vector_2D%x_min_ & + ,x_max = vector_2D%x_max_ & + ,order = vector_2D%order_ & + ) & + ,divergence_operator_1D_t( & + k = vector_2D%order_ & + ,dx = (vector_2D%x_max_ - vector_2D%x_min_)/vector_2D%cells_ & + ,cells = vector_2D%cells_ & + ) ) + end associate + + end procedure + + module procedure vector_2D_premultiply_scalar + scalar_x_vector = vector_2D * scalar_2D + end procedure + module procedure vector_2D_dot_vector call_julienne_assert(lhs%conformable(rhs)) associate( & - lhs_ => lhs%co_located_components() & - ,rhs_ => rhs%co_located_components() & + lhs_ => lhs%to_centers_extended() & + ,rhs_ => rhs%to_centers_extended() & ) call_julienne_assert(.all. (shape(lhs_) .equalsExpected. shape(rhs_))) @@ -217,7 +250,7 @@ pure function description(coordinate, component) result(point_cloud) header => [string_t("x, y, " // name)] & ,x => cell_centers_extended_1D(self%x_min_(x_dir), self%x_max_(x_dir), self%cells_(x_dir)) & ,y => cell_centers_extended_1D(self%x_min_(y_dir), self%x_max_(y_dir), self%cells_(y_dir)) & - ,vectors => self%co_located_components() & + ,vectors => self%to_centers_extended() & ) allocate(lines(size(header) + size(x)*size(y))) diff --git a/test/scalar_2D_test_m.F90 b/test/scalar_2D_test_m.F90 index 8da5816..ef9fec8 100644 --- a/test/scalar_2D_test_m.F90 +++ b/test/scalar_2D_test_m.F90 @@ -77,7 +77,7 @@ function check_gradient() result(test_diagnosis) ,expected_gradient => vector_2D_t(expected_gradient_initializer, mold=scalar_2D) & ) test_diagnosis = test_diagnosis .also. & - .all. (grad_scalar%co_located_components() .approximates. expected_gradient%co_located_components() .within. tolerance) & + .all. (grad_scalar%to_centers_extended() .approximates. expected_gradient%to_centers_extended() .within. tolerance) & // string_t(" for order ") // string_t(order) end associate end associate