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 stridesu(3),stridesqd(3) 71 integer erestrictx,erestrictu,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 130 call ceedelemrestrictioncreate(ceed,imode,nelem,p*p,ndofs,1,& 131 & ceed_mem_host,ceed_use_pointer,indx,erestrictu,err) 132 stridesu=[1,q*q,q*q] 133 call ceedelemrestrictioncreatestrided(ceed,nelem,q*q,nqpts,& 134 & 1,stridesu,erestrictui,err) 135 136 stridesqd=[1,q*q,q*q*d*(d+1)/2] 137 call ceedelemrestrictioncreatestrided(ceed,nelem,q*q,nqpts,& 138 & d*(d+1)/2,stridesqd,erestrictqi,err) 139 140! Bases 141 call ceedbasiscreatetensorh1lagrange(ceed,d,d,p,q,ceed_gauss,bx,err) 142 call ceedbasiscreatetensorh1lagrange(ceed,d,1,p,q,ceed_gauss,bu,err) 143 144! QFunction - setup mass 145 call ceedqfunctioncreateinterior(ceed,1,setup_mass,& 146 &SOURCE_DIR& 147 &//'t532-operator.h:setup_mass'//char(0),qf_setup_mass,err) 148 call ceedqfunctionaddinput(qf_setup_mass,'dx',d*d,ceed_eval_grad,err) 149 call ceedqfunctionaddinput(qf_setup_mass,'_weight',1,ceed_eval_weight,err) 150 call ceedqfunctionaddoutput(qf_setup_mass,'qdata',1,ceed_eval_none,err) 151 152! Operator - setup mass 153 call ceedoperatorcreate(ceed,qf_setup_mass,ceed_qfunction_none,& 154 & ceed_qfunction_none,op_setup_mass,err) 155 call ceedoperatorsetfield(op_setup_mass,'dx',erestrictx,& 156 & bx,ceed_vector_active,err) 157 call ceedoperatorsetfield(op_setup_mass,'_weight',& 158 & ceed_elemrestriction_none,bx,ceed_vector_none,err) 159 call ceedoperatorsetfield(op_setup_mass,'qdata',erestrictui,& 160 & ceed_basis_collocated,ceed_vector_active,err) 161 162! QFunction - setup diff 163 call ceedqfunctioncreateinterior(ceed,1,setup_diff,& 164 &SOURCE_DIR& 165 &//'t532-operator.h:setup_diff'//char(0),qf_setup_diff,err) 166 call ceedqfunctionaddinput(qf_setup_diff,'dx',d*d,ceed_eval_grad,err) 167 call ceedqfunctionaddinput(qf_setup_diff,'_weight',1,ceed_eval_weight,err) 168 call ceedqfunctionaddoutput(qf_setup_diff,'qdata',& 169 & d*(d+1)/2,ceed_eval_none,err) 170 171! Operator - setup diff 172 call ceedoperatorcreate(ceed,qf_setup_diff,ceed_qfunction_none,& 173 & ceed_qfunction_none,op_setup_diff,err) 174 call ceedoperatorsetfield(op_setup_diff,'dx',erestrictx,& 175 & bx,ceed_vector_active,err) 176 call ceedoperatorsetfield(op_setup_diff,'_weight',& 177 & ceed_elemrestriction_none,bx,ceed_vector_none,err) 178 call ceedoperatorsetfield(op_setup_diff,'qdata',erestrictqi,& 179 & ceed_basis_collocated,ceed_vector_active,err) 180 181! Apply Setup Operators 182 call ceedoperatorapply(op_setup_mass,x,qdata_mass,& 183 & ceed_request_immediate,err) 184 call ceedoperatorapply(op_setup_diff,x,qdata_diff,& 185 & ceed_request_immediate,err) 186 187! QFunction - apply 188 call ceedqfunctioncreateinterior(ceed,1,apply,& 189 &SOURCE_DIR& 190 &//'t532-operator.h:apply'//char(0),qf_apply,err) 191 call ceedqfunctionaddinput(qf_apply,'du',d,ceed_eval_grad,err) 192 call ceedqfunctionaddinput(qf_apply,'qdata_mass',1,ceed_eval_none,err) 193 call ceedqfunctionaddinput(qf_apply,'qdata_diff',& 194 & d*(d+1)/2,ceed_eval_none,err) 195 call ceedqfunctionaddinput(qf_apply,'u',1,ceed_eval_interp,err) 196 call ceedqfunctionaddoutput(qf_apply,'v',1,ceed_eval_interp,err) 197 call ceedqfunctionaddoutput(qf_apply,'dv',d,ceed_eval_grad,err) 198 199! Operator - apply 200 call ceedoperatorcreate(ceed,qf_apply,ceed_qfunction_none,& 201 & ceed_qfunction_none,op_apply,err) 202 call ceedoperatorsetfield(op_apply,'du',erestrictu,& 203 & bu,ceed_vector_active,err) 204 call ceedoperatorsetfield(op_apply,'qdata_mass',erestrictui,& 205 & ceed_basis_collocated,qdata_mass,err) 206 call ceedoperatorsetfield(op_apply,'qdata_diff',erestrictqi,& 207 & ceed_basis_collocated,qdata_diff,err) 208 call ceedoperatorsetfield(op_apply,'u',erestrictu,& 209 & bu,ceed_vector_active,err) 210 call ceedoperatorsetfield(op_apply,'v',erestrictu,& 211 & bu,ceed_vector_active,err) 212 call ceedoperatorsetfield(op_apply,'dv',erestrictu,& 213 & bu,ceed_vector_active,err) 214 215! Assemble Diagonal 216 call ceedoperatorassemblelineardiagonal(op_apply,a,& 217 & ceed_request_immediate,err) 218 219! Manually assemble diagonal 220 call ceedvectorcreate(ceed,ndofs,u,err) 221 call ceedvectorsetvalue(u,0.d0,err) 222 call ceedvectorcreate(ceed,ndofs,v,err) 223 do i=1,ndofs 224 call ceedvectorgetarray(u,ceed_mem_host,uu,uoffset,err) 225 uu(i+uoffset)=1.d0 226 if (i>1) then 227 uu(i-1+uoffset)=0.d0 228 endif 229 call ceedvectorrestorearray(u,uu,uoffset,err) 230 231 call ceedoperatorapply(op_apply,u,v,ceed_request_immediate,err) 232 233 call ceedvectorgetarrayread(v,ceed_mem_host,vv,voffset,err) 234 atrue(i)=vv(voffset+i) 235 call ceedvectorrestorearrayread(v,vv,voffset,err) 236 enddo 237 238! Check Output 239 call ceedvectorgetarrayread(a,ceed_mem_host,aa,aoffset,err) 240 do i=1,ndofs 241 if (abs(aa(aoffset+i)-atrue(i))>1.0d-14) then 242! LCOV_EXCL_START 243 write(*,*) '[',i,'] Error in assembly: ',aa(aoffset+i),' != ',& 244 & atrue(i) 245! LCOV_EXCL_STOP 246 endif 247 enddo 248 call ceedvectorrestorearrayread(a,aa,aoffset,err) 249 250! Cleanup 251 call ceedqfunctiondestroy(qf_setup_mass,err) 252 call ceedqfunctiondestroy(qf_setup_diff,err) 253 call ceedqfunctiondestroy(qf_apply,err) 254 call ceedoperatordestroy(op_setup_mass,err) 255 call ceedoperatordestroy(op_setup_diff,err) 256 call ceedoperatordestroy(op_apply,err) 257 call ceedelemrestrictiondestroy(erestrictu,err) 258 call ceedelemrestrictiondestroy(erestrictx,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