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