formal/example/print-assembled-1D-operators.F90 at main · BerkeleyLab/formal · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
! Copyright (c) 2026, The Regents of the University of California
! Terms of use are as specified in LICENSE.txt

program print_assembled_1D_operators
  !! Print the full 2nd- and 4th-order mimetic gradient and divergence
  !! operator matrices as comma-separated rows for 1D grids with the
  !! minimum number of cells (16) required for computing the 4th-order,
  !! gradient quadrature weight matrix Q of Corbino and Castillo (2022).
  use julienne_m, only : operator(.csv.), string_t, command_line_t
  use formal_m, only : gradient_operator_1D_t, divergence_operator_1D_t
  implicit none

  type(command_line_t) command_line
  integer row

  command_line_settings: &
  associate( &
     gradient   => command_line%argument_present(["--grad" ]) &
    ,divergence => command_line%argument_present(["--div"  ]) &
    ,order      => command_line%flag_value("--order") &
  )

    if (command_line%argument_present([character(len=len("--help")) :: ("--help"), "-h"])) then
      stop                                 new_line('') // new_line('') &
        // 'Usage:'                                     // new_line('') &
        // '  fpm run \'                                // new_line('') &
        // '  --example print-assembled-1D-operators \' // new_line('') &
        // '  --compiler flang-new \'                   // new_line('') &
        // '  --flag "-O3" \'                           // new_line('') &
        // '  -- [--help|-h] | [--grad]  [--div] [--order <integer>]' // new_line('') // new_line('') &
        // 'where square brackets indicate optional arguments and angular brackets indicate user input values.' // new_line('')
    end if

    default_usage: &
    associate(print_all => .not. any([gradient, divergence, len(order)/=0]))

      if (print_all .or. (gradient   .and. len(order)==0) .or. (gradient   .and. order=="2")) call print_gradient_operator(  k=2, dx=1D0, m=16)
      if (print_all .or. (divergence .and. len(order)==0) .or. (divergence .and. order=="2")) call print_divergence_operator(k=2, dx=1D0, m=16)
      if (print_all .or. (gradient   .and. len(order)==0) .or. (gradient   .and. order=="4")) call print_gradient_operator(  k=4, dx=1D0, m=16)
      if (print_all .or. (divergence .and. len(order)==0) .or. (divergence .and. order=="4")) call print_divergence_operator(k=4, dx=1D0, m=16)

    end associate default_usage
#ifdef __GFORTRAN__
    stop
#endif
  end associate command_line_settings

contains

  subroutine print_gradient_operator(k, dx, m)
    integer, intent(in) :: k, m
    double precision, intent(in) :: dx

    print *, new_line(""), "Gradient operator: order = ", k, " | cells = ", m, " | dx = ", dx

    associate(grad_op => gradient_operator_1D_t(k, dx, cells=m))
      associate(G => grad_op%assemble())
        do row = 1, size(G,1)
          associate(csv_row => .csv. string_t(G(row,:)))
            print '(a)', csv_row%string()
          end associate
        end do
      end associate
    end associate

  end subroutine

  subroutine print_divergence_operator(k, dx, m)
    integer, intent(in) :: k, m
    double precision, intent(in) :: dx

    print *, new_line(""), "Divergence operator: order = ", k, " | cells = ", m, " | dx = ", dx

    associate(div_op => divergence_operator_1D_t(k, dx, cells=m))
      associate(D => div_op%assemble())
        do row = 1, size(D,1)
          associate(csv_row => .csv. string_t(D(row,:)))
            print '(a)', csv_row%string()
          end associate
        end do
      end associate
    end associate

  end subroutine

end program