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