1!----------------------------------------------------------------------- 2! 3! Header with QFunctions 4! 5 include 't535-operator-f.h' 6!----------------------------------------------------------------------- 7 program test 8 implicit none 9 include 'ceed/fortran.h' 10 11 integer ceed,err,i,j,k 12 integer stridesu(3),stridesqd(3) 13 integer erestrictx,erestrictu,erestrictui,erestrictqi 14 integer bx,bu 15 integer qf_setup_mass,qf_setup_diff,qf_apply 16 integer op_setup_mass,op_setup_diff,op_apply 17 integer qdata_mass,qdata_diff,x,a,u,v 18 integer nelem,p,q,d 19 integer row,col,offset 20 parameter(nelem=6) 21 parameter(p=3) 22 parameter(q=4) 23 parameter(d=2) 24 integer ndofs,nqpts,nx,ny 25 parameter(nx=3) 26 parameter(ny=2) 27 parameter(ndofs=(nx*2+1)*(ny*2+1)) 28 parameter(nqpts=nelem*q*q) 29 integer indx(nelem*p*p) 30 real*8 arrx(d*ndofs),aa(nqpts),uu(ndofs),vv(ndofs),atrue(ndofs) 31 integer*8 xoffset,aoffset,uoffset,voffset 32 33 character arg*32 34 35 external setup_mass,setup_diff,apply,apply_lin 36 37 call getarg(1,arg) 38 39 call ceedinit(trim(arg)//char(0),ceed,err) 40 41! DoF Coordinates 42 do i=0,nx*2 43 do j=0,ny*2 44 arrx(i+j*(nx*2+1)+0*ndofs+1)=1.d0*i/(2*nx) 45 arrx(i+j*(nx*2+1)+1*ndofs+1)=1.d0*j/(2*ny) 46 enddo 47 enddo 48 call ceedvectorcreate(ceed,d*ndofs,x,err) 49 xoffset=0 50 call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,arrx,xoffset,err) 51 52! Qdata Vector 53 call ceedvectorcreate(ceed,nqpts,qdata_mass,err) 54 call ceedvectorcreate(ceed,nqpts*d*(d+1)/2,qdata_diff,err) 55 56! Element Setup 57 do i=0,nelem-1 58 col=mod(i,nx) 59 row=i/nx 60 offset=col*(p-1)+row*(nx*2+1)*(p-1) 61 do j=0,p-1 62 do k=0,p-1 63 indx(p*(p*i+k)+j+1)=offset+k*(nx*2+1)+j 64 enddo 65 enddo 66 enddo 67 68! Restrictions 69 call ceedelemrestrictioncreate(ceed,nelem,p*p,d,ndofs,d*ndofs,& 70 & ceed_mem_host,ceed_use_pointer,indx,erestrictx,err) 71 72 call ceedelemrestrictioncreate(ceed,nelem,p*p,1,1,ndofs,& 73 & ceed_mem_host,ceed_use_pointer,indx,erestrictu,err) 74 stridesu=[1,q*q,q*q] 75 call ceedelemrestrictioncreatestrided(ceed,nelem,q*q,1,nqpts,& 76 & stridesu,erestrictui,err) 77 78 stridesqd=[1,q*q,q*q*d*(d+1)/2] 79 call ceedelemrestrictioncreatestrided(ceed,nelem,q*q,d*(d+1)/2,& 80 & d*(d+1)/2*nqpts,stridesqd,erestrictqi,err) 81 82! Bases 83 call ceedbasiscreatetensorh1lagrange(ceed,d,d,p,q,ceed_gauss,bx,err) 84 call ceedbasiscreatetensorh1lagrange(ceed,d,1,p,q,ceed_gauss,bu,err) 85 86! QFunction - setup mass 87 call ceedqfunctioncreateinterior(ceed,1,setup_mass,& 88 &SOURCE_DIR& 89 &//'t532-operator.h:setup_mass'//char(0),qf_setup_mass,err) 90 call ceedqfunctionaddinput(qf_setup_mass,'dx',d*d,ceed_eval_grad,err) 91 call ceedqfunctionaddinput(qf_setup_mass,'weight',1,ceed_eval_weight,err) 92 call ceedqfunctionaddoutput(qf_setup_mass,'qdata',1,ceed_eval_none,err) 93 94! Operator - setup mass 95 call ceedoperatorcreate(ceed,qf_setup_mass,ceed_qfunction_none,& 96 & ceed_qfunction_none,op_setup_mass,err) 97 call ceedoperatorsetfield(op_setup_mass,'dx',erestrictx,& 98 & bx,ceed_vector_active,err) 99 call ceedoperatorsetfield(op_setup_mass,'weight',& 100 & ceed_elemrestriction_none,bx,ceed_vector_none,err) 101 call ceedoperatorsetfield(op_setup_mass,'qdata',erestrictui,& 102 ceed_basis_none,ceed_vector_active,err) 103 104! QFunction - setup diff 105 call ceedqfunctioncreateinterior(ceed,1,setup_diff,& 106 &SOURCE_DIR& 107 &//'t532-operator.h:setup_diff'//char(0),qf_setup_diff,err) 108 call ceedqfunctionaddinput(qf_setup_diff,'dx',d*d,ceed_eval_grad,err) 109 call ceedqfunctionaddinput(qf_setup_diff,'weight',1,ceed_eval_weight,err) 110 call ceedqfunctionaddoutput(qf_setup_diff,'qdata',& 111 & d*(d+1)/2,ceed_eval_none,err) 112 113! Operator - setup diff 114 call ceedoperatorcreate(ceed,qf_setup_diff,ceed_qfunction_none,& 115 & ceed_qfunction_none,op_setup_diff,err) 116 call ceedoperatorsetfield(op_setup_diff,'dx',erestrictx,& 117 & bx,ceed_vector_active,err) 118 call ceedoperatorsetfield(op_setup_diff,'weight',& 119 & ceed_elemrestriction_none,bx,ceed_vector_none,err) 120 call ceedoperatorsetfield(op_setup_diff,'qdata',erestrictqi,& 121 ceed_basis_none,ceed_vector_active,err) 122 123! Apply Setup Operators 124 call ceedoperatorapply(op_setup_mass,x,qdata_mass,& 125 & ceed_request_immediate,err) 126 call ceedoperatorapply(op_setup_diff,x,qdata_diff,& 127 & ceed_request_immediate,err) 128 129! QFunction - apply 130 call ceedqfunctioncreateinterior(ceed,1,apply,& 131 &SOURCE_DIR& 132 &//'t532-operator.h:apply'//char(0),qf_apply,err) 133 call ceedqfunctionaddinput(qf_apply,'du',d,ceed_eval_grad,err) 134 call ceedqfunctionaddinput(qf_apply,'mass qdata',1,ceed_eval_none,err) 135 call ceedqfunctionaddinput(qf_apply,'diff qdata',& 136 & d*(d+1)/2,ceed_eval_none,err) 137 call ceedqfunctionaddinput(qf_apply,'u',1,ceed_eval_interp,err) 138 call ceedqfunctionaddoutput(qf_apply,'v',1,ceed_eval_interp,err) 139 call ceedqfunctionaddoutput(qf_apply,'dv',d,ceed_eval_grad,err) 140 141! Operator - apply 142 call ceedoperatorcreate(ceed,qf_apply,ceed_qfunction_none,& 143 & ceed_qfunction_none,op_apply,err) 144 call ceedoperatorsetfield(op_apply,'du',erestrictu,& 145 & bu,ceed_vector_active,err) 146 call ceedoperatorsetfield(op_apply,'mass qdata',erestrictui,& 147 ceed_basis_none,qdata_mass,err) 148 call ceedoperatorsetfield(op_apply,'diff qdata',erestrictqi,& 149 ceed_basis_none,qdata_diff,err) 150 call ceedoperatorsetfield(op_apply,'u',erestrictu,& 151 & bu,ceed_vector_active,err) 152 call ceedoperatorsetfield(op_apply,'v',erestrictu,& 153 & bu,ceed_vector_active,err) 154 call ceedoperatorsetfield(op_apply,'dv',erestrictu,& 155 & bu,ceed_vector_active,err) 156 157! Assemble Diagonal 158 call ceedvectorcreate(ceed,ndofs,a,err) 159 call ceedoperatorlinearassemblediagonal(op_apply,a,& 160 & ceed_request_immediate,err) 161 162! Manually assemble diagonal 163 call ceedvectorcreate(ceed,ndofs,u,err) 164 call ceedvectorsetvalue(u,0.d0,err) 165 call ceedvectorcreate(ceed,ndofs,v,err) 166 do i=1,ndofs 167 call ceedvectorgetarray(u,ceed_mem_host,uu,uoffset,err) 168 uu(i+uoffset)=1.d0 169 if (i>1) then 170 uu(i-1+uoffset)=0.d0 171 endif 172 call ceedvectorrestorearray(u,uu,uoffset,err) 173 174 call ceedoperatorapply(op_apply,u,v,ceed_request_immediate,err) 175 176 call ceedvectorgetarrayread(v,ceed_mem_host,vv,voffset,err) 177 atrue(i)=vv(voffset+i) 178 call ceedvectorrestorearrayread(v,vv,voffset,err) 179 enddo 180 181! Check Output 182 call ceedvectorgetarrayread(a,ceed_mem_host,aa,aoffset,err) 183 do i=1,ndofs 184 if (abs(aa(aoffset+i)-atrue(i))>1.0d-14) then 185! LCOV_EXCL_START 186 write(*,*) '[',i,'] Error in assembly: ',aa(aoffset+i),' != ',& 187 & atrue(i) 188! LCOV_EXCL_STOP 189 endif 190 enddo 191 call ceedvectorrestorearrayread(a,aa,aoffset,err) 192 193! Cleanup 194 call ceedqfunctiondestroy(qf_setup_mass,err) 195 call ceedqfunctiondestroy(qf_setup_diff,err) 196 call ceedqfunctiondestroy(qf_apply,err) 197 call ceedoperatordestroy(op_setup_mass,err) 198 call ceedoperatordestroy(op_setup_diff,err) 199 call ceedoperatordestroy(op_apply,err) 200 call ceedelemrestrictiondestroy(erestrictu,err) 201 call ceedelemrestrictiondestroy(erestrictx,err) 202 call ceedelemrestrictiondestroy(erestrictui,err) 203 call ceedelemrestrictiondestroy(erestrictqi,err) 204 call ceedbasisdestroy(bu,err) 205 call ceedbasisdestroy(bx,err) 206 call ceedvectordestroy(x,err) 207 call ceedvectordestroy(a,err) 208 call ceedvectordestroy(u,err) 209 call ceedvectordestroy(v,err) 210 call ceedvectordestroy(qdata_mass,err) 211 call ceedvectordestroy(qdata_diff,err) 212 call ceeddestroy(ceed,err) 213 end 214!----------------------------------------------------------------------- 215