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 erestrictx,erestrictu,erestrictxi,erestrictui,erestrictqi 74 integer bx,bu 75 integer qf_setup_mass,qf_setup_diff,qf_apply 76 integer op_setup_mass,op_setup_diff,op_apply 77 integer qdata_mass,qdata_diff,x,a,u,v 78 integer nelem,p,q,d 79 integer row,col,offset 80 parameter(nelem=12) 81 parameter(p=6) 82 parameter(q=4) 83 parameter(d=2) 84 integer ndofs,nqpts,nx,ny 85 parameter(nx=3) 86 parameter(ny=2) 87 parameter(ndofs=(nx*2+1)*(ny*2+1)) 88 parameter(nqpts=nelem*q*q) 89 integer indx(nelem*p*p) 90 real*8 arrx(d*ndofs),aa(nqpts),uu(ndofs),vv(ndofs),atrue(ndofs) 91 integer*8 xoffset,aoffset,uoffset,voffset 92 93 real*8 qref(d*q) 94 real*8 qweight(q) 95 real*8 interp(p*q) 96 real*8 grad(d*p*q) 97 98 character arg*32 99 100 external setup_mass,setup_diff,apply,apply_lin 101 102 call getarg(1,arg) 103 104 call ceedinit(trim(arg)//char(0),ceed,err) 105 106! DoF Coordinates 107 do i=0,ndofs-1 108 arrx(i+1)=mod(i,(nx*2+1)) 109 arrx(i+1)=arrx(i+1)*(1.d0/(nx*2.d0)) 110 val=(i/(nx*2+1)) 111 arrx(i+1+ndofs)=val*(1.d0/(ny*2.d0)) 112 enddo 113 call ceedvectorcreate(ceed,d*ndofs,x,err) 114 xoffset=0 115 call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,arrx,xoffset,err) 116 117! Qdata Vector 118 call ceedvectorcreate(ceed,nqpts,qdata_mass,err) 119 call ceedvectorcreate(ceed,nqpts*d*(d+1)/2,qdata_diff,err) 120 121! Element Setup 122 do i=0,5 123 col=mod(i,nx) 124 row=i/nx 125 offset=col*2+row*(nx*2+1)*2 126 127 indx(i*2*p+1)=2+offset 128 indx(i*2*p+2)=9+offset 129 indx(i*2*p+3)=16+offset 130 indx(i*2*p+4)=1+offset 131 indx(i*2*p+5)=8+offset 132 indx(i*2*p+6)=0+offset 133 134 indx(i*2*p+7)=14+offset 135 indx(i*2*p+8)=7+offset 136 indx(i*2*p+9)=0+offset 137 indx(i*2*p+10)=15+offset 138 indx(i*2*p+11)=8+offset 139 indx(i*2*p+12)=16+offset 140 enddo 141 142! Restrictions 143 call ceedelemrestrictioncreate(ceed,nelem,p,ndofs,d,ceed_mem_host,& 144 & ceed_use_pointer,indx,erestrictx,err) 145 call ceedelemrestrictioncreateidentity(ceed,nelem,p,nelem*p,d,& 146 & erestrictxi,err) 147 148 call ceedelemrestrictioncreate(ceed,nelem,p,ndofs,1,ceed_mem_host,& 149 & ceed_use_pointer,indx,erestrictu,err) 150 call ceedelemrestrictioncreateidentity(ceed,nelem,q,nqpts,1,& 151 & erestrictui,err) 152 153 call ceedelemrestrictioncreateidentity(ceed,nelem,q,nqpts,d*(d+1)/2,& 154 & erestrictqi,err) 155 156! Bases 157 call buildmats(qref,qweight,interp,grad) 158 call ceedbasiscreateh1(ceed,ceed_triangle,d,p,q,interp,grad,qref,qweight,& 159 & bx,err) 160 call buildmats(qref,qweight,interp,grad) 161 call ceedbasiscreateh1(ceed,ceed_triangle,1,p,q,interp,grad,qref,qweight,& 162 & bu,err) 163 164! QFunction - setup mass 165 call ceedqfunctioncreateinterior(ceed,1,setup_mass,& 166 &SOURCE_DIR& 167 &//'t532-operator.h:setup_mass'//char(0),qf_setup_mass,err) 168 call ceedqfunctionaddinput(qf_setup_mass,'dx',d*d,ceed_eval_grad,err) 169 call ceedqfunctionaddinput(qf_setup_mass,'_weight',1,ceed_eval_weight,err) 170 call ceedqfunctionaddoutput(qf_setup_mass,'qdata',1,ceed_eval_none,err) 171 172! Operator - setup mass 173 call ceedoperatorcreate(ceed,qf_setup_mass,ceed_null,ceed_null,& 174 & op_setup_mass,err) 175 call ceedoperatorsetfield(op_setup_mass,'dx',erestrictx,& 176 & ceed_notranspose,bx,ceed_vector_active,err) 177 call ceedoperatorsetfield(op_setup_mass,'_weight',erestrictxi,& 178 & ceed_notranspose,bx,ceed_vector_none,err) 179 call ceedoperatorsetfield(op_setup_mass,'qdata',erestrictui,& 180 & ceed_notranspose,ceed_basis_collocated,ceed_vector_active,err) 181 182! QFunction - setup diff 183 call ceedqfunctioncreateinterior(ceed,1,setup_diff,& 184 &SOURCE_DIR& 185 &//'t532-operator.h:setup_diff'//char(0),qf_setup_diff,err) 186 call ceedqfunctionaddinput(qf_setup_diff,'dx',d*d,ceed_eval_grad,err) 187 call ceedqfunctionaddinput(qf_setup_diff,'_weight',1,ceed_eval_weight,err) 188 call ceedqfunctionaddoutput(qf_setup_diff,'qdata',& 189 & d*(d+1)/2,ceed_eval_none,err) 190 191! Operator - setup diff 192 call ceedoperatorcreate(ceed,qf_setup_diff,ceed_null,ceed_null,& 193 & op_setup_diff,err) 194 call ceedoperatorsetfield(op_setup_diff,'dx',erestrictx,& 195 & ceed_notranspose,bx,ceed_vector_active,err) 196 call ceedoperatorsetfield(op_setup_diff,'_weight',erestrictxi,& 197 & ceed_notranspose,bx,ceed_vector_none,err) 198 call ceedoperatorsetfield(op_setup_diff,'qdata',erestrictqi,& 199 & ceed_notranspose,ceed_basis_collocated,ceed_vector_active,err) 200 201! Apply Setup Operators 202 call ceedoperatorapply(op_setup_mass,x,qdata_mass,& 203 & ceed_request_immediate,err) 204 call ceedoperatorapply(op_setup_diff,x,qdata_diff,& 205 & ceed_request_immediate,err) 206 207! QFunction - apply 208 call ceedqfunctioncreateinterior(ceed,1,apply,& 209 &SOURCE_DIR& 210 &//'t532-operator.h:apply'//char(0),qf_apply,err) 211 call ceedqfunctionaddinput(qf_apply,'du',d,ceed_eval_grad,err) 212 call ceedqfunctionaddinput(qf_apply,'qdata_mass',1,ceed_eval_none,err) 213 call ceedqfunctionaddinput(qf_apply,'qdata_diff',& 214 & d*(d+1)/2,ceed_eval_none,err) 215 call ceedqfunctionaddinput(qf_apply,'u',1,ceed_eval_interp,err) 216 call ceedqfunctionaddoutput(qf_apply,'v',1,ceed_eval_interp,err) 217 call ceedqfunctionaddoutput(qf_apply,'dv',d,ceed_eval_grad,err) 218 219! Operator - apply 220 call ceedoperatorcreate(ceed,qf_apply,ceed_null,ceed_null,op_apply,& 221 & err) 222 call ceedoperatorsetfield(op_apply,'du',erestrictu,& 223 & ceed_notranspose,bu,ceed_vector_active,err) 224 call ceedoperatorsetfield(op_apply,'qdata_mass',erestrictui,& 225 & ceed_notranspose,ceed_basis_collocated,qdata_mass,err) 226 call ceedoperatorsetfield(op_apply,'qdata_diff',erestrictqi,& 227 & ceed_notranspose,ceed_basis_collocated,qdata_diff,err) 228 call ceedoperatorsetfield(op_apply,'u',erestrictu,& 229 & ceed_notranspose,bu,ceed_vector_active,err) 230 call ceedoperatorsetfield(op_apply,'v',erestrictu,& 231 & ceed_notranspose,bu,ceed_vector_active,err) 232 call ceedoperatorsetfield(op_apply,'dv',erestrictu,& 233 & ceed_notranspose,bu,ceed_vector_active,err) 234 235! Assemble Diagonal 236 call ceedoperatorassemblelineardiagonal(op_apply,a,& 237 & ceed_request_immediate,err) 238 239! Manually assemble diagonal 240 call ceedvectorcreate(ceed,ndofs,u,err) 241 call ceedvectorsetvalue(u,0.d0,err) 242 call ceedvectorcreate(ceed,ndofs,v,err) 243 do i=1,ndofs 244 call ceedvectorgetarray(u,ceed_mem_host,uu,uoffset,err) 245 uu(i+uoffset)=1.d0 246 if (i>1) then 247 uu(i-1+uoffset)=0.d0 248 endif 249 call ceedvectorrestorearray(u,uu,uoffset,err) 250 251 call ceedoperatorapply(op_apply,u,v,ceed_request_immediate,err) 252 253 call ceedvectorgetarrayread(v,ceed_mem_host,vv,voffset,err) 254 atrue(i)=vv(voffset+i) 255 call ceedvectorrestorearrayread(v,vv,voffset,err) 256 enddo 257 258! Check Output 259 call ceedvectorgetarrayread(a,ceed_mem_host,aa,aoffset,err) 260 do i=1,ndofs 261 if (abs(aa(aoffset+i)-atrue(i))>1.0d-14) then 262! LCOV_EXCL_START 263 write(*,*) '[',i,'] Error in assembly: ',aa(aoffset+i),' != ',& 264 & atrue(i) 265! LCOV_EXCL_STOP 266 endif 267 enddo 268 call ceedvectorrestorearrayread(a,aa,aoffset,err) 269 270! Cleanup 271 call ceedqfunctiondestroy(qf_setup_mass,err) 272 call ceedqfunctiondestroy(qf_setup_diff,err) 273 call ceedqfunctiondestroy(qf_apply,err) 274 call ceedoperatordestroy(op_setup_mass,err) 275 call ceedoperatordestroy(op_setup_diff,err) 276 call ceedoperatordestroy(op_apply,err) 277 call ceedelemrestrictiondestroy(erestrictu,err) 278 call ceedelemrestrictiondestroy(erestrictx,err) 279 call ceedelemrestrictiondestroy(erestrictxi,err) 280 call ceedelemrestrictiondestroy(erestrictui,err) 281 call ceedelemrestrictiondestroy(erestrictqi,err) 282 call ceedbasisdestroy(bu,err) 283 call ceedbasisdestroy(bx,err) 284 call ceedvectordestroy(x,err) 285 call ceedvectordestroy(a,err) 286 call ceedvectordestroy(u,err) 287 call ceedvectordestroy(v,err) 288 call ceedvectordestroy(qdata_mass,err) 289 call ceedvectordestroy(qdata_diff,err) 290 call ceeddestroy(ceed,err) 291 end 292!----------------------------------------------------------------------- 293