1!----------------------------------------------------------------------- 2! 3! Header with common subroutine 4! 5 include 't320-basis-f.h' 6!----------------------------------------------------------------------- 7 subroutine setup(ctx,q,u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,u14,& 8& u15,u16,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,v16,ierr) 9 real*8 ctx 10 real*8 u1(1) 11 real*8 u2(1) 12 real*8 v1(1) 13 real*8 w 14 integer q,ierr 15 16 do i=1,q 17 w = u1(i)/(u2(i+0*q)*u2(i+3*q)-u2(i+1*q)*u2(i+2*q)) 18 v1(i+0*q)=w*(u2(i+2*q)*u2(i+2*q)+u2(i+3*q)*u2(i+3*q)) 19 v1(i+1*q)=w*(u2(i+0*q)*u2(i+0*q)+u2(i+1*q)*u2(i+1*q)) 20 v1(i+2*q)=w*(u2(i+0*q)*u2(i+2*q)+u2(i+1*q)*u2(i+3*q)) 21 enddo 22 23 ierr=0 24 end 25!----------------------------------------------------------------------- 26 subroutine diff(ctx,q,u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,u14,& 27& u15,u16,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,v16,ierr) 28 real*8 ctx 29 real*8 u1(1) 30 real*8 u2(1) 31 real*8 v1(1) 32 integer q,ierr 33 34 do i=1,q 35 v1(i+0*q)=u1(i+0*q)*u2(i+0*q)+u1(i+2*q)*u2(i+1*q) 36 v1(i+1*q)=u1(i+2*q)*u2(i+0*q)+u1(i+1*q)*u2(i+1*q) 37 enddo 38 39 ierr=0 40 end 41!----------------------------------------------------------------------- 42 program test 43 44 include 'ceedf.h' 45 46 integer ceed,err,i,j,k 47 integer imode 48 parameter(imode=ceed_noninterlaced) 49 integer stridesxtet(3),stridesqdtet(3),stridesxhex(3),stridesqdhex(3) 50 integer erestrictxtet,erestrictutet,erestrictxitet,erestrictqditet,& 51& erestrictxhex,erestrictuhex,erestrictxihex,erestrictqdihex 52 integer bxtet,butet,bxhex,buhex 53 integer qf_setuptet,qf_difftet,qf_setuphex,qf_diffhex 54 integer op_setuptet,op_difftet,op_setuphex,op_diffhex,op_setup,op_diff 55 integer qdatatet,qdatahex,x,u,v 56 integer nelemtet,nelemhex,ptet,phex,qtet,qhex,d 57 integer row,col,offset 58 parameter(nelemtet=6) 59 parameter(ptet=6) 60 parameter(qtet=4) 61 parameter(nelemhex=6) 62 parameter(phex=3) 63 parameter(qhex=4) 64 parameter(d=2) 65 integer ndofs,nqptstet,nqptshex,nqpts,nx,ny,nxtet,nytet,nxhex 66 parameter(nx=3) 67 parameter(ny=3) 68 parameter(nxtet=3) 69 parameter(nytet=1) 70 parameter(nxhex=3) 71 parameter(ndofs=(nx*2+1)*(ny*2+1)) 72 parameter(nqptstet=nelemtet*qtet) 73 parameter(nqptshex=nelemhex*qhex*qhex) 74 parameter(nqpts=nqptstet+nqptshex) 75 integer indxtet(nelemtet*ptet),indxhex(nelemhex*phex*phex) 76 real*8 arrx(d*ndofs) 77 integer*8 voffset,xoffset 78 79 real*8 qref(d*qtet) 80 real*8 qweight(qtet) 81 real*8 interp(ptet*qtet) 82 real*8 grad(d*ptet*qtet) 83 84 real*8 hv(ndofs) 85 real*8 total 86 87 character arg*32 88 89 external setup,diff 90 91 call getarg(1,arg) 92 93 call ceedinit(trim(arg)//char(0),ceed,err) 94 95! DoF Coordinates 96 do i=0,ny*2 97 do j=0,nx*2 98 arrx(i+j*(ny*2+1)+0*ndofs+1)=1.d0*i/(2*ny) 99 arrx(i+j*(ny*2+1)+1*ndofs+1)=1.d0*j/(2*nx) 100 enddo 101 enddo 102 103 call ceedvectorcreate(ceed,d*ndofs,x,err) 104 xoffset=0 105 call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,arrx,xoffset,err) 106 107! Qdata Vectors 108 call ceedvectorcreate(ceed,nqptstet*d*(d+1)/2,qdatatet,err) 109 call ceedvectorcreate(ceed,nqptshex*d*(d+1)/2,qdatahex,err) 110 111! Tet Elements 112 do i=0,2 113 col=mod(i,nx) 114 row=i/nx 115 offset=col*2+row*(nx*2+1)*2 116 117 indxtet(i*2*ptet+1)=2+offset 118 indxtet(i*2*ptet+2)=9+offset 119 indxtet(i*2*ptet+3)=16+offset 120 indxtet(i*2*ptet+4)=1+offset 121 indxtet(i*2*ptet+5)=8+offset 122 indxtet(i*2*ptet+6)=0+offset 123 124 indxtet(i*2*ptet+7)=14+offset 125 indxtet(i*2*ptet+8)=7+offset 126 indxtet(i*2*ptet+9)=0+offset 127 indxtet(i*2*ptet+10)=15+offset 128 indxtet(i*2*ptet+11)=8+offset 129 indxtet(i*2*ptet+12)=16+offset 130 enddo 131 132! -- Restrictions 133 call ceedelemrestrictioncreate(ceed,imode,nelemtet,ptet,ndofs,d,& 134 & ceed_mem_host,ceed_use_pointer,indxtet,erestrictxtet,err) 135 stridesxtet=[1,ptet,d*ptet] 136 call ceedelemrestrictioncreatestrided(ceed,nelemtet,ptet,& 137 & nelemtet*ptet,d,stridesxtet,erestrictxitet,err) 138 139 call ceedelemrestrictioncreate(ceed,imode,nelemtet,ptet,ndofs,1,& 140 & ceed_mem_host,ceed_use_pointer,indxtet,erestrictutet,err) 141 stridesqdtet=[1,qtet,qtet*d*(d+1)/2] 142 call ceedelemrestrictioncreatestrided(ceed,nelemtet,qtet,nqptstet,& 143 & d*(d+1)/2,stridesqdtet,erestrictqditet,err) 144 145! -- Bases 146 call buildmats(qref,qweight,interp,grad) 147 call ceedbasiscreateh1(ceed,ceed_triangle,d,ptet,qtet,interp,grad,qref,& 148 & qweight,bxtet,err) 149 call buildmats(qref,qweight,interp,grad) 150 call ceedbasiscreateh1(ceed,ceed_triangle,1,ptet,qtet,interp,grad,qref,& 151 & qweight,butet,err) 152 153! -- QFunctions 154 call ceedqfunctioncreateinterior(ceed,1,setup,& 155 &SOURCE_DIR& 156 &//'t522-operator.h:setup'//char(0),qf_setuptet,err) 157 call ceedqfunctionaddinput(qf_setuptet,'_weight',1,ceed_eval_weight,err) 158 call ceedqfunctionaddinput(qf_setuptet,'dx',d*d,ceed_eval_grad,err) 159 call ceedqfunctionaddoutput(qf_setuptet,'rho',d*(d+1)/2,ceed_eval_none,& 160 & err) 161 162 call ceedqfunctioncreateinterior(ceed,1,diff,& 163 &SOURCE_DIR& 164 &//'t522-operator.h:diff'//char(0),qf_difftet,err) 165 call ceedqfunctionaddinput(qf_difftet,'rho',d*(d+1)/2,ceed_eval_none,err) 166 call ceedqfunctionaddinput(qf_difftet,'u',d,ceed_eval_grad,err) 167 call ceedqfunctionaddoutput(qf_difftet,'v',d,ceed_eval_grad,err) 168 169! -- Operators 170! ---- Setup Tet 171 call ceedoperatorcreate(ceed,qf_setuptet,ceed_qfunction_none,& 172 & ceed_qfunction_none,op_setuptet,err) 173 call ceedoperatorsetfield(op_setuptet,'_weight',erestrictxitet,& 174 & bxtet,ceed_vector_none,err) 175 call ceedoperatorsetfield(op_setuptet,'dx',erestrictxtet,& 176 & bxtet,ceed_vector_active,err) 177 call ceedoperatorsetfield(op_setuptet,'rho',erestrictqditet,& 178 & ceed_basis_collocated,qdatatet,err) 179! ---- diff Tet 180 call ceedoperatorcreate(ceed,qf_difftet,ceed_qfunction_none,& 181 & ceed_qfunction_none,op_difftet,err) 182 call ceedoperatorsetfield(op_difftet,'rho',erestrictqditet,& 183 & ceed_basis_collocated,qdatatet,err) 184 call ceedoperatorsetfield(op_difftet,'u',erestrictutet,& 185 & butet,ceed_vector_active,err) 186 call ceedoperatorsetfield(op_difftet,'v',erestrictutet,& 187 & butet,ceed_vector_active,err) 188 189! Hex Elements 190 do i=0,nelemhex-1 191 col=mod(i,nx) 192 row=i/nx 193 offset=(nxtet*2+1)*(nytet*2)*(1+row)+col*2 194 do j=0,phex-1 195 do k=0,phex-1 196 indxhex(phex*(phex*i+k)+j+1)=offset+k*(nxhex*2+1)+j 197 enddo 198 enddo 199 enddo 200 201! -- Restrictions 202 call ceedelemrestrictioncreate(ceed,imode,nelemhex,phex*phex,ndofs,d,& 203 & ceed_mem_host,ceed_use_pointer,indxhex,erestrictxhex,err) 204 stridesxhex=[1,phex*phex,phex*phex*d] 205 call ceedelemrestrictioncreatestrided(ceed,nelemhex,phex*phex,& 206 & nelemhex*phex*phex,d,stridesxhex,erestrictxihex,err) 207 208 call ceedelemrestrictioncreate(ceed,imode,nelemhex,phex*phex,ndofs,1,& 209 & ceed_mem_host,ceed_use_pointer,indxhex,erestrictuhex,err) 210 stridesqdhex=[1,qhex*qhex,qhex*qhex*d*(d+1)/2] 211 call ceedelemrestrictioncreatestrided(ceed,nelemhex,qhex*qhex,& 212 & nqptshex,d*(d+1)/2,stridesqdhex,erestrictqdihex,err) 213 214! -- Bases 215 call ceedbasiscreatetensorh1lagrange(ceed,d,d,phex,qhex,ceed_gauss,& 216 & bxhex,err) 217 call ceedbasiscreatetensorh1lagrange(ceed,d,1,phex,qhex,ceed_gauss,& 218 & buhex,err) 219 220! -- QFunctions 221 call ceedqfunctioncreateinterior(ceed,1,setup,& 222 &SOURCE_DIR& 223 &//'t522-operator.h:setup'//char(0),qf_setuphex,err) 224 call ceedqfunctionaddinput(qf_setuphex,'_weight',1,ceed_eval_weight,err) 225 call ceedqfunctionaddinput(qf_setuphex,'dx',d*d,ceed_eval_grad,err) 226 call ceedqfunctionaddoutput(qf_setuphex,'rho',d*(d+1)/2,ceed_eval_none,& 227 & err) 228 229 call ceedqfunctioncreateinterior(ceed,1,diff,& 230 &SOURCE_DIR& 231 &//'t522-operator.h:diff'//char(0),qf_diffhex,err) 232 call ceedqfunctionaddinput(qf_diffhex,'rho',d*(d+1)/2,ceed_eval_none,err) 233 call ceedqfunctionaddinput(qf_diffhex,'u',d,ceed_eval_grad,err) 234 call ceedqfunctionaddoutput(qf_diffhex,'v',d,ceed_eval_grad,err) 235 236! -- Operators 237! ---- Setup Hex 238 call ceedoperatorcreate(ceed,qf_setuphex,ceed_qfunction_none,& 239 & ceed_qfunction_none,op_setuphex,err) 240 call ceedoperatorsetfield(op_setuphex,'_weight',erestrictxihex,& 241 & bxhex,ceed_vector_none,err) 242 call ceedoperatorsetfield(op_setuphex,'dx',erestrictxhex,& 243 & bxhex,ceed_vector_active,err) 244 call ceedoperatorsetfield(op_setuphex,'rho',erestrictqdihex,& 245 & ceed_basis_collocated,qdatahex,err) 246! ---- diff Hex 247 call ceedoperatorcreate(ceed,qf_diffhex,ceed_qfunction_none,& 248 & ceed_qfunction_none,op_diffhex,err) 249 call ceedoperatorsetfield(op_diffhex,'rho',erestrictqdihex,& 250 & ceed_basis_collocated,qdatahex,err) 251 call ceedoperatorsetfield(op_diffhex,'u',erestrictuhex,& 252 & buhex,ceed_vector_active,err) 253 call ceedoperatorsetfield(op_diffhex,'v',erestrictuhex,& 254 & buhex,ceed_vector_active,err) 255 256! Composite Operators 257 call ceedcompositeoperatorcreate(ceed,op_setup,err) 258 call ceedcompositeoperatoraddsub(op_setup,op_setuptet,err) 259 call ceedcompositeoperatoraddsub(op_setup,op_setuphex,err) 260 261 call ceedcompositeoperatorcreate(ceed,op_diff,err) 262 call ceedcompositeoperatoraddsub(op_diff,op_difftet,err) 263 call ceedcompositeoperatoraddsub(op_diff,op_diffhex,err) 264 265! Apply Setup Operator 266 call ceedoperatorapply(op_setup,x,ceed_vector_none,& 267 & ceed_request_immediate,err) 268 269! Apply diff Operator 270 call ceedvectorcreate(ceed,ndofs,u,err) 271 call ceedvectorsetvalue(u,1.d0,err) 272 call ceedvectorcreate(ceed,ndofs,v,err) 273 274 call ceedoperatorapply(op_diff,u,v,ceed_request_immediate,err) 275 276! Check Output 277 call ceedvectorgetarrayread(v,ceed_mem_host,hv,voffset,err) 278 do i=1,ndofs 279 if (abs(hv(voffset+i))>1.0d-14) then 280! LCOV_EXCL_START 281 write(*,*) 'Computed: ',total,' != True: 0.0' 282! LCOV_EXCL_STOP 283 endif 284 enddo 285 call ceedvectorrestorearrayread(v,hv,voffset,err) 286 287! Cleanup 288 call ceedqfunctiondestroy(qf_setuptet,err) 289 call ceedqfunctiondestroy(qf_difftet,err) 290 call ceedoperatordestroy(op_setuptet,err) 291 call ceedoperatordestroy(op_difftet,err) 292 call ceedqfunctiondestroy(qf_setuphex,err) 293 call ceedqfunctiondestroy(qf_diffhex,err) 294 call ceedoperatordestroy(op_setuphex,err) 295 call ceedoperatordestroy(op_diffhex,err) 296 call ceedoperatordestroy(op_setup,err) 297 call ceedoperatordestroy(op_diff,err) 298 call ceedelemrestrictiondestroy(erestrictutet,err) 299 call ceedelemrestrictiondestroy(erestrictxtet,err) 300 call ceedelemrestrictiondestroy(erestrictqditet,err) 301 call ceedelemrestrictiondestroy(erestrictxitet,err) 302 call ceedelemrestrictiondestroy(erestrictuhex,err) 303 call ceedelemrestrictiondestroy(erestrictxhex,err) 304 call ceedelemrestrictiondestroy(erestrictqdihex,err) 305 call ceedelemrestrictiondestroy(erestrictxihex,err) 306 call ceedbasisdestroy(butet,err) 307 call ceedbasisdestroy(bxtet,err) 308 call ceedbasisdestroy(buhex,err) 309 call ceedbasisdestroy(bxhex,err) 310 call ceedvectordestroy(x,err) 311 call ceedvectordestroy(u,err) 312 call ceedvectordestroy(v,err) 313 call ceedvectordestroy(qdatatet,err) 314 call ceedvectordestroy(qdatahex,err) 315 call ceeddestroy(ceed,err) 316 end 317!----------------------------------------------------------------------- 318