Skip to content

Commit f9d08bb

Browse files
rousonrouson
committed
feat(interpolate): .div.dyad maps centers to faces
This commit adds a centers-to-faces interpolator and uses that interpolator when computing the divergence of a dyad.
1 parent a5ea34c commit f9d08bb

5 files changed

Lines changed: 89 additions & 36 deletions

File tree

‎example/burgers-1D.F90‎

Lines changed: 12 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -35,14 +35,13 @@ program burgers_1D
3535
stop new_line('') // new_line('') &
3636
// 'Usage:' // new_line('') &
3737
// ' fpm run \' // new_line('') &
38-
// ' --example burgers-1D \' // new_line('') &
39-
// ' --compiler flang-new \' // new_line('') &
40-
// ' --flag "-O3" \' // new_line('') &
41-
// ' -- [--help|-h] | [--order <integer>]' // new_line('') // new_line('') &
38+
// ' --example burgers-1D \' // new_line('') &
39+
// ' --compiler flang \' // new_line('') &
40+
// ' --flag "-O3" \' // new_line('') &
41+
// ' [--help|-h] | [--order <integer>]' // new_line('') // new_line('') &
4242
// 'where square brackets indicate optional arguments and angular brackets indicate user input values.' // new_line('')
4343
end if
4444

45-
4645
print *, new_line('')
4746
print *," Initial condition"
4847
print *," ================="
@@ -60,8 +59,14 @@ program burgers_1D
6059
print *, "order = ", order
6160

6261
u = vector_1D_t(vector_1D_initializer, order, x_min=0D0, x_max=1D0, cells=50)
63-
associate(div_uu_2 => .div. (u*u/2)) ! result is at cell centers; To Do: interpolate to cell faces
64-
end associate
62+
63+
block
64+
double precision dt
65+
dt = 1D0
66+
!associate(dt_grad_u2_2 => dt * (.div. ((u*u))))
67+
associate(dt_div_uu => .div. (u*u))
68+
end associate
69+
end block
6570

6671
#ifdef __GFORTRAN__
6772
stop

‎src/formal/dyad_1D_s.F90‎

Lines changed: 20 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,7 @@
55

66
submodule(tensors_1D_m) dyad_1D_s
77
use julienne_m, only : call_julienne_assert_
8+
use interpolator_1D_m, only : centers_to_faces_1D_t
89
implicit none
910
contains
1011

@@ -28,18 +29,26 @@
2829
associate(D => (self%divergence_operator_1D_))
2930
#endif
3031
call_julienne_assert(size(self%values_))
31-
associate(Dv => D .x. self%values_)
32-
divergence_1D%tensor_1D_t = tensor_1D_t(Dv(2:size(Dv)-1), self%x_min_, self%x_max_, self%cells_, self%order_)
33-
error stop "To Do: interploate the dyad's divergence to the cell faces and construct a vector_1D_t"
32+
associate( &
33+
Dv => D .x. self%values_ &
34+
,dx => (self%x_max_ - self%x_min_)/self%cells_ &
35+
)
36+
associate(interpolator => centers_to_faces_1D_t(order=self%order_, cells=self%cells_, dx=dx))
37+
vector_1D = vector_1D_t( &
38+
tensor_1D_t(interpolator .center. Dv(2:size(Dv)-1), self%x_min_, self%x_max_, self%cells_, self%order_) &
39+
,divergence_operator_1D_t(k=self%order_, dx=dx, cells=self%cells_) &
40+
)
41+
end associate
3442
#if ASSERTIONS
35-
associate( &
36-
q => divergence_1D%weights() &
37-
,dx => (self%x_max_ - self%x_min_)/self%cells_ &
38-
,b => [-1D0, [(0D0, center = 1, self%cells_-1)], 1D0] &
39-
)
40-
call_julienne_assert(.all. ([size(Dv), size(q)] .equalsExpected. self%cells_+2))
41-
call_julienne_assert((.all. (matmul(transpose(D%assemble()), q) .approximates. b/dx .within. double_equivalence)))
42-
! Check D^T * a = b_{m+1}, Eq. (19), Corbino & Castillo (2020)
43+
associate(divergence_1D => divergence_1D_t(tensor_1D_t(Dv(2:size(Dv)-1), self%x_min_, self%x_max_, self%cells_, self%order_)))
44+
associate( &
45+
q => divergence_1D%weights() &
46+
,b => [-1D0, [(0D0, center = 1, self%cells_-1)], 1D0] &
47+
)
48+
call_julienne_assert(.all. ([size(Dv), size(q)] .equalsExpected. self%cells_+2))
49+
call_julienne_assert((.all. (matmul(transpose(D%assemble()), q) .approximates. b/dx .within. double_equivalence)))
50+
! Check D^T * a = b_{m+1}, Eq. (19), Corbino & Castillo (2020)
51+
end associate
4352
end associate
4453
#endif
4554
end associate

src/formal/interpolation_operators_1D_m.F90 renamed to src/formal/interpolator_1D_m.F90

Lines changed: 26 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -1,35 +1,37 @@
11
! Copyright (c) 2026, The Regents of the University of California
22
! Terms of use are as specified in LICENSE.txt
33

4-
module interpolation_operators_1D_m
4+
module interpolator_1D_m
55
!! Define a sparse matrix storage format tailored to the staggered-grid interpolation matrix
66
!! operators detailed by Dumett & Castillo (2022) https://doi.org/10.13140/RG.2.2.26630.14400
77
implicit none
88

99
private
10-
public :: cells_to_faces_t
11-
public :: faces_to_cells_t
10+
public :: centers_to_faces_1D_t
1211

13-
type interpolation_operator_1D_t
12+
type, abstract :: interpolator_1D_t
1413
!! Encapsulate a staggered-grid interpolation matrix with a corresponding matrix-vector product operator
1514
private
16-
integer order
15+
integer order_, cells_, dx_
1716
double precision first_
1817
double precision, allocatable :: upper_(:,:)
1918
double precision, allocatable :: inner_(:)
2019
double precision, allocatable :: lower_(:,:)
2120
double precision final_
2221
end type
2322

24-
type, extends(interpolation_operator_1D_t) :: centers_to_faces_1D_t
23+
type, extends(interpolator_1D_t) :: centers_to_faces_1D_t
24+
contains
25+
generic :: operator(.center.) => interpolate_to_faces
26+
procedure, non_overridable, private :: interpolate_to_faces
2527
end type
2628

27-
type, extends(interpolation_operator_1D_t) :: faces_to_centers_1D_t
29+
type, extends(interpolator_1D_t) :: faces_to_centers_1D_t
2830
end type
2931

3032
interface centers_to_faces_1D_t
3133

32-
pure module function c2f_component_constructor(order, cells, dx) result(centers_to_faces_1D)
34+
pure module function c2f_constructor(order, cells, dx) result(centers_to_faces_1D)
3335
!! Construct centers-to-faces interpolation operator
3436
implicit none
3537
integer, intent(in) :: order, cells
@@ -41,14 +43,26 @@ pure module function c2f_component_constructor(order, cells, dx) result(centers_
4143

4244
interface faces_to_centers_1D_t
4345

44-
pure module function f2c_component_constructor(order, cells, dx) result(faces_to_centers_1D)
45-
!! Construct faces-to-centers interpolation operator
46+
pure module function f2c_constructor(order, cells, dx) result(faces_to_centers_1D)
47+
!! Construct centers-to-faces interpolation operator
4648
implicit none
4749
integer, intent(in) :: order, cells
4850
double precision, intent(in) :: dx
49-
type(faces_to_centers_1D_t) faces_to_centers_1D
51+
type(centers_to_faces_1D_t) faces_to_centers_1D
52+
end function
53+
54+
end interface
55+
56+
interface
57+
58+
pure module function interpolate_to_faces(self, centers) result(faces)
59+
!! Interpolate cell-centered values to face-centered values
60+
implicit none
61+
class(centers_to_faces_1D_t), intent(in) :: self
62+
double precision, intent(in) :: centers(:)
63+
double precision, allocatable :: faces(:)
5064
end function
5165

5266
end interface
5367

54-
end module interpolation_operators_1D_m
68+
end module interpolator_1D_m

src/formal/interpolation_operators_1D_s.F90 renamed to src/formal/interpolator_1D_s.F90

Lines changed: 29 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,14 +1,19 @@
11
! Copyright (c) 2026, The Regents of the University of California
22
! Terms of use are as specified in LICENSE.txt
33

4-
submodule(interpolation_operators_1D_m) interpolation_operators_1D_s
4+
#include "julienne-assert-macros.h"
5+
6+
submodule(interpolator_1D_m) interpolator_1D_s
7+
use julienne_m, only : call_julienne_assert_, operator(.equalsExpected.)
58
implicit none
69

710
contains
811

9-
module procedure c2f_component_constructor
12+
module procedure c2f_constructor
1013

14+
centers_to_faces_1D%order_ = order
1115
centers_to_faces_1D%cells_ = cells
16+
centers_to_faces_1D%dx_ = dx
1217

1318
select case(order)
1419
case(2)
@@ -27,11 +32,14 @@
2732
error stop "c2f_component_constructor: unsupported order"
2833
end select
2934

35+
call_julienne_assert(shape(self%lower_) .equalsExpected. shape(self%upper_))
3036
end procedure
3137

32-
module procedure f2c_component_constructor
38+
module procedure f2c_constructor
3339

40+
faces_to_centers_1D%order_ = order
3441
faces_to_centers_1D%cells_ = cells
42+
faces_to_centers_1D%dx_ = dx
3543

3644
select case(order)
3745
case(2)
@@ -50,6 +58,23 @@
5058
error stop "f2c_component_constructor: unsupported order"
5159
end select
5260

61+
call_julienne_assert(shape(self%lower_) .equalsExpected. shape(self%upper_))
62+
end procedure
63+
64+
module procedure interpolate_to_faces
65+
integer row
66+
associate( &
67+
N => size(centers) , inner_cols => size(self%inner_) &
68+
,upper_rows => size(self%upper_,1), upper_cols => size(self%upper_,2) &
69+
,lower_rows => size(self%lower_,1), lower_cols => size(self%lower_,2) &
70+
)
71+
faces = [ self%first_ * centers(1) &
72+
, matmul(self%upper_ , centers(1:upper_cols)) &
73+
,[(dot_product(self%inner_ , centers(row-upper_rows+1 : row-upper_rows+inner_cols)), row = upper_rows+1, N-lower_rows-1)] &
74+
, matmul(self%lower_ , centers(N-lower_cols+1:N)) &
75+
, self%final_ * centers(N) &
76+
]
77+
end associate
5378
end procedure
5479

55-
end submodule interpolation_operators_1D_s
80+
end submodule interpolator_1D_s

‎src/formal/tensors_1D_m.F90‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -311,8 +311,8 @@ pure module function reduced_order_boundary_depth(self) result(num_nodes)
311311
integer num_nodes
312312
end function
313313

314-
pure module function div_dyad(self) result(divergence_1D)
315-
!! Result is mimetic divergence of the dyad_1D_t "self"
314+
pure module function div_dyad(self) result(vector_1D)
315+
!! Result is the vector divergence of the dyad_1D_t "self"
316316
implicit none
317317
class(dyad_1D_t), intent(in) :: self
318318
type(vector_1D_t) vector_1D

0 commit comments

Comments
 (0)