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 stridesu(3),stridesqd(3) 76 integer erestrictx,erestrictu,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 149 call ceedelemrestrictioncreate(ceed,imode,nelem,p,ndofs,1,& 150 & ceed_mem_host,ceed_use_pointer,indx,erestrictu,err) 151 stridesu=[1,q,q] 152 call ceedelemrestrictioncreatestrided(ceed,nelem,q,nqpts,& 153 & 1,stridesu,erestrictui,err) 154 155 stridesqd=[1,q,q*d*(d+1)/2] 156 call ceedelemrestrictioncreatestrided(ceed,nelem,q,nqpts,& 157 & d*(d+1)/2,stridesqd,erestrictqi,err) 158 159! Bases 160 call buildmats(qref,qweight,interp,grad) 161 call ceedbasiscreateh1(ceed,ceed_triangle,d,p,q,interp,grad,qref,qweight,& 162 & bx,err) 163 call buildmats(qref,qweight,interp,grad) 164 call ceedbasiscreateh1(ceed,ceed_triangle,1,p,q,interp,grad,qref,qweight,& 165 & bu,err) 166 167! QFunction - setup mass 168 call ceedqfunctioncreateinterior(ceed,1,setup_mass,& 169 &SOURCE_DIR& 170 &//'t532-operator.h:setup_mass'//char(0),qf_setup_mass,err) 171 call ceedqfunctionaddinput(qf_setup_mass,'dx',d*d,ceed_eval_grad,err) 172 call ceedqfunctionaddinput(qf_setup_mass,'_weight',1,ceed_eval_weight,err) 173 call ceedqfunctionaddoutput(qf_setup_mass,'qdata',1,ceed_eval_none,err) 174 175! Operator - setup mass 176 call ceedoperatorcreate(ceed,qf_setup_mass,ceed_qfunction_none,& 177 & ceed_qfunction_none,op_setup_mass,err) 178 call ceedoperatorsetfield(op_setup_mass,'dx',erestrictx,& 179 & bx,ceed_vector_active,err) 180 call ceedoperatorsetfield(op_setup_mass,'_weight',& 181 & ceed_elemrestriction_none,bx,ceed_vector_none,err) 182 call ceedoperatorsetfield(op_setup_mass,'qdata',erestrictui,& 183 & ceed_basis_collocated,ceed_vector_active,err) 184 185! QFunction - setup diff 186 call ceedqfunctioncreateinterior(ceed,1,setup_diff,& 187 &SOURCE_DIR& 188 &//'t532-operator.h:setup_diff'//char(0),qf_setup_diff,err) 189 call ceedqfunctionaddinput(qf_setup_diff,'dx',d*d,ceed_eval_grad,err) 190 call ceedqfunctionaddinput(qf_setup_diff,'_weight',1,ceed_eval_weight,err) 191 call ceedqfunctionaddoutput(qf_setup_diff,'qdata',& 192 & d*(d+1)/2,ceed_eval_none,err) 193 194! Operator - setup diff 195 call ceedoperatorcreate(ceed,qf_setup_diff,ceed_qfunction_none,& 196 & ceed_qfunction_none,op_setup_diff,err) 197 call ceedoperatorsetfield(op_setup_diff,'dx',erestrictx,& 198 & bx,ceed_vector_active,err) 199 call ceedoperatorsetfield(op_setup_diff,'_weight',& 200 & ceed_elemrestriction_none,bx,ceed_vector_none,err) 201 call ceedoperatorsetfield(op_setup_diff,'qdata',erestrictqi,& 202 & ceed_basis_collocated,ceed_vector_active,err) 203 204! Apply Setup Operators 205 call ceedoperatorapply(op_setup_mass,x,qdata_mass,& 206 & ceed_request_immediate,err) 207 call ceedoperatorapply(op_setup_diff,x,qdata_diff,& 208 & ceed_request_immediate,err) 209 210! QFunction - apply 211 call ceedqfunctioncreateinterior(ceed,1,apply,& 212 &SOURCE_DIR& 213 &//'t532-operator.h:apply'//char(0),qf_apply,err) 214 call ceedqfunctionaddinput(qf_apply,'du',d,ceed_eval_grad,err) 215 call ceedqfunctionaddinput(qf_apply,'qdata_mass',1,ceed_eval_none,err) 216 call ceedqfunctionaddinput(qf_apply,'qdata_diff',& 217 & d*(d+1)/2,ceed_eval_none,err) 218 call ceedqfunctionaddinput(qf_apply,'u',1,ceed_eval_interp,err) 219 call ceedqfunctionaddoutput(qf_apply,'v',1,ceed_eval_interp,err) 220 call ceedqfunctionaddoutput(qf_apply,'dv',d,ceed_eval_grad,err) 221 222! Operator - apply 223 call ceedoperatorcreate(ceed,qf_apply,ceed_qfunction_none,& 224 & ceed_qfunction_none,op_apply,err) 225 call ceedoperatorsetfield(op_apply,'du',erestrictu,& 226 & bu,ceed_vector_active,err) 227 call ceedoperatorsetfield(op_apply,'qdata_mass',erestrictui,& 228 & ceed_basis_collocated,qdata_mass,err) 229 call ceedoperatorsetfield(op_apply,'qdata_diff',erestrictqi,& 230 & ceed_basis_collocated,qdata_diff,err) 231 call ceedoperatorsetfield(op_apply,'u',erestrictu,& 232 & bu,ceed_vector_active,err) 233 call ceedoperatorsetfield(op_apply,'v',erestrictu,& 234 & bu,ceed_vector_active,err) 235 call ceedoperatorsetfield(op_apply,'dv',erestrictu,& 236 & bu,ceed_vector_active,err) 237 238! Assemble Diagonal 239 call ceedoperatorassemblelineardiagonal(op_apply,a,& 240 & ceed_request_immediate,err) 241 242! Manually assemble diagonal 243 call ceedvectorcreate(ceed,ndofs,u,err) 244 call ceedvectorsetvalue(u,0.d0,err) 245 call ceedvectorcreate(ceed,ndofs,v,err) 246 do i=1,ndofs 247 call ceedvectorgetarray(u,ceed_mem_host,uu,uoffset,err) 248 uu(i+uoffset)=1.d0 249 if (i>1) then 250 uu(i-1+uoffset)=0.d0 251 endif 252 call ceedvectorrestorearray(u,uu,uoffset,err) 253 254 call ceedoperatorapply(op_apply,u,v,ceed_request_immediate,err) 255 256 call ceedvectorgetarrayread(v,ceed_mem_host,vv,voffset,err) 257 atrue(i)=vv(voffset+i) 258 call ceedvectorrestorearrayread(v,vv,voffset,err) 259 enddo 260 261! Check Output 262 call ceedvectorgetarrayread(a,ceed_mem_host,aa,aoffset,err) 263 do i=1,ndofs 264 if (abs(aa(aoffset+i)-atrue(i))>1.0d-14) then 265! LCOV_EXCL_START 266 write(*,*) '[',i,'] Error in assembly: ',aa(aoffset+i),' != ',& 267 & atrue(i) 268! LCOV_EXCL_STOP 269 endif 270 enddo 271 call ceedvectorrestorearrayread(a,aa,aoffset,err) 272 273! Cleanup 274 call ceedqfunctiondestroy(qf_setup_mass,err) 275 call ceedqfunctiondestroy(qf_setup_diff,err) 276 call ceedqfunctiondestroy(qf_apply,err) 277 call ceedoperatordestroy(op_setup_mass,err) 278 call ceedoperatordestroy(op_setup_diff,err) 279 call ceedoperatordestroy(op_apply,err) 280 call ceedelemrestrictiondestroy(erestrictu,err) 281 call ceedelemrestrictiondestroy(erestrictx,err) 282 call ceedelemrestrictiondestroy(erestrictui,err) 283 call ceedelemrestrictiondestroy(erestrictqi,err) 284 call ceedbasisdestroy(bu,err) 285 call ceedbasisdestroy(bx,err) 286 call ceedvectordestroy(x,err) 287 call ceedvectordestroy(a,err) 288 call ceedvectordestroy(u,err) 289 call ceedvectordestroy(v,err) 290 call ceedvectordestroy(qdata_mass,err) 291 call ceedvectordestroy(qdata_diff,err) 292 call ceeddestroy(ceed,err) 293 end 294!----------------------------------------------------------------------- 295