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