1!----------------------------------------------------------------------- 2 subroutine setup(ctx,q,u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,u14,& 3& u15,u16,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,v16,ierr) 4 real*8 ctx 5 real*8 u1(1) 6 real*8 u2(1) 7 real*8 v1(1) 8 integer q,ierr 9 10 do i=1,q 11 v1(i)=u1(i)*u2(i) 12 enddo 13 14 ierr=0 15 end 16!----------------------------------------------------------------------- 17 subroutine mass(ctx,q,u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,u14,& 18& u15,u16,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,v16,ierr) 19 real*8 ctx 20 real*8 u1(1) 21 real*8 u2(1) 22 real*8 v1(1) 23 integer q,ierr 24 25 do i=1,q 26 v1(i)=u2(i)*u1(i) 27 v1(q+i)=u2(q+i)*u1(i) 28 enddo 29 30 ierr=0 31 end 32!----------------------------------------------------------------------- 33 program test 34 35 include 'ceedf.h' 36 37 integer ceed,err,i,j 38 integer stridesu(3) 39 integer erestrictx,erestrictui 40 integer erestrictucoarse,erestrictufine 41 integer bx,bufine,bucoarse,bctof 42 integer qf_setup,qf_mass 43 integer op_setup,op_masscoarse,op_massfine 44 integer op_prolong,op_restrict 45 integer qdata,x,ucoarse,ufine,vcoarse,vfine,pmultfine 46 integer nelem,pcoarse,pfine,q,ncomp 47 parameter(ncomp=2) 48 parameter(nelem=15) 49 parameter(pcoarse=3) 50 parameter(pfine=5) 51 parameter(q=8) 52 integer nx,nucoarse,nufine 53 parameter(nx=nelem+1) 54 parameter(nucoarse=nelem*(pcoarse-1)+1) 55 parameter(nufine=nelem*(pfine-1)+1) 56 integer indx(nelem*2) 57 integer inducoarse(nelem*pcoarse) 58 integer indufine(nelem*pfine) 59 real*8 arrx(nx) 60 integer*8 voffset,xoffset,ioffset 61 real*8 val 62 integer interpsize 63 parameter(interpsize=pcoarse*pfine); 64 real*8 interp1d(interpsize),interpctof(interpsize) 65 66 real*8 hv(nufine*ncomp) 67 real*8 total 68 69 character arg*32 70 71 external setup,mass 72 73 call getarg(1,arg) 74 call ceedinit(trim(arg)//char(0),ceed,err) 75 76 do i=0,nx-1 77 arrx(i+1)=i/(nx-1.d0) 78 enddo 79 do i=0,nelem-1 80 indx(2*i+1)=i 81 indx(2*i+2)=i+1 82 enddo 83 call ceedelemrestrictioncreate(ceed,nelem,2,1,1,nx,ceed_mem_host,& 84 & ceed_use_pointer,indx,erestrictx,err) 85 86 do i=0,nelem-1 87 do j=0,pcoarse-1 88 inducoarse(pcoarse*i+j+1)=i*(pcoarse-1)+j 89 enddo 90 enddo 91 call ceedelemrestrictioncreate(ceed,nelem,pcoarse,ncomp,nucoarse,& 92 & ncomp*nucoarse,ceed_mem_host,ceed_use_pointer,inducoarse,& 93 & erestrictucoarse,err) 94 95 do i=0,nelem-1 96 do j=0,pfine-1 97 indufine(pfine*i+j+1)=i*(pfine-1)+j 98 enddo 99 enddo 100 call ceedelemrestrictioncreate(ceed,nelem,pfine,ncomp,nufine,& 101 & ncomp*nufine,ceed_mem_host,ceed_use_pointer,indufine,& 102 & erestrictufine,err) 103 104 stridesu=[1,q,q] 105 call ceedelemrestrictioncreatestrided(ceed,nelem,q,1,q*nelem,stridesu,& 106 & erestrictui,err) 107 108 call ceedbasiscreatetensorh1lagrange(ceed,1,1,2,q,ceed_gauss,bx,err) 109 call ceedbasiscreatetensorh1lagrange(ceed,1,ncomp,pfine,q,ceed_gauss,& 110 & bufine,err) 111 112 call ceedqfunctioncreateinterior(ceed,1,setup,& 113 &SOURCE_DIR& 114 &//'t502-operator.h:setup'//char(0),qf_setup,err) 115 call ceedqfunctionaddinput(qf_setup,'weights',1,ceed_eval_weight,err) 116 call ceedqfunctionaddinput(qf_setup,'dx',1,ceed_eval_grad,err) 117 call ceedqfunctionaddoutput(qf_setup,'qdata',1,ceed_eval_none,err) 118 119 call ceedqfunctioncreateinterior(ceed,1,mass,& 120 &SOURCE_DIR& 121 &//'t502-operator.h:mass'//char(0),qf_mass,err) 122 call ceedqfunctionaddinput(qf_mass,'qdata',1,ceed_eval_none,err) 123 call ceedqfunctionaddinput(qf_mass,'u',ncomp,ceed_eval_interp,err) 124 call ceedqfunctionaddoutput(qf_mass,'v',ncomp,ceed_eval_interp,err) 125 126 call ceedoperatorcreate(ceed,qf_setup,ceed_qfunction_none,& 127 & ceed_qfunction_none,op_setup,err) 128 call ceedoperatorcreate(ceed,qf_mass,ceed_qfunction_none,& 129 & ceed_qfunction_none,op_massfine,err) 130 131 call ceedvectorcreate(ceed,nx,x,err) 132 xoffset=0 133 call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,arrx,xoffset,err) 134 call ceedvectorcreate(ceed,nelem*q,qdata,err) 135 136 call ceedoperatorsetfield(op_setup,'weights',ceed_elemrestriction_none,& 137 & bx,ceed_vector_none,err) 138 call ceedoperatorsetfield(op_setup,'dx',erestrictx,bx,& 139 & ceed_vector_active,err) 140 call ceedoperatorsetfield(op_setup,'qdata',erestrictui,& 141 & ceed_basis_collocated,ceed_vector_active,err) 142 call ceedoperatorsetfield(op_massfine,'qdata',erestrictui,& 143 & ceed_basis_collocated,qdata,err) 144 call ceedoperatorsetfield(op_massfine,'u',erestrictufine,bufine,& 145 & ceed_vector_active,err) 146 call ceedoperatorsetfield(op_massfine,'v',erestrictufine,bufine,& 147 & ceed_vector_active,err) 148 149 call ceedoperatorapply(op_setup,x,qdata,ceed_request_immediate,err) 150 151! Create multigrid level 152 call ceedvectorcreate(ceed,ncomp*nufine,pmultfine,err) 153 val=1.0 154 call ceedvectorsetvalue(pmultfine,val,err) 155 call ceedbasiscreatetensorh1lagrange(ceed,1,ncomp,pcoarse,q,& 156 & ceed_gauss,bucoarse,err) 157 call ceedbasiscreatetensorh1lagrange(ceed,1,ncomp,pcoarse,pfine,& 158 & ceed_gauss_lobatto,bctof,err) 159 call ceedbasisgetinterp1d(bctof,interp1d,ioffset,err) 160 do i=1,interpsize 161 interpctof(i)=interp1d(ioffset+i) 162 enddo 163 call ceedoperatormultigridlevelcreateh1(op_massfine,pmultfine,& 164 & erestrictucoarse,bucoarse,interpctof,op_masscoarse,& 165 & op_prolong,op_restrict,err) 166 167! Coarse problem 168 call ceedvectorcreate(ceed,ncomp*nucoarse,ucoarse,err) 169 val=1.0 170 call ceedvectorsetvalue(ucoarse,val,err) 171 call ceedvectorcreate(ceed,ncomp*nucoarse,vcoarse,err) 172 call ceedoperatorapply(op_masscoarse,ucoarse,vcoarse,& 173 & ceed_request_immediate,err) 174 175! Check output 176 call ceedvectorgetarrayread(vcoarse,ceed_mem_host,hv,voffset,err) 177 total=0. 178 do i=1,nucoarse*ncomp 179 total=total+hv(voffset+i) 180 enddo 181 if (abs(total-2.)>1.0d-10) then 182! LCOV_EXCL_START 183 write(*,*) 'Computed Area: ',total,' != True Area: 1.0' 184! LCOV_EXCL_STOP 185 endif 186 call ceedvectorrestorearrayread(vcoarse,hv,voffset,err) 187 188! Prolong coarse u 189 call ceedvectorcreate(ceed,ncomp*nufine,ufine,err) 190 call ceedoperatorapply(op_prolong,ucoarse,ufine,& 191 & ceed_request_immediate,err) 192 193! Fine problem 194 call ceedvectorcreate(ceed,ncomp*nufine,vfine,err) 195 call ceedoperatorapply(op_massfine,ufine,vfine,& 196 & ceed_request_immediate,err) 197 198! Check output 199 call ceedvectorgetarrayread(vfine,ceed_mem_host,hv,voffset,err) 200 total=0. 201 do i=1,nufine*ncomp 202 total=total+hv(voffset+i) 203 enddo 204 if (abs(total-2.)>1.0d-10) then 205! LCOV_EXCL_START 206 write(*,*) 'Computed Area Fine Grid: ',total,' != True Area: 1.0' 207! LCOV_EXCL_STOP 208 endif 209 call ceedvectorrestorearrayread(vfine,hv,voffset,err) 210 211! Restrict state to coarse grid 212 call ceedoperatorapply(op_restrict,vfine,vcoarse,& 213 & ceed_request_immediate,err) 214 215! Check output 216 call ceedvectorgetarrayread(vcoarse,ceed_mem_host,hv,voffset,err) 217 total=0. 218 do i=1,nucoarse*ncomp 219 total=total+hv(voffset+i) 220 enddo 221 if (abs(total-2.)>1.0d-10) then 222! LCOV_EXCL_START 223 write(*,*) 'Computed Area Coarse Grid: ',total,' != True Area: 1.0' 224! LCOV_EXCL_STOP 225 endif 226 call ceedvectorrestorearrayread(vcoarse,hv,voffset,err) 227 228 call ceedvectordestroy(qdata,err) 229 call ceedvectordestroy(x,err) 230 call ceedvectordestroy(ucoarse,err) 231 call ceedvectordestroy(ufine,err) 232 call ceedvectordestroy(vcoarse,err) 233 call ceedvectordestroy(vfine,err) 234 call ceedvectordestroy(pmultfine,err) 235 call ceedoperatordestroy(op_masscoarse,err) 236 call ceedoperatordestroy(op_massfine,err) 237 call ceedoperatordestroy(op_prolong,err) 238 call ceedoperatordestroy(op_restrict,err) 239 call ceedoperatordestroy(op_setup,err) 240 call ceedqfunctiondestroy(qf_mass,err) 241 call ceedqfunctiondestroy(qf_setup,err) 242 call ceedbasisdestroy(bufine,err) 243 call ceedbasisdestroy(bx,err) 244 call ceedelemrestrictiondestroy(erestrictucoarse,err) 245 call ceedelemrestrictiondestroy(erestrictufine,err) 246 call ceedelemrestrictiondestroy(erestrictx,err) 247 call ceedelemrestrictiondestroy(erestrictui,err) 248 call ceeddestroy(ceed,err) 249 end 250!----------------------------------------------------------------------- 251