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 ceedvectorcreate(ceed,ndofs,a,err) 238 call ceedoperatorlinearassemblediagonal(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(erestrictui,err) 282 call ceedelemrestrictiondestroy(erestrictqi,err) 283 call ceedbasisdestroy(bu,err) 284 call ceedbasisdestroy(bx,err) 285 call ceedvectordestroy(x,err) 286 call ceedvectordestroy(a,err) 287 call ceedvectordestroy(u,err) 288 call ceedvectordestroy(v,err) 289 call ceedvectordestroy(qdata_mass,err) 290 call ceedvectordestroy(qdata_diff,err) 291 call ceeddestroy(ceed,err) 292 end 293!----------------------------------------------------------------------- 294