1! Copyright (c) 2017-2026, Lawrence Livermore National Security, LLC and other CEED contributors. 2! All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3 4! SPDX-License-Identifier: BSD-2-Clause 5 6! This file is part of CEED: http:Cgithub.com/ceed 7 8! libCEED Example 1 9 10! This example illustrates a simple usage of libCEED to compute the volume of a 3D body using matrix-free application of a mass operator. 11! Arbitrary mesh and solution degrees in 1D, 2D and 3D are supported from the same code. 12 13! The example has no dependencies, and is designed to be self-contained. 14! For additional examples that use external discretization libraries (MFEM, PETSc, etc.) see the subdirectories in libceed/examples. 15 16! All libCEED objects use a Ceed device object constructed based on a command line argument (-ceed). 17 18! Build with: 19 20! make ex1-volume [CEED_DIR = </path/to/libceed>] 21 22! Sample runs: 23 24! ./ex1-volume-f 25! ./ex1-volume-f -ceed /cpu/self 26! ./ex1-volume-f -ceed /gpu/cuda 27 28! Test in 1D-3D 29! TESTARGS(name = "1D User QFunction") -ceed {ceed_resource} -d 1 -t 30! TESTARGS(name = "2D User QFunction") -ceed {ceed_resource} -d 2 -t 31! TESTARGS(name = "3D User QFunction") -ceed {ceed_resource} -d 3 -t 32! TESTARGS(name = "1D Gallery QFunction") -ceed {ceed_resource} -d 1 -t -g 33! TESTARGS(name = "2D Gallery QFunction") -ceed {ceed_resource} -d 2 -t -g 34! TESTARGS(name = "3D Gallery QFunction") -ceed {ceed_resource} -d 3 -t -g 35 36!> @file 37!> libCEED example using mass operator to compute volume 38 39 include 'ex1-volume-f.h' 40 41!----------------------------------------------------------------------- 42subroutine getcartesianmeshsize(fe_dim, degree, prob_size, num_xyz) 43 implicit none 44 integer fe_dim 45 integer degree 46 integer prob_size 47 integer num_xyz(3) 48 49 integer num_elem 50 integer s, r, d, sd 51 num_elem = prob_size/(degree**fe_dim) 52 s = 0 53 54! Use the approximate formula: 55! prob_size ~ num_elem * degree^dim 56! find s: num_elem/2 < 2^s <= num_elem 57 58 do while (num_elem > 1) 59 num_elem = num_elem/2 60 s = s + 1 61 end do 62 r = mod(s, fe_dim) 63 64 do d = 1, fe_dim 65 sd = s/fe_dim 66 if (r > 0) then 67 sd = sd + 1 68 r = r - 1 69 end if 70 num_xyz(d) = ISHFT(1, sd) 71 end do 72end 73 74!----------------------------------------------------------------------- 75subroutine buildcartesianrestriction(ceed, fe_dim, num_xyz, degree, num_comp, mesh_size, num_qpts, restriction,& 76& q_data_restriction, err) 77 implicit none 78 include 'ceed/fortran.h' 79 80 integer ceed 81 integer fe_dim 82 integer num_xyz(3) 83 integer degree 84 integer num_comp 85 integer mesh_size 86 integer num_qpts 87 integer restriction 88 integer q_data_restriction 89 integer err 90 91 integer p 92 integer num_nodes 93 integer elem_qpts 94 integer num_elem 95 integer scalar_size 96 integer nd(3) 97 integer elem_nodes_size 98 integer e_xyz(3), re 99 integer g_nodes, g_nodes_stride, r_nodes 100 integer, dimension (:), allocatable :: elem_nodes 101 102 integer i, j, k 103 104 p = degree + 1 105 num_nodes = p**fe_dim 106 elem_qpts = num_qpts**fe_dim 107 num_elem = 1 108 scalar_size = 1 109 110 do i = 1, fe_dim 111 num_elem = num_elem * num_xyz(i) 112 nd(i) = num_xyz(i) * (p - 1) + 1 113 scalar_size = scalar_size*nd(i) 114 end do 115 mesh_size = scalar_size*num_comp 116! elem: 0 1 n-1 117! |---*-...-*---|---*-...-*---|- ... -|--...--| 118! num_nodes: 0 1 p-1 p p+1 2*p n*p 119 elem_nodes_size = num_elem*num_nodes 120 allocate (elem_nodes(elem_nodes_size)) 121 122 do i = 1, num_elem 123 e_xyz(1) = 1 124 e_xyz(2) = 1 125 e_xyz(3) = 1 126 re = i - 1 127 128 do j = 1, fe_dim 129 e_xyz(j) = mod(re, num_xyz(j)) 130 re = re/num_xyz(j) 131 end do 132 133 do j = 1, num_nodes 134 g_nodes = 0 135 g_nodes_stride = 1 136 r_nodes = j - 1 137 138 do k = 1, fe_dim 139 g_nodes = g_nodes + (e_xyz(k) * (p - 1) + mod(r_nodes, p)) * g_nodes_stride 140 g_nodes_stride = g_nodes_stride * nd(k) 141 r_nodes = r_nodes/p 142 end do 143 elem_nodes((i - 1) * num_nodes + j) = g_nodes 144 end do 145 end do 146 147 call ceedelemrestrictioncreate(ceed, num_elem, num_nodes, num_comp, scalar_size, mesh_size, ceed_mem_host,& 148 &ceed_copy_values, elem_nodes, restriction, err) 149 if (q_data_restriction /= ceed_qfunction_none) then 150 call ceedelemrestrictioncreatestrided(ceed, num_elem, elem_qpts, num_comp, num_comp * elem_qpts * num_elem,& 151 &ceed_strides_backend, q_data_restriction, err) 152 end if 153 deallocate (elem_nodes) 154end 155 156!----------------------------------------------------------------------- 157subroutine transformmeshcoords(fe_dim, mesh_size, coords, exact_volume, err) 158 implicit none 159 160 integer fe_dim 161 integer mesh_size, scalar_size 162 real*8 coords(mesh_size) 163 real*8 exact_volume 164 real*8 m_pi, m_pi_2 165 parameter(m_pi = 3.14159265358979323846d0) 166 parameter(m_pi_2 = 1.57079632679489661923d0) 167 integer err 168 169 integer i 170 real*8 u, v 171 172 scalar_size = mesh_size/fe_dim 173 select case (fe_dim) 174 case (1) 175 do i = 1, scalar_size 176 coords(i) = 0.5d0 + (1.d0/sqrt(3.d0)) * sin((2.d0/3.d0) * m_pi * (coords(i) - 0.5d0)) 177 end do 178 exact_volume = 1.d0 179 180 case (2, 3) 181 do i = 1, scalar_size 182 u = 1.d0 + coords(i) 183 v = m_pi_2 * coords(i + scalar_size) 184 185 coords(i) = u * cos(v) 186 coords(i + scalar_size) = u * sin(v) 187 end do 188 exact_volume = 3.d0/4.d0 * m_pi 189 end select 190end 191 192!----------------------------------------------------------------------- 193subroutine setcartesianmeshcoords(fe_dim, num_xyz, mesh_degree, mesh_coords, exact_volume, err) 194 implicit none 195 include 'ceed/fortran.h' 196 197 integer fe_dim 198 integer num_xyz(3) 199 integer mesh_degree 200 integer mesh_coords 201 real*8 exact_volume 202 integer err 203 204 integer p 205 integer scalar_size 206 integer coords_size 207 integer r_nodes 208 integer d_1d 209 integer nd(3) 210 real*8, dimension (:), allocatable :: nodes, qpts 211 real*8, dimension (:), allocatable :: coords 212 integer*8 offset 213 integer i, j 214 p = mesh_degree + 1 215 scalar_size = 1 216 217 do i = 1, fe_dim 218 nd(i) = num_xyz(i) * (p - 1) + 1 219 scalar_size = scalar_size * nd(i) 220 end do 221 222 coords_size = scalar_size * fe_dim 223 allocate (coords(coords_size)) 224 225! The H1 basis uses Lobatto quadrature points as nodes 226 allocate (nodes(p)) 227 allocate (qpts(p)) 228 call ceedlobattoquadrature(p, nodes, qpts, err) 229 deallocate(qpts) 230 do i = 1, p 231 nodes(i) = 0.5 + 0.5 * nodes(i) 232 end do 233 234 do i = 1, scalar_size 235 r_nodes = i - 1 236 237 do j = 1, fe_dim 238 d_1d = mod(r_nodes, nd(j)) 239 coords(scalar_size * (j - 1) + i) = ((d_1d/(p - 1)) + nodes(mod(d_1d, p - 1) + 1))/num_xyz(j) 240 r_nodes = r_nodes/nd(j) 241 end do 242 end do 243 deallocate(nodes) 244 245 call transformmeshcoords(fe_dim, coords_size, coords, exact_volume, err) 246 247 offset = 0 248 call ceedvectorsetarray(mesh_coords, ceed_mem_host, ceed_copy_values, coords, offset, err) 249 deallocate(coords) 250end 251 252!----------------------------------------------------------------------- 253program main 254 implicit none 255 include 'ceed/fortran.h' 256 257 character ceed_spec*32 258 integer fe_dim, num_comp_x, mesh_degree, sol_degree, num_qpts 259 integer num_elem, num_xyz(3), elem_qpts 260 integer prob_size, mesh_size, sol_size 261 integer help, test, gallery, benchmark 262 integer i, num_args, err 263 character arg*32, arg_value*32 264 real*8 exact_volume, computed_volume 265 266 integer ceed 267 real*8, dimension (:), allocatable :: u_array, v_array 268 integer mesh_coords, q_data, u, v 269 integer mesh_restriction, sol_restriction, q_data_restriction 270 integer mesh_basis, sol_basis 271 integer*8 offset 272 integer build_ctx 273 integer build_ctx_size 274 parameter(build_ctx_size = 2) 275 integer*8 build_ctx_data(build_ctx_size) 276 integer qf_build, qf_apply 277 integer op_build, op_apply 278 279 external build_mass, apply_mass 280 281! Initial values 282 ceed_spec = '/cpu/self' 283 fe_dim = 3 284 num_comp_x = 3 285 mesh_degree = 4 286 sol_degree = 4 287 num_qpts = mesh_degree + 2 288 prob_size = -1 289 help = 0 290 test = 0 291 gallery = 0 292 benchmark = 0 293 294! Process command line arguments 295 296 num_args = command_argument_count() 297 do i = 1, num_args 298 call get_command_argument(i, arg) 299 300 select case (arg) 301! LCOV_EXCL_START 302 case ('-h') 303 help = 1 304 305 case ('-c', '-ceed') 306 call get_command_argument(i + 1, ceed_spec) 307 308 case ('-d') 309 call get_command_argument(i + 1, arg_value) 310 read(arg_value, '(I10)') fe_dim 311 num_comp_x = fe_dim 312 313 case ('-m') 314 call get_command_argument(i + 1, arg_value) 315 read(arg_value, '(I10)') mesh_degree 316 317 case ('-p') 318 call get_command_argument(i + 1, arg_value) 319 read(arg_value, '(I10)') sol_degree 320 321 case ('-q') 322 call get_command_argument(i + 1, arg_value) 323 read(arg_value, '(I10)') num_qpts 324 325 case ('-s') 326 call get_command_argument(i + 1, arg_value) 327 read(arg_value, '(I10)') prob_size 328 329 case ('-b') 330 call get_command_argument(i + 1, arg_value) 331 read(arg_value, '(I10)') benchmark 332! LCOV_EXCL_STOP 333 334 case ('-t') 335 test = 1 336 337 case ('-g') 338 gallery = 1 339 end select 340 end do 341 342 if (prob_size < 0) then 343 if (test == 1) then 344 prob_size = 8 * 16 345 else 346 prob_size = 256 * 1024 347 end if 348 end if 349 350! Print options 351 if ((test /= 1) .OR. (help == 1)) then 352! LCOV_EXCL_START 353 write (*, *) 'Selected options: [command line option] : <current value>' 354 write (*, *) ' Ceed specification [-c] : ', ceed_spec 355 write (*, *) ' Mesh dimension [-d] : ', fe_dim 356 write (*, *) ' Mesh degree [-m] : ', mesh_degree 357 write (*, *) ' Solution degree [-p] : ', sol_degree 358 write (*, *) ' Num. 1D quadrature pts [-q] : ', num_qpts 359 write (*, *) ' Approx. # unknowns [-s] : ', prob_size 360 if (gallery == 1) then 361 write (*, *) ' QFunction source [-g] : gallery' 362 else 363 write (*, *) ' QFunction source [-g] : header' 364 end if 365 if (help == 1) then 366 if (test == 0) then 367 write (*, *) 'Test/quiet mode is OFF (use -t to enable)' 368 else 369 write (*, *) 'Test/quiet mode is ON' 370 end if 371 end if 372! LCOV_EXCL_STOP 373 end if 374 375! Select appropriate backend and logical device based on the (-ceed) command line argument 376 call ceedinit(trim(ceed_spec)//char(0), ceed, err) 377 378! Construct the mesh and solution bases 379 call ceedbasiscreatetensorh1lagrange(ceed, fe_dim, num_comp_x, mesh_degree + 1, num_qpts, ceed_gauss, mesh_basis,& 380 &err) 381 call ceedbasiscreatetensorh1lagrange(ceed, fe_dim, 1, sol_degree + 1, num_qpts, ceed_gauss, sol_basis, err) 382 383! Determine the mesh size based on the given approximate problem size 384 call getcartesianmeshsize(fe_dim, sol_degree, prob_size, num_xyz) 385 if (test == 0) then 386! LCOV_EXCL_START 387 write (*, '(A16, I8)', advance='no') 'Mesh size: nx = ', num_xyz(1) 388 if (num_comp_x > 1) then 389 write (*, '(A7, I8)', advance='no') ', ny = ', num_xyz(2) 390 end if 391 if (num_comp_x > 2) then 392 write (*, '(A7, I8)', advance='no') ', nz = ', num_xyz(3) 393 end if 394 write (*, *) 395! LCOV_EXCL_STOP 396 endif 397 398! Build CeedElemRestriction objects describing the mesh and solution discrete representation 399 call buildcartesianrestriction(ceed, fe_dim, num_xyz, mesh_degree, num_comp_x, mesh_size, num_qpts,& 400 &mesh_restriction, ceed_qfunction_none, err) 401 call buildcartesianrestriction(ceed, fe_dim, num_xyz, sol_degree, 1, sol_size, num_qpts, sol_restriction,& 402 &q_data_restriction, err) 403 404 if (test == 0) then 405! LCOV_EXCL_START 406 write (*, *) 'Number of mesh nodes : ', mesh_size/fe_dim 407 write (*, *) 'Number of solution nodes : ', sol_size 408! LCOV_EXCL_STOP 409 end if 410 411! Create a CeedVector with the mesh coordinates 412! Apply a transformation to the mesh 413 call ceedvectorcreate(ceed, mesh_size, mesh_coords, err) 414 call setcartesianmeshcoords(fe_dim, num_xyz, mesh_degree, mesh_coords, exact_volume, err) 415 416! Context data to be passed to the 'build_mass' QFunction 417 build_ctx_data(1) = fe_dim 418 build_ctx_data(2) = num_comp_x 419 call ceedqfunctioncontextcreate(ceed, build_ctx, err) 420! Note: The context technically only takes arrays of double precision values, but we can pass arrays of ints of the same length 421 offset = 0 422 call ceedqfunctioncontextsetdata(build_ctx, ceed_mem_host, ceed_use_pointer, build_ctx_size, build_ctx_data,& 423 &offset, err) 424 425! Create the QFunction that builds the mass operator (i.e. computes its quadrature data) and set its context data 426 if (gallery == 1) then 427 select case (fe_dim) 428 case (1) 429 call ceedqfunctioncreateinteriorbyname(ceed, 'Mass1DBuild', qf_build, err) 430 431 case (2) 432 call ceedqfunctioncreateinteriorbyname(ceed, 'Mass2DBuild', qf_build, err) 433 434 case (3) 435 call ceedqfunctioncreateinteriorbyname(ceed, 'Mass3DBuild', qf_build, err) 436 end select 437 else 438 call ceedqfunctioncreateinterior(ceed, 1, build_mass,& 439 &SOURCE_DIR& 440 &//'ex1-volume-f-c.h:build_mass'//char(0), qf_build, err) 441 call ceedqfunctionaddinput(qf_build, 'dx', num_comp_x * fe_dim, ceed_eval_grad, err) 442 call ceedqfunctionaddinput(qf_build, 'weights', 1, ceed_eval_weight, err) 443 call ceedqfunctionaddoutput(qf_build, 'qdata', 1, ceed_eval_none, err) 444 call ceedqfunctionsetcontext(qf_build, build_ctx, err) 445 end if 446 447! Create the operator that builds the quadrature data for the mass operator 448 call ceedoperatorcreate(ceed, qf_build, ceed_qfunction_none, ceed_qfunction_none, op_build, err) 449 call ceedoperatorsetfield(op_build, 'dx', mesh_restriction, mesh_basis, ceed_vector_active, err) 450 call ceedoperatorsetfield(op_build, 'weights', ceed_elemrestriction_none, mesh_basis, ceed_vector_none, err) 451 call ceedoperatorsetfield(op_build, 'qdata', q_data_restriction, ceed_basis_none, ceed_vector_active, err) 452 453! Compute the quadrature data for the mass operator 454 num_elem = 1 455 elem_qpts = num_qpts**fe_dim 456 do i = 1, fe_dim 457 num_elem = num_elem * num_xyz(i) 458 end do 459 call ceedvectorcreate(ceed, num_elem * elem_qpts, q_data, err) 460 call ceedoperatorapply(op_build, mesh_coords, q_data, ceed_request_immediate, err) 461 462! Create the QFunction that defines the action of the mass operator 463 if (gallery == 1) then 464 call ceedqfunctioncreateinteriorbyname(ceed, 'MassApply', qf_apply, err) 465 else 466 call ceedqfunctioncreateinterior(ceed, 1, apply_mass,& 467 &SOURCE_DIR& 468 &//'ex1-volume-f-c.h:apply_mass'//char(0), qf_apply, err) 469 call ceedqfunctionaddinput(qf_apply, 'u', 1, ceed_eval_interp, err) 470 call ceedqfunctionaddinput(qf_apply, 'qdata', 1, ceed_eval_none, err) 471 call ceedqfunctionaddoutput(qf_apply, 'v', 1, ceed_eval_interp, err) 472 end if 473 474! Create the mass operator 475 call ceedoperatorcreate(ceed, qf_apply, ceed_qfunction_none, ceed_qfunction_none, op_apply, err) 476 call ceedoperatorsetfield(op_apply, 'u', sol_restriction, sol_basis, ceed_vector_active, err) 477 call ceedoperatorsetfield(op_apply, 'qdata', q_data_restriction, ceed_basis_none, q_data, err) 478 call ceedoperatorsetfield(op_apply, 'v', sol_restriction, sol_basis, ceed_vector_active, err) 479 480! Create auxiliary solution-size vectors 481 allocate (u_array(sol_size)) 482 allocate (v_array(sol_size)) 483 484 call ceedvectorcreate(ceed, sol_size, u, err) 485 offset = 0 486 call ceedvectorsetarray(u, ceed_mem_host, ceed_use_pointer, u_array, offset, err) 487 call ceedvectorcreate(ceed, sol_size, v, err) 488 offset = 0 489 call ceedvectorsetarray(v, ceed_mem_host, ceed_use_pointer, v_array, offset, err) 490 491! Initialize 'u' with ones 492 call ceedvectorsetvalue(u, 1.d0, err) 493 494! Compute the mesh volume using the mass operator: volume = 1^T \cdot M \cdot 1 495 call ceedoperatorapply(op_apply, u, v, ceed_request_immediate, err) 496 497! Benchmark runs 498 if (test /= 1 .AND. benchmark /= 0) then 499! LCOV_EXCL_START 500 write (*, *) ' Executing ', benchmark, ' benchmarking runs...' 501! LCOV_EXCL_STOP 502 end if 503 do i = 1, benchmark 504! LCOV_EXCL_START 505 call ceedoperatorapply(op_apply, u, v, ceed_request_immediate, err) 506! LCOV_EXCL_STOP 507 end do 508 509! Compute and print the sum of the entries of 'v' giving the mesh volume 510 computed_volume = 0.d0 511 512 call ceedvectorgetarrayread(v, ceed_mem_host, v_array, offset, err) 513 do i = 1, sol_size 514 computed_volume = computed_volume + v_array(offset + i) 515 end do 516 call ceedvectorrestorearrayread(v, v_array, offset, err) 517 518 if (test /= 1) then 519! LCOV_EXCL_START 520 write (*, *) ' done.' 521 write (*, *) 'Exact mesh volume :', exact_volume 522 write (*, *) 'Computed mesh volume :', computed_volume 523 write (*, *) 'Volume error :', (exact_volume - computed_volume) 524! LCOV_EXCL_STOP 525 else 526 if (fe_dim == 1) then 527 if (abs(exact_volume - computed_volume) > 200.d0 * 1e-15) then 528! LCOV_EXCL_START 529 write (*, *) 'Volume error : ', (exact_volume - computed_volume) 530! LCOV_EXCL_STOP 531 end if 532 else 533 if (abs(exact_volume - computed_volume) > 1e-5) then 534! LCOV_EXCL_START 535 write (*, *) 'Volume error : ', (exact_volume - computed_volume) 536! LCOV_EXCL_STOP 537 end if 538 end if 539 end if 540 541! Free dynamically allocated memory 542 call ceedvectordestroy(mesh_coords, err) 543 call ceedvectordestroy(q_data, err) 544 call ceedvectordestroy(u, err) 545 call ceedvectordestroy(v, err) 546 deallocate (u_array) 547 deallocate (v_array) 548 call ceedbasisdestroy(sol_basis, err) 549 call ceedbasisdestroy(mesh_basis, err) 550 call ceedqfunctioncontextdestroy(build_ctx, err) 551 call ceedqfunctiondestroy(qf_build, err) 552 call ceedqfunctiondestroy(qf_apply, err) 553 call ceedoperatordestroy(op_build, err) 554 call ceedoperatordestroy(op_apply, err) 555 call ceeddestroy(ceed, err) 556end 557!----------------------------------------------------------------------- 558