1!----------------------------------------------------------------------- 2 subroutine setup_mass(ctx,q,u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,& 3& u14,u15,u16,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,& 4& v16,ierr) 5 real*8 ctx 6 real*8 u1(1) 7 real*8 u2(1) 8 real*8 v1(1) 9 integer q,ierr 10 11 do i=1,q 12 v1(i)=u2(i)*(u1(i+q*0)*u1(i+q*3)-u1(i+q*1)*u1(i+q*2)) 13 enddo 14 15 ierr=0 16 end 17!----------------------------------------------------------------------- 18 subroutine setup_diff(ctx,q,u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,& 19& u14,u15,u16,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,& 20& v16,ierr) 21 real*8 ctx 22 real*8 u1(1) 23 real*8 u2(1) 24 real*8 v1(1) 25 real*8 w 26 integer q,ierr 27 28 do i=1,q 29 w=u2(i)/(u1(i+q*0)*u1(i+q*3)-u1(i+q*1)*u1(i+q*2)) 30 v1(i+q*0)=w*(u1(i+q*2)*u1(i+q*2)+u1(i+q*3)*u1(i+q*3)) 31 v1(i+q*1)=w*(u1(i+q*0)*u1(i+q*0)+u1(i+q*1)*u1(i+q*1)) 32 v1(i+q*2)=-w*(u1(i+q*0)*u1(i+q*2)+u1(i+q*2)*u1(i+q*3)) 33 enddo 34 35 ierr=0 36 end 37!----------------------------------------------------------------------- 38 subroutine apply(ctx,q,u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,u14,& 39& u15,u16,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,v16,ierr) 40 real*8 ctx 41 real*8 u1(1) 42 real*8 u2(1) 43 real*8 u3(1) 44 real*8 u4(1) 45 real*8 v1(1) 46 real*8 v2(1) 47 real*8 du0,du1 48 integer q,ierr 49 50 do i=1,q 51! mass 52 v1(i) = u2(i)*u4(i) 53! diff 54 du0=u1(i+q*0) 55 du1=u1(i+q*1) 56 v2(i+q*0)=u3(i+q*0)*du0+u3(i+q*2)*du1 57 v2(i+q*1)=u3(i+q*2)*du0+u3(i+q*1)*du1 58 enddo 59 60 ierr=0 61 end 62!----------------------------------------------------------------------- 63 program test 64 65 include 'ceedf.h' 66 67 integer ceed,err,i,j,k 68 integer imode 69 parameter(imode=ceed_noninterlaced) 70 integer stridesx(3),stridesu(3),stridesqd(3) 71 integer erestrictx,erestrictu,erestrictxi,erestrictui,erestrictqi 72 integer bx,bu 73 integer qf_setup_mass,qf_setup_diff,qf_apply 74 integer op_setup_mass,op_setup_diff,op_apply 75 integer qdata_mass,qdata_diff,x,a,u,v 76 integer nelem,p,q,d 77 integer row,col,offset 78 parameter(nelem=6) 79 parameter(p=3) 80 parameter(q=4) 81 parameter(d=2) 82 integer ndofs,nqpts,nx,ny 83 parameter(nx=3) 84 parameter(ny=2) 85 parameter(ndofs=(nx*2+1)*(ny*2+1)) 86 parameter(nqpts=nelem*q*q) 87 integer indx(nelem*p*p) 88 real*8 arrx(d*ndofs),aa(nqpts),uu(ndofs),vv(ndofs),atrue(ndofs) 89 integer*8 xoffset,aoffset,uoffset,voffset 90 91 character arg*32 92 93 external setup_mass,setup_diff,apply,apply_lin 94 95 call getarg(1,arg) 96 97 call ceedinit(trim(arg)//char(0),ceed,err) 98 99! DoF Coordinates 100 do i=0,nx*2 101 do j=0,ny*2 102 arrx(i+j*(nx*2+1)+0*ndofs+1)=1.d0*i/(2*nx) 103 arrx(i+j*(nx*2+1)+1*ndofs+1)=1.d0*j/(2*ny) 104 enddo 105 enddo 106 call ceedvectorcreate(ceed,d*ndofs,x,err) 107 xoffset=0 108 call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,arrx,xoffset,err) 109 110! Qdata Vector 111 call ceedvectorcreate(ceed,nqpts,qdata_mass,err) 112 call ceedvectorcreate(ceed,nqpts*d*(d+1)/2,qdata_diff,err) 113 114! Element Setup 115 do i=0,nelem-1 116 col=mod(i,nx) 117 row=i/nx 118 offset=col*(p-1)+row*(nx*2+1)*(p-1) 119 do j=0,p-1 120 do k=0,p-1 121 indx(p*(p*i+k)+j+1)=offset+k*(nx*2+1)+j 122 enddo 123 enddo 124 enddo 125 126! Restrictions 127 call ceedelemrestrictioncreate(ceed,imode,nelem,p*p,ndofs,d,& 128 & ceed_mem_host,ceed_use_pointer,indx,erestrictx,err) 129 stridesx=[1,p*p,p*p*d] 130 call ceedelemrestrictioncreatestrided(ceed,nelem,p*p,& 131 & nelem*p*p,d,stridesx,erestrictxi,err) 132 133 call ceedelemrestrictioncreate(ceed,imode,nelem,p*p,ndofs,1,& 134 & ceed_mem_host,ceed_use_pointer,indx,erestrictu,err) 135 stridesu=[1,q*q,q*q] 136 call ceedelemrestrictioncreatestrided(ceed,nelem,q*q,nqpts,& 137 & 1,stridesu,erestrictui,err) 138 139 stridesqd=[1,q*q,q*q*d*(d+1)/2] 140 call ceedelemrestrictioncreatestrided(ceed,nelem,q*q,nqpts,& 141 & d*(d+1)/2,stridesqd,erestrictqi,err) 142 143! Bases 144 call ceedbasiscreatetensorh1lagrange(ceed,d,d,p,q,ceed_gauss,bx,err) 145 call ceedbasiscreatetensorh1lagrange(ceed,d,1,p,q,ceed_gauss,bu,err) 146 147! QFunction - setup mass 148 call ceedqfunctioncreateinterior(ceed,1,setup_mass,& 149 &SOURCE_DIR& 150 &//'t532-operator.h:setup_mass'//char(0),qf_setup_mass,err) 151 call ceedqfunctionaddinput(qf_setup_mass,'dx',d*d,ceed_eval_grad,err) 152 call ceedqfunctionaddinput(qf_setup_mass,'_weight',1,ceed_eval_weight,err) 153 call ceedqfunctionaddoutput(qf_setup_mass,'qdata',1,ceed_eval_none,err) 154 155! Operator - setup mass 156 call ceedoperatorcreate(ceed,qf_setup_mass,ceed_qfunction_none,& 157 & ceed_qfunction_none,op_setup_mass,err) 158 call ceedoperatorsetfield(op_setup_mass,'dx',erestrictx,& 159 & bx,ceed_vector_active,err) 160 call ceedoperatorsetfield(op_setup_mass,'_weight',erestrictxi,& 161 & bx,ceed_vector_none,err) 162 call ceedoperatorsetfield(op_setup_mass,'qdata',erestrictui,& 163 & ceed_basis_collocated,ceed_vector_active,err) 164 165! QFunction - setup diff 166 call ceedqfunctioncreateinterior(ceed,1,setup_diff,& 167 &SOURCE_DIR& 168 &//'t532-operator.h:setup_diff'//char(0),qf_setup_diff,err) 169 call ceedqfunctionaddinput(qf_setup_diff,'dx',d*d,ceed_eval_grad,err) 170 call ceedqfunctionaddinput(qf_setup_diff,'_weight',1,ceed_eval_weight,err) 171 call ceedqfunctionaddoutput(qf_setup_diff,'qdata',& 172 & d*(d+1)/2,ceed_eval_none,err) 173 174! Operator - setup diff 175 call ceedoperatorcreate(ceed,qf_setup_diff,ceed_qfunction_none,& 176 & ceed_qfunction_none,op_setup_diff,err) 177 call ceedoperatorsetfield(op_setup_diff,'dx',erestrictx,& 178 & bx,ceed_vector_active,err) 179 call ceedoperatorsetfield(op_setup_diff,'_weight',erestrictxi,& 180 & bx,ceed_vector_none,err) 181 call ceedoperatorsetfield(op_setup_diff,'qdata',erestrictqi,& 182 & ceed_basis_collocated,ceed_vector_active,err) 183 184! Apply Setup Operators 185 call ceedoperatorapply(op_setup_mass,x,qdata_mass,& 186 & ceed_request_immediate,err) 187 call ceedoperatorapply(op_setup_diff,x,qdata_diff,& 188 & ceed_request_immediate,err) 189 190! QFunction - apply 191 call ceedqfunctioncreateinterior(ceed,1,apply,& 192 &SOURCE_DIR& 193 &//'t532-operator.h:apply'//char(0),qf_apply,err) 194 call ceedqfunctionaddinput(qf_apply,'du',d,ceed_eval_grad,err) 195 call ceedqfunctionaddinput(qf_apply,'qdata_mass',1,ceed_eval_none,err) 196 call ceedqfunctionaddinput(qf_apply,'qdata_diff',& 197 & d*(d+1)/2,ceed_eval_none,err) 198 call ceedqfunctionaddinput(qf_apply,'u',1,ceed_eval_interp,err) 199 call ceedqfunctionaddoutput(qf_apply,'v',1,ceed_eval_interp,err) 200 call ceedqfunctionaddoutput(qf_apply,'dv',d,ceed_eval_grad,err) 201 202! Operator - apply 203 call ceedoperatorcreate(ceed,qf_apply,ceed_qfunction_none,& 204 & ceed_qfunction_none,op_apply,err) 205 call ceedoperatorsetfield(op_apply,'du',erestrictu,& 206 & bu,ceed_vector_active,err) 207 call ceedoperatorsetfield(op_apply,'qdata_mass',erestrictui,& 208 & ceed_basis_collocated,qdata_mass,err) 209 call ceedoperatorsetfield(op_apply,'qdata_diff',erestrictqi,& 210 & ceed_basis_collocated,qdata_diff,err) 211 call ceedoperatorsetfield(op_apply,'u',erestrictu,& 212 & bu,ceed_vector_active,err) 213 call ceedoperatorsetfield(op_apply,'v',erestrictu,& 214 & bu,ceed_vector_active,err) 215 call ceedoperatorsetfield(op_apply,'dv',erestrictu,& 216 & bu,ceed_vector_active,err) 217 218! Assemble Diagonal 219 call ceedoperatorassemblelineardiagonal(op_apply,a,& 220 & ceed_request_immediate,err) 221 222! Manually assemble diagonal 223 call ceedvectorcreate(ceed,ndofs,u,err) 224 call ceedvectorsetvalue(u,0.d0,err) 225 call ceedvectorcreate(ceed,ndofs,v,err) 226 do i=1,ndofs 227 call ceedvectorgetarray(u,ceed_mem_host,uu,uoffset,err) 228 uu(i+uoffset)=1.d0 229 if (i>1) then 230 uu(i-1+uoffset)=0.d0 231 endif 232 call ceedvectorrestorearray(u,uu,uoffset,err) 233 234 call ceedoperatorapply(op_apply,u,v,ceed_request_immediate,err) 235 236 call ceedvectorgetarrayread(v,ceed_mem_host,vv,voffset,err) 237 atrue(i)=vv(voffset+i) 238 call ceedvectorrestorearrayread(v,vv,voffset,err) 239 enddo 240 241! Check Output 242 call ceedvectorgetarrayread(a,ceed_mem_host,aa,aoffset,err) 243 do i=1,ndofs 244 if (abs(aa(aoffset+i)-atrue(i))>1.0d-14) then 245! LCOV_EXCL_START 246 write(*,*) '[',i,'] Error in assembly: ',aa(aoffset+i),' != ',& 247 & atrue(i) 248! LCOV_EXCL_STOP 249 endif 250 enddo 251 call ceedvectorrestorearrayread(a,aa,aoffset,err) 252 253! Cleanup 254 call ceedqfunctiondestroy(qf_setup_mass,err) 255 call ceedqfunctiondestroy(qf_setup_diff,err) 256 call ceedqfunctiondestroy(qf_apply,err) 257 call ceedoperatordestroy(op_setup_mass,err) 258 call ceedoperatordestroy(op_setup_diff,err) 259 call ceedoperatordestroy(op_apply,err) 260 call ceedelemrestrictiondestroy(erestrictu,err) 261 call ceedelemrestrictiondestroy(erestrictx,err) 262 call ceedelemrestrictiondestroy(erestrictxi,err) 263 call ceedelemrestrictiondestroy(erestrictui,err) 264 call ceedelemrestrictiondestroy(erestrictqi,err) 265 call ceedbasisdestroy(bu,err) 266 call ceedbasisdestroy(bx,err) 267 call ceedvectordestroy(x,err) 268 call ceedvectordestroy(a,err) 269 call ceedvectordestroy(u,err) 270 call ceedvectordestroy(v,err) 271 call ceedvectordestroy(qdata_mass,err) 272 call ceedvectordestroy(qdata_diff,err) 273 call ceeddestroy(ceed,err) 274 end 275!----------------------------------------------------------------------- 276