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