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