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