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 ceedvectorcreate(ceed,ndofs,a,err) 215 call ceedoperatorlinearassemblediagonal(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(erestrictui,err) 259 call ceedelemrestrictiondestroy(erestrictqi,err) 260 call ceedbasisdestroy(bu,err) 261 call ceedbasisdestroy(bx,err) 262 call ceedvectordestroy(x,err) 263 call ceedvectordestroy(a,err) 264 call ceedvectordestroy(u,err) 265 call ceedvectordestroy(v,err) 266 call ceedvectordestroy(qdata_mass,err) 267 call ceedvectordestroy(qdata_diff,err) 268 call ceeddestroy(ceed,err) 269 end 270!----------------------------------------------------------------------- 271