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