1!----------------------------------------------------------------------- 2! 3! Header with common subroutine 4! 5 include 't320-basis-f.h' 6!----------------------------------------------------------------------- 7 subroutine setup_mass(ctx,q,u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,& 8& u14,u15,u16,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,& 9& v16,ierr) 10 real*8 ctx 11 real*8 u1(1) 12 real*8 u2(1) 13 real*8 v1(1) 14 integer q,ierr 15 16 do i=1,q 17 v1(i)=u2(i)*(u1(i+q*0)*u1(i+q*3)-u1(i+q*1)*u1(i+q*2)) 18 enddo 19 20 ierr=0 21 end 22!----------------------------------------------------------------------- 23 subroutine setup_diff(ctx,q,u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,& 24& u14,u15,u16,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,& 25& v16,ierr) 26 real*8 ctx 27 real*8 u1(1) 28 real*8 u2(1) 29 real*8 v1(1) 30 real*8 w 31 integer q,ierr 32 33 do i=1,q 34 w=u2(i)/(u1(i+q*0)*u1(i+q*3)-u1(i+q*1)*u1(i+q*2)) 35 v1(i+q*0)=w*(u1(i+q*2)*u1(i+q*2)+u1(i+q*3)*u1(i+q*3)) 36 v1(i+q*1)=w*(u1(i+q*0)*u1(i+q*0)+u1(i+q*1)*u1(i+q*1)) 37 v1(i+q*2)=-w*(u1(i+q*0)*u1(i+q*2)+u1(i+q*2)*u1(i+q*3)) 38 enddo 39 40 ierr=0 41 end 42!----------------------------------------------------------------------- 43 subroutine apply(ctx,q,u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,u14,& 44& u15,u16,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,v16,ierr) 45 real*8 ctx 46 real*8 u1(1) 47 real*8 u2(1) 48 real*8 u3(1) 49 real*8 u4(1) 50 real*8 v1(1) 51 real*8 v2(1) 52 real*8 du0,du1 53 integer q,ierr 54 55 do i=1,q 56! mass 57 v1(i) = u2(i)*u4(i) 58! diff 59 du0=u1(i+q*0) 60 du1=u1(i+q*1) 61 v2(i+q*0)=u3(i+q*0)*du0+u3(i+q*2)*du1 62 v2(i+q*1)=u3(i+q*2)*du0+u3(i+q*1)*du1 63 enddo 64 65 ierr=0 66 end 67!----------------------------------------------------------------------- 68 program test 69 70 include 'ceedf.h' 71 72 integer ceed,err,i 73 integer imode 74 parameter(imode=ceed_noninterlaced) 75 integer stridesx(3),stridesu(3),stridesqd(3) 76 integer erestrictx,erestrictu,erestrictxi,erestrictui,erestrictqi 77 integer bx,bu 78 integer qf_setup_mass,qf_setup_diff,qf_apply 79 integer op_setup_mass,op_setup_diff,op_apply 80 integer qdata_mass,qdata_diff,x,a,u,v 81 integer nelem,p,q,d 82 integer row,col,offset 83 parameter(nelem=12) 84 parameter(p=6) 85 parameter(q=4) 86 parameter(d=2) 87 integer ndofs,nqpts,nx,ny 88 parameter(nx=3) 89 parameter(ny=2) 90 parameter(ndofs=(nx*2+1)*(ny*2+1)) 91 parameter(nqpts=nelem*q*q) 92 integer indx(nelem*p*p) 93 real*8 arrx(d*ndofs),aa(nqpts),uu(ndofs),vv(ndofs),atrue(ndofs) 94 integer*8 xoffset,aoffset,uoffset,voffset 95 96 real*8 qref(d*q) 97 real*8 qweight(q) 98 real*8 interp(p*q) 99 real*8 grad(d*p*q) 100 101 character arg*32 102 103 external setup_mass,setup_diff,apply 104 105 call getarg(1,arg) 106 107 call ceedinit(trim(arg)//char(0),ceed,err) 108 109! DoF Coordinates 110 do i=0,ndofs-1 111 arrx(i+1)=mod(i,(nx*2+1)) 112 arrx(i+1)=arrx(i+1)*(1.d0/(nx*2.d0)) 113 val=(i/(nx*2+1)) 114 arrx(i+1+ndofs)=val*(1.d0/(ny*2.d0)) 115 enddo 116 call ceedvectorcreate(ceed,d*ndofs,x,err) 117 xoffset=0 118 call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,arrx,xoffset,err) 119 120! Qdata Vector 121 call ceedvectorcreate(ceed,nqpts,qdata_mass,err) 122 call ceedvectorcreate(ceed,nqpts*d*(d+1)/2,qdata_diff,err) 123 124! Element Setup 125 do i=0,5 126 col=mod(i,nx) 127 row=i/nx 128 offset=col*2+row*(nx*2+1)*2 129 130 indx(i*2*p+1)=2+offset 131 indx(i*2*p+2)=9+offset 132 indx(i*2*p+3)=16+offset 133 indx(i*2*p+4)=1+offset 134 indx(i*2*p+5)=8+offset 135 indx(i*2*p+6)=0+offset 136 137 indx(i*2*p+7)=14+offset 138 indx(i*2*p+8)=7+offset 139 indx(i*2*p+9)=0+offset 140 indx(i*2*p+10)=15+offset 141 indx(i*2*p+11)=8+offset 142 indx(i*2*p+12)=16+offset 143 enddo 144 145! Restrictions 146 call ceedelemrestrictioncreate(ceed,imode,nelem,p,ndofs,d,& 147 & ceed_mem_host,ceed_use_pointer,indx,erestrictx,err) 148 stridesx=[1,p,p*d] 149 call ceedelemrestrictioncreatestrided(ceed,nelem,p,& 150 & nelem*p,d,stridesx,erestrictxi,err) 151 152 call ceedelemrestrictioncreate(ceed,imode,nelem,p,ndofs,1,& 153 & ceed_mem_host,ceed_use_pointer,indx,erestrictu,err) 154 stridesu=[1,q,q] 155 call ceedelemrestrictioncreatestrided(ceed,nelem,q,nqpts,& 156 & 1,stridesu,erestrictui,err) 157 158 stridesqd=[1,q,q*d*(d+1)/2] 159 call ceedelemrestrictioncreatestrided(ceed,nelem,q,nqpts,& 160 & d*(d+1)/2,stridesqd,erestrictqi,err) 161 162! Bases 163 call buildmats(qref,qweight,interp,grad) 164 call ceedbasiscreateh1(ceed,ceed_triangle,d,p,q,interp,grad,qref,qweight,& 165 & bx,err) 166 call buildmats(qref,qweight,interp,grad) 167 call ceedbasiscreateh1(ceed,ceed_triangle,1,p,q,interp,grad,qref,qweight,& 168 & bu,err) 169 170! QFunction - setup mass 171 call ceedqfunctioncreateinterior(ceed,1,setup_mass,& 172 &SOURCE_DIR& 173 &//'t532-operator.h:setup_mass'//char(0),qf_setup_mass,err) 174 call ceedqfunctionaddinput(qf_setup_mass,'dx',d*d,ceed_eval_grad,err) 175 call ceedqfunctionaddinput(qf_setup_mass,'_weight',1,ceed_eval_weight,err) 176 call ceedqfunctionaddoutput(qf_setup_mass,'qdata',1,ceed_eval_none,err) 177 178! Operator - setup mass 179 call ceedoperatorcreate(ceed,qf_setup_mass,ceed_qfunction_none,& 180 & ceed_qfunction_none,op_setup_mass,err) 181 call ceedoperatorsetfield(op_setup_mass,'dx',erestrictx,& 182 & bx,ceed_vector_active,err) 183 call ceedoperatorsetfield(op_setup_mass,'_weight',erestrictxi,& 184 & bx,ceed_vector_none,err) 185 call ceedoperatorsetfield(op_setup_mass,'qdata',erestrictui,& 186 & ceed_basis_collocated,ceed_vector_active,err) 187 188! QFunction - setup diff 189 call ceedqfunctioncreateinterior(ceed,1,setup_diff,& 190 &SOURCE_DIR& 191 &//'t532-operator.h:setup_diff'//char(0),qf_setup_diff,err) 192 call ceedqfunctionaddinput(qf_setup_diff,'dx',d*d,ceed_eval_grad,err) 193 call ceedqfunctionaddinput(qf_setup_diff,'_weight',1,ceed_eval_weight,err) 194 call ceedqfunctionaddoutput(qf_setup_diff,'qdata',& 195 & d*(d+1)/2,ceed_eval_none,err) 196 197! Operator - setup diff 198 call ceedoperatorcreate(ceed,qf_setup_diff,ceed_qfunction_none,& 199 & ceed_qfunction_none,op_setup_diff,err) 200 call ceedoperatorsetfield(op_setup_diff,'dx',erestrictx,& 201 & bx,ceed_vector_active,err) 202 call ceedoperatorsetfield(op_setup_diff,'_weight',erestrictxi,& 203 & bx,ceed_vector_none,err) 204 call ceedoperatorsetfield(op_setup_diff,'qdata',erestrictqi,& 205 & ceed_basis_collocated,ceed_vector_active,err) 206 207! Apply Setup Operators 208 call ceedoperatorapply(op_setup_mass,x,qdata_mass,& 209 & ceed_request_immediate,err) 210 call ceedoperatorapply(op_setup_diff,x,qdata_diff,& 211 & ceed_request_immediate,err) 212 213! QFunction - apply 214 call ceedqfunctioncreateinterior(ceed,1,apply,& 215 &SOURCE_DIR& 216 &//'t532-operator.h:apply'//char(0),qf_apply,err) 217 call ceedqfunctionaddinput(qf_apply,'du',d,ceed_eval_grad,err) 218 call ceedqfunctionaddinput(qf_apply,'qdata_mass',1,ceed_eval_none,err) 219 call ceedqfunctionaddinput(qf_apply,'qdata_diff',& 220 & d*(d+1)/2,ceed_eval_none,err) 221 call ceedqfunctionaddinput(qf_apply,'u',1,ceed_eval_interp,err) 222 call ceedqfunctionaddoutput(qf_apply,'v',1,ceed_eval_interp,err) 223 call ceedqfunctionaddoutput(qf_apply,'dv',d,ceed_eval_grad,err) 224 225! Operator - apply 226 call ceedoperatorcreate(ceed,qf_apply,ceed_qfunction_none,& 227 & ceed_qfunction_none,op_apply,err) 228 call ceedoperatorsetfield(op_apply,'du',erestrictu,& 229 & bu,ceed_vector_active,err) 230 call ceedoperatorsetfield(op_apply,'qdata_mass',erestrictui,& 231 & ceed_basis_collocated,qdata_mass,err) 232 call ceedoperatorsetfield(op_apply,'qdata_diff',erestrictqi,& 233 & ceed_basis_collocated,qdata_diff,err) 234 call ceedoperatorsetfield(op_apply,'u',erestrictu,& 235 & bu,ceed_vector_active,err) 236 call ceedoperatorsetfield(op_apply,'v',erestrictu,& 237 & bu,ceed_vector_active,err) 238 call ceedoperatorsetfield(op_apply,'dv',erestrictu,& 239 & bu,ceed_vector_active,err) 240 241! Assemble Diagonal 242 call ceedoperatorassemblelineardiagonal(op_apply,a,& 243 & ceed_request_immediate,err) 244 245! Manually assemble diagonal 246 call ceedvectorcreate(ceed,ndofs,u,err) 247 call ceedvectorsetvalue(u,0.d0,err) 248 call ceedvectorcreate(ceed,ndofs,v,err) 249 do i=1,ndofs 250 call ceedvectorgetarray(u,ceed_mem_host,uu,uoffset,err) 251 uu(i+uoffset)=1.d0 252 if (i>1) then 253 uu(i-1+uoffset)=0.d0 254 endif 255 call ceedvectorrestorearray(u,uu,uoffset,err) 256 257 call ceedoperatorapply(op_apply,u,v,ceed_request_immediate,err) 258 259 call ceedvectorgetarrayread(v,ceed_mem_host,vv,voffset,err) 260 atrue(i)=vv(voffset+i) 261 call ceedvectorrestorearrayread(v,vv,voffset,err) 262 enddo 263 264! Check Output 265 call ceedvectorgetarrayread(a,ceed_mem_host,aa,aoffset,err) 266 do i=1,ndofs 267 if (abs(aa(aoffset+i)-atrue(i))>1.0d-14) then 268! LCOV_EXCL_START 269 write(*,*) '[',i,'] Error in assembly: ',aa(aoffset+i),' != ',& 270 & atrue(i) 271! LCOV_EXCL_STOP 272 endif 273 enddo 274 call ceedvectorrestorearrayread(a,aa,aoffset,err) 275 276! Cleanup 277 call ceedqfunctiondestroy(qf_setup_mass,err) 278 call ceedqfunctiondestroy(qf_setup_diff,err) 279 call ceedqfunctiondestroy(qf_apply,err) 280 call ceedoperatordestroy(op_setup_mass,err) 281 call ceedoperatordestroy(op_setup_diff,err) 282 call ceedoperatordestroy(op_apply,err) 283 call ceedelemrestrictiondestroy(erestrictu,err) 284 call ceedelemrestrictiondestroy(erestrictx,err) 285 call ceedelemrestrictiondestroy(erestrictxi,err) 286 call ceedelemrestrictiondestroy(erestrictui,err) 287 call ceedelemrestrictiondestroy(erestrictqi,err) 288 call ceedbasisdestroy(bu,err) 289 call ceedbasisdestroy(bx,err) 290 call ceedvectordestroy(x,err) 291 call ceedvectordestroy(a,err) 292 call ceedvectordestroy(u,err) 293 call ceedvectordestroy(v,err) 294 call ceedvectordestroy(qdata_mass,err) 295 call ceedvectordestroy(qdata_diff,err) 296 call ceeddestroy(ceed,err) 297 end 298!----------------------------------------------------------------------- 299