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 stridesqdtet(3),stridesqdhex(3) 50 integer erestrictxtet,erestrictutet,erestrictqditet,& 51& erestrictxhex,erestrictuhex,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 136 call ceedelemrestrictioncreate(ceed,imode,nelemtet,ptet,ndofs,1,& 137 & ceed_mem_host,ceed_use_pointer,indxtet,erestrictutet,err) 138 stridesqdtet=[1,qtet,qtet*d*(d+1)/2] 139 call ceedelemrestrictioncreatestrided(ceed,nelemtet,qtet,nqptstet,& 140 & d*(d+1)/2,stridesqdtet,erestrictqditet,err) 141 142! -- Bases 143 call buildmats(qref,qweight,interp,grad) 144 call ceedbasiscreateh1(ceed,ceed_triangle,d,ptet,qtet,interp,grad,qref,& 145 & qweight,bxtet,err) 146 call buildmats(qref,qweight,interp,grad) 147 call ceedbasiscreateh1(ceed,ceed_triangle,1,ptet,qtet,interp,grad,qref,& 148 & qweight,butet,err) 149 150! -- QFunctions 151 call ceedqfunctioncreateinterior(ceed,1,setup,& 152 &SOURCE_DIR& 153 &//'t522-operator.h:setup'//char(0),qf_setuptet,err) 154 call ceedqfunctionaddinput(qf_setuptet,'_weight',1,ceed_eval_weight,err) 155 call ceedqfunctionaddinput(qf_setuptet,'dx',d*d,ceed_eval_grad,err) 156 call ceedqfunctionaddoutput(qf_setuptet,'rho',d*(d+1)/2,ceed_eval_none,& 157 & err) 158 159 call ceedqfunctioncreateinterior(ceed,1,diff,& 160 &SOURCE_DIR& 161 &//'t522-operator.h:diff'//char(0),qf_difftet,err) 162 call ceedqfunctionaddinput(qf_difftet,'rho',d*(d+1)/2,ceed_eval_none,err) 163 call ceedqfunctionaddinput(qf_difftet,'u',d,ceed_eval_grad,err) 164 call ceedqfunctionaddoutput(qf_difftet,'v',d,ceed_eval_grad,err) 165 166! -- Operators 167! ---- Setup Tet 168 call ceedoperatorcreate(ceed,qf_setuptet,ceed_qfunction_none,& 169 & ceed_qfunction_none,op_setuptet,err) 170 call ceedoperatorsetfield(op_setuptet,'_weight',& 171 & ceed_elemrestriction_none,bxtet,ceed_vector_none,err) 172 call ceedoperatorsetfield(op_setuptet,'dx',erestrictxtet,& 173 & bxtet,ceed_vector_active,err) 174 call ceedoperatorsetfield(op_setuptet,'rho',erestrictqditet,& 175 & ceed_basis_collocated,qdatatet,err) 176! ---- diff Tet 177 call ceedoperatorcreate(ceed,qf_difftet,ceed_qfunction_none,& 178 & ceed_qfunction_none,op_difftet,err) 179 call ceedoperatorsetfield(op_difftet,'rho',erestrictqditet,& 180 & ceed_basis_collocated,qdatatet,err) 181 call ceedoperatorsetfield(op_difftet,'u',erestrictutet,& 182 & butet,ceed_vector_active,err) 183 call ceedoperatorsetfield(op_difftet,'v',erestrictutet,& 184 & butet,ceed_vector_active,err) 185 186! Hex Elements 187 do i=0,nelemhex-1 188 col=mod(i,nx) 189 row=i/nx 190 offset=(nxtet*2+1)*(nytet*2)*(1+row)+col*2 191 do j=0,phex-1 192 do k=0,phex-1 193 indxhex(phex*(phex*i+k)+j+1)=offset+k*(nxhex*2+1)+j 194 enddo 195 enddo 196 enddo 197 198! -- Restrictions 199 call ceedelemrestrictioncreate(ceed,imode,nelemhex,phex*phex,ndofs,d,& 200 & ceed_mem_host,ceed_use_pointer,indxhex,erestrictxhex,err) 201 202 call ceedelemrestrictioncreate(ceed,imode,nelemhex,phex*phex,ndofs,1,& 203 & ceed_mem_host,ceed_use_pointer,indxhex,erestrictuhex,err) 204 stridesqdhex=[1,qhex*qhex,qhex*qhex*d*(d+1)/2] 205 call ceedelemrestrictioncreatestrided(ceed,nelemhex,qhex*qhex,& 206 & nqptshex,d*(d+1)/2,stridesqdhex,erestrictqdihex,err) 207 208! -- Bases 209 call ceedbasiscreatetensorh1lagrange(ceed,d,d,phex,qhex,ceed_gauss,& 210 & bxhex,err) 211 call ceedbasiscreatetensorh1lagrange(ceed,d,1,phex,qhex,ceed_gauss,& 212 & buhex,err) 213 214! -- QFunctions 215 call ceedqfunctioncreateinterior(ceed,1,setup,& 216 &SOURCE_DIR& 217 &//'t522-operator.h:setup'//char(0),qf_setuphex,err) 218 call ceedqfunctionaddinput(qf_setuphex,'_weight',1,ceed_eval_weight,err) 219 call ceedqfunctionaddinput(qf_setuphex,'dx',d*d,ceed_eval_grad,err) 220 call ceedqfunctionaddoutput(qf_setuphex,'rho',d*(d+1)/2,ceed_eval_none,& 221 & err) 222 223 call ceedqfunctioncreateinterior(ceed,1,diff,& 224 &SOURCE_DIR& 225 &//'t522-operator.h:diff'//char(0),qf_diffhex,err) 226 call ceedqfunctionaddinput(qf_diffhex,'rho',d*(d+1)/2,ceed_eval_none,err) 227 call ceedqfunctionaddinput(qf_diffhex,'u',d,ceed_eval_grad,err) 228 call ceedqfunctionaddoutput(qf_diffhex,'v',d,ceed_eval_grad,err) 229 230! -- Operators 231! ---- Setup Hex 232 call ceedoperatorcreate(ceed,qf_setuphex,ceed_qfunction_none,& 233 & ceed_qfunction_none,op_setuphex,err) 234 call ceedoperatorsetfield(op_setuphex,'_weight',& 235 & ceed_elemrestriction_none,bxhex,ceed_vector_none,err) 236 call ceedoperatorsetfield(op_setuphex,'dx',erestrictxhex,& 237 & bxhex,ceed_vector_active,err) 238 call ceedoperatorsetfield(op_setuphex,'rho',erestrictqdihex,& 239 & ceed_basis_collocated,qdatahex,err) 240! ---- diff Hex 241 call ceedoperatorcreate(ceed,qf_diffhex,ceed_qfunction_none,& 242 & ceed_qfunction_none,op_diffhex,err) 243 call ceedoperatorsetfield(op_diffhex,'rho',erestrictqdihex,& 244 & ceed_basis_collocated,qdatahex,err) 245 call ceedoperatorsetfield(op_diffhex,'u',erestrictuhex,& 246 & buhex,ceed_vector_active,err) 247 call ceedoperatorsetfield(op_diffhex,'v',erestrictuhex,& 248 & buhex,ceed_vector_active,err) 249 250! Composite Operators 251 call ceedcompositeoperatorcreate(ceed,op_setup,err) 252 call ceedcompositeoperatoraddsub(op_setup,op_setuptet,err) 253 call ceedcompositeoperatoraddsub(op_setup,op_setuphex,err) 254 255 call ceedcompositeoperatorcreate(ceed,op_diff,err) 256 call ceedcompositeoperatoraddsub(op_diff,op_difftet,err) 257 call ceedcompositeoperatoraddsub(op_diff,op_diffhex,err) 258 259! Apply Setup Operator 260 call ceedoperatorapply(op_setup,x,ceed_vector_none,& 261 & ceed_request_immediate,err) 262 263! Apply diff Operator 264 call ceedvectorcreate(ceed,ndofs,u,err) 265 call ceedvectorsetvalue(u,1.d0,err) 266 call ceedvectorcreate(ceed,ndofs,v,err) 267 268 call ceedoperatorapply(op_diff,u,v,ceed_request_immediate,err) 269 270! Check Output 271 call ceedvectorgetarrayread(v,ceed_mem_host,hv,voffset,err) 272 do i=1,ndofs 273 if (abs(hv(voffset+i))>1.0d-14) then 274! LCOV_EXCL_START 275 write(*,*) 'Computed: ',total,' != True: 0.0' 276! LCOV_EXCL_STOP 277 endif 278 enddo 279 call ceedvectorrestorearrayread(v,hv,voffset,err) 280 281! Cleanup 282 call ceedqfunctiondestroy(qf_setuptet,err) 283 call ceedqfunctiondestroy(qf_difftet,err) 284 call ceedoperatordestroy(op_setuptet,err) 285 call ceedoperatordestroy(op_difftet,err) 286 call ceedqfunctiondestroy(qf_setuphex,err) 287 call ceedqfunctiondestroy(qf_diffhex,err) 288 call ceedoperatordestroy(op_setuphex,err) 289 call ceedoperatordestroy(op_diffhex,err) 290 call ceedoperatordestroy(op_setup,err) 291 call ceedoperatordestroy(op_diff,err) 292 call ceedelemrestrictiondestroy(erestrictutet,err) 293 call ceedelemrestrictiondestroy(erestrictxtet,err) 294 call ceedelemrestrictiondestroy(erestrictqditet,err) 295 call ceedelemrestrictiondestroy(erestrictuhex,err) 296 call ceedelemrestrictiondestroy(erestrictxhex,err) 297 call ceedelemrestrictiondestroy(erestrictqdihex,err) 298 call ceedbasisdestroy(butet,err) 299 call ceedbasisdestroy(bxtet,err) 300 call ceedbasisdestroy(buhex,err) 301 call ceedbasisdestroy(bxhex,err) 302 call ceedvectordestroy(x,err) 303 call ceedvectordestroy(u,err) 304 call ceedvectordestroy(v,err) 305 call ceedvectordestroy(qdatatet,err) 306 call ceedvectordestroy(qdatahex,err) 307 call ceeddestroy(ceed,err) 308 end 309!----------------------------------------------------------------------- 310