Lines Matching +full:format +full:- +full:c
1 C Copyright (c) 2017-2026, Lawrence Livermore National Security, LLC and other CEED contributors
2 C All Rights Reserved. See the top-level COPYRIGHT and NOTICE files for details.
3 C
4 C SPDX-License-Identifier: (BSD-2-Clause)
5 C
6 C This file is part of CEED: http://github.com/ceed
8 C> @file
9 C> Mass and diffusion operators examples using Nek5000
10 C_TESTARGS(name="BP1") -c {ceed_resource} -e bp1 -n 1 -b 4 -test
11 C_TESTARGS(name="BP3") -c {ceed_resource} -e bp3 -n 1 -b 4 -test
13 C-----------------------------------------------------------------------
17 C Set up mass operator
18 C Input: u1,u2,u3,q Output: v1,v2,ierr
30 C Quadrature Point Loop
44 g11 = (a22*a33-a23*a32)
45 g12 = (a13*a32-a33*a12)
46 g13 = (a12*a23-a22*a13)
48 g21 = (a23*a31-a21*a33)
49 g22 = (a11*a33-a31*a13)
50 g23 = (a13*a21-a23*a11)
52 g31 = (a21*a32-a22*a31)
53 g32 = (a12*a31-a32*a11)
54 g33 = (a11*a22-a21*a12)
58 C Rho
61 C RHS
70 C-----------------------------------------------------------------------
74 C Apply mass operator
75 C Input: u1,u2,q Output: v1,ierr
82 C Quadrature Point Loop
89 C-----------------------------------------------------------------------
93 C Set up diffusion operator
94 C Input: u1,u2,u3,q Output: v1,v2,ierr
105 real*8 c(3),k(3)
107 C Quadrature Point Loop
111 c(1)=0.
112 c(2)=1.
113 c(3)=2.
130 g11 = (a22*a33-a23*a32)
131 g12 = (a13*a32-a33*a12)
132 g13 = (a12*a23-a22*a13)
134 g21 = (a23*a31-a21*a33)
135 g22 = (a11*a33-a31*a13)
136 g23 = (a13*a21-a23*a11)
138 g31 = (a21*a32-a22*a31)
139 g32 = (a12*a31-a32*a11)
140 g33 = (a11*a22-a21*a12)
146 C Geometric factors
147 C Stored in Voigt convention
148 C 0 5 4
149 C 5 1 3
150 C 4 3 2
158 C RHS
160 $ *dsin(pi*(c(1)+k(1)*u1(i+0*q)))
161 $ *dsin(pi*(c(2)+k(2)*u1(i+1*q)))
162 $ *dsin(pi*(c(3)+k(3)*u1(i+2*q)))
169 C-----------------------------------------------------------------------
173 C Apply diffusion operator
174 C Input: u1,u2,q Output: v1,ierr
181 C Quadrature Point Loop
193 C-----------------------------------------------------------------------
195 C Set h2 as rhoJac
196 C Input: bmq,nxq Output: h2
217 C-----------------------------------------------------------------------
219 C Set distance initial condition for BP1
220 C Input: Output: e
240 C-----------------------------------------------------------------------
242 C Set sine initial condition for BP3
243 C Input: Output: e
248 real*8 c(3),k(3)
253 c(1)=0.
254 c(2)=1.
255 c(3)=2.
265 e(i) = dsin(pi*(c(1)+k(1)*x))
266 & *dsin(pi*(c(2)+k(2)*y))
267 & *dsin(pi*(c(3)+k(3)*z))
275 C-----------------------------------------------------------------------
281 C e = gllel(eg)
288 C-----------------------------------------------------------------------
290 C
291 C Note: this is an acceleration term, NOT a force!
292 C Thus, ffx will subsequently be multiplied by rho(x,t).
293 C
298 C e = gllel(eg)
306 C-----------------------------------------------------------------------
318 C-----------------------------------------------------------------------
320 C NOTE ::: This subroutine MAY NOT be called by every process
333 C-----------------------------------------------------------------------
347 C-----------------------------------------------------------------------
354 C-----------------------------------------------------------------------
369 C-----------------------------------------------------------------------
376 C-----------------------------------------------------------------------
378 C Apply zero Dirichlet boundary conditions
379 C Input: h2,nel Output: r1
391 C-----------------------------------------------------------------------
393 C Compute Linfty norm of (x-y)
394 C Input: x,y Output: n
402 diff=abs(x(i)-y(i))
415 glrdif = -dmx ! Negative indicates something strange
420 C-----------------------------------------------------------------------
422 C 3D transpose of local gradient
423 C Input: u,N,D,Dt Output: ur,us,ut
439 c-----------------------------------------------------------------------
441 C 3D transpose of local gradient
442 C Input: ur,us,ut,N,D,Dt Output: u
462 C-----------------------------------------------------------------------
464 C Routine to generate elemental geometric matrices on mesh 1
465 C (Gauss-Legendre Lobatto mesh).
469 parameter (lg=3+3*(ldim-2),lzq=lx1+1,lxyd=lzq**ldim)
500 nxqm1 = lzq-1
503 call intp_rstd (tmp,xm1(1,1,1,e),lx1,lzq,if3d,0) ! 0-->Fwd interpolation
525 g11 = (a22*a33-a23*a32)
526 g12 = (a13*a32-a33*a12)
527 g13 = (a12*a23-a22*a13)
529 g21 = (a23*a31-a21*a33)
530 g22 = (a11*a33-a31*a13)
531 g23 = (a13*a21-a23*a11)
533 g31 = (a21*a32-a22*a31)
534 g32 = (a12*a31-a32*a11)
535 g33 = (a11*a22-a21*a12)
554 C-----------------------------------------------------------------------
556 C Generate diagonal preconditioner for Helmholtz operator
557 C Input: h1,h2 Output: d
561 parameter (lxyz=lx1*ly1*lz1,lg=3+3*(ldim-2))
595 C
596 C Add cross terms if element is deformed.
597 C
600 do i2=1,ly1,ly1-1
601 do i1=1,lx1,lx1-1
631 C If axisymmetric, add a diagonal term in the radial direction (ISD=2)
658 1 format(i9,1p4e12.4,' diag prec')
662 C-----------------------------------------------------------------------
664 C Generate dummy diagonal preconditioner for Helmholtz operator
665 C Input: h1,h2 Output: d
676 C-----------------------------------------------------------------------
695 C-----------------------------------------------------------------------
697 C Solution to BP1 using libCEED
705 parameter (lx=lx1*ly1*lz1,lg=3+3*(ldim-2),lq=lzq**ldim)
738 C Set up coordinates
740 do j=0,nelt-1
752 C Init ceed library
758 C Set up Nek geometry data
762 C Set up true soln
766 C Set up solver parameters
767 tol = 1e-10
773 C Create ceed basis for mesh and computation
783 C Create ceed element restrictions for mesh and computation
795 C Create ceed vectors
806 C Create ceed qfunctions for masssetupf and massf
831 C Create ceed operators
854 C Compute setup data
860 C Set up true RHS
863 C Set up algebraic RHS with libCEED
872 C Set up algebraic RHS with Nek5000
878 C Solve true RHS
885 C Output
886 telaps = (tstop-tstart)
895 if (dabs(er1)>5e-3) then
900 nx = nx1-1
905 dofps = nnode/telaps ! DOF/sec - scalar form
912 C Solve libCEED algebraic RHS
920 C Output
921 telaps = (tstop-tstart)
930 if (dabs(er1)>1e-5) then
935 nx = nx1-1
940 dofps = nnode/telaps ! DOF/sec - scalar form
947 C Solve Nek5000 algebraic RHS
955 C Output
956 telaps = (tstop-tstart)
965 if (dabs(er1)>1e-5) then
970 nx = nx1-1
975 dofps = nnode/telaps ! DOF/sec - scalar form
982 1 format(a12,i7,i3,i7,i10,i14,i10,i4,1p4e13.5)
983 3 format(i3,i9,e12.4,1x,a8,i9)
985 C Destroy ceed handles
1004 C-----------------------------------------------------------------------
1006 C Solution to BP3 using libCEED
1014 parameter (lx=lx1*ly1*lz1,lg=3+3*(ldim-2),lq=lzq**ldim)
1048 C Set up coordinates and mask
1050 do j=0,nelt-1
1069 C Init ceed library
1075 C Set up Nek geometry data
1078 C Set up true soln
1083 C Set up solver parameters
1084 tol = 1e-10
1090 C Create ceed basis for mesh and computation
1100 C Create ceed element restrictions for mesh and computation
1114 C Create ceed vectors
1125 C Create ceed qfunctions for diffsetupf and diffusionf
1150 C Create ceed operators
1173 C Compute setup data
1179 C Set up true RHS
1183 C Set up algebraic RHS with libCEED
1193 C Set up algebraic RHS with Nek5000
1200 C Solve true RHS
1207 C Output
1208 telaps = (tstop-tstart)
1217 if (dabs(er1)>1e-3) then
1222 nx = nx1-1
1227 dofps = nnode/telaps ! DOF/sec - scalar form
1234 C Solve libCEED algebraic RHS
1242 C Output
1243 telaps = (tstop-tstart)
1252 if (dabs(er1)>1e-9) then
1257 nx = nx1-1
1262 dofps = nnode/telaps ! DOF/sec - scalar form
1269 C Solve Nek5000 algebraic RHS
1277 C Output
1278 telaps = (tstop-tstart)
1287 if (dabs(er1)>1e-9) then
1292 nx = nx1-1
1297 dofps = nnode/telaps ! DOF/sec - scalar form
1304 1 format(a12,i7,i3,i7,i10,i14,i10,i4,1p4e13.5)
1305 3 format(i3,i9,e12.4,1x,a8,i9)
1307 C Destroy ceed handles
1326 C-----------------------------------------------------------------------
1329 C Scalar conjugate gradient iteration for solution of uncoupled
1330 C Helmholtz equations
1331 C Input: r1,h1,h2,rmult,binv,tin,ceed,ceed_op,vec_p1,vec_ap1,bpname
1332 C Output: u1,maxit
1339 C INPUT: rhs1 - rhs
1340 C h1 - exact solution
1357 nx = nx1-1 ! Polynomial order (just for i/o)
1362 call setprecn_bp3(dpc,h1,h2) ! Set up diagional pre-conidtioner
1364 call setprecn_bp1(dpc,h1,h2) ! Set up diagional pre-conidtioner
1371 s=rmult(i) ! -1
1389 r1(i)=r1(i)-alph(1)*ap1(i)
1392 C tolerance check here
1401 C if (nio.eq.0) write(6,1) ifield,istep,iter,nx,(wv(k),k=1,1)
1402 1 format(i2,i9,i5,i4,1p1e12.4,' cggos')
1410 C if (nio.eq.0) write(6,2) iter,enorm,alph(1),pap(1),'alpha'
1411 2 format(i5,1p3e12.4,2x,a5)
1424 iter = iter-1
1430 3000 format(i12,1x,'cggo scalar:',i6,1p5e13.4)
1431 3001 format(2i6,' Unconverged cggo scalar: rbnorm =',1p2e13.6)
1435 C-----------------------------------------------------------------------
1438 C Vector conjugate gradient matvec for solution of uncoupled
1439 C Helmholtz equations
1440 C Input: pap,p1,h1,h2,bpname,ceed,ceed_op,vec_ap1,vec_p1
1441 C Output: ap1
1446 parameter (lx=lx1*ly1*lz1,lg=3+3*(ldim-2))
1480 C-----------------------------------------------------------------------
1482 C Local matrix-vector for solution of BP3 (stiffness matrix)
1483 C Input: u,g,h1,h2,b,ju,us,ut Output: w
1487 parameter (lxyz=lx1*ly1*lz1,lg=3+3*(ldim-2))
1495 call intp_rstd (ju,u,lx1,nxq,if3d,0) ! 0 --> Fwd interpolation
1499 call intp_rstd (w,ju,lx1,nxq,if3d,1) ! 1 --> ju-->u
1503 C-----------------------------------------------------------------------
1505 C Vector conjugate gradient matvec for solution of BP1 (mass matrix)
1506 C Input: pap,p1,h1,h2 Output: ap1
1510 parameter (lx=lx1*ly1*lz1,lg=3+3*(ldim-2))
1542 C-----------------------------------------------------------------------
1544 C Local matrix-vector for solution of BP3 (stiffness matrix)
1545 C Input: u,g,ur,us,ut,wk Output: w
1557 n = lzq-1
1559 call intp_rstd (wk,u,lx1,lzq,if3d,0) ! 0 --> Fwd interpolation
1572 call intp_rstd (w,wk,lx1,lzq,if3d,1) ! 1 --> ju-->u
1576 C-----------------------------------------------------------------------
1578 C Vector conjugate gradient matvec for solution of BP3 (stiffness matrix)
1579 C Input: pap,p1,h1,h2 Output: ap1
1584 parameter (lx=lx1*ly1*lz1,lg=3+3*(ldim-2),lq=lzq**ldim)
1610 C-----------------------------------------------------------------------
1612 C Vector conjugate gradient matvec for solution of uncoupled
1613 C Helmholtz equations
1614 C Input: pap,p1,h1,h2,bpname Output: ap1
1639 C-----------------------------------------------------------------------
1641 C Get BP to run
1642 C Input: Output: bp
1658 C-----------------------------------------------------------------------
1660 C Get CEED backend specification
1661 C Input: Output: spec
1672 C-----------------------------------------------------------------------
1674 C Get test mode flag
1675 C Input: Output: test
1689 C-----------------------------------------------------------------------