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