xref: /libCEED/tests/t550-operator-f.f90 (revision d99fa3c5cd91a1690aedf0679cbf290d44fec74c)
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(0*q+i)=u2(0*q+i)*u1(i)
27        v1(1*q+i)=u2(1*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,bucoarse,bufine
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
61      real*8 val
62
63      real*8 hv(nufine*ncomp)
64      real*8 total
65
66      character arg*32
67
68      external setup,mass
69
70      call getarg(1,arg)
71      call ceedinit(trim(arg)//char(0),ceed,err)
72
73      do i=0,nx-1
74        arrx(i+1)=i/(nx-1.d0)
75      enddo
76      do i=0,nelem-1
77        indx(2*i+1)=i
78        indx(2*i+2)=i+1
79      enddo
80      call ceedelemrestrictioncreate(ceed,nelem,2,1,1,nx,ceed_mem_host,&
81     & ceed_use_pointer,indx,erestrictx,err)
82
83      do i=0,nelem-1
84        do j=0,pcoarse-1
85          inducoarse(pcoarse*i+j+1)=i*(pcoarse-1)+j
86        enddo
87      enddo
88      call ceedelemrestrictioncreate(ceed,nelem,pcoarse,ncomp,nucoarse,&
89     & ncomp*nucoarse,ceed_mem_host,ceed_use_pointer,inducoarse,&
90     & erestrictucoarse,err)
91
92      do i=0,nelem-1
93        do j=0,pfine-1
94          indufine(pfine*i+j+1)=i*(pfine-1)+j
95        enddo
96      enddo
97      call ceedelemrestrictioncreate(ceed,nelem,pfine,ncomp,nufine,&
98     & ncomp*nufine,ceed_mem_host,ceed_use_pointer,indufine,&
99     & erestrictufine,err)
100
101     stridesu=[1,q,q]
102      call ceedelemrestrictioncreatestrided(ceed,nelem,q,1,q*nelem,stridesu,&
103     & erestrictui,err)
104
105      call ceedbasiscreatetensorh1lagrange(ceed,1,1,2,q,ceed_gauss,bx,err)
106      call ceedbasiscreatetensorh1lagrange(ceed,1,ncomp,pfine,q,ceed_gauss,&
107     & bufine,err)
108      call ceedbasiscreatetensorh1lagrange(ceed,1,ncomp,pcoarse,q,ceed_gauss,&
109     & bucoarse,err)
110
111     call ceedqfunctioncreateinterior(ceed,1,setup,&
112    &SOURCE_DIR&
113    &//'t502-operator.h:setup'//char(0),qf_setup,err)
114     call ceedqfunctionaddinput(qf_setup,'weights',1,ceed_eval_weight,err)
115     call ceedqfunctionaddinput(qf_setup,'dx',1,ceed_eval_grad,err)
116     call ceedqfunctionaddoutput(qf_setup,'qdata',1,ceed_eval_none,err)
117
118     call ceedqfunctioncreateinterior(ceed,1,mass,&
119    &SOURCE_DIR&
120    &//'t502-operator.h:mass'//char(0),qf_mass,err)
121     call ceedqfunctionaddinput(qf_mass,'qdata',1,ceed_eval_none,err)
122     call ceedqfunctionaddinput(qf_mass,'u',ncomp,ceed_eval_interp,err)
123     call ceedqfunctionaddoutput(qf_mass,'v',ncomp,ceed_eval_interp,err)
124
125      call ceedoperatorcreate(ceed,qf_setup,ceed_qfunction_none,&
126     & ceed_qfunction_none,op_setup,err)
127      call ceedoperatorcreate(ceed,qf_mass,ceed_qfunction_none,&
128     & ceed_qfunction_none,op_massfine,err)
129
130      call ceedvectorcreate(ceed,nx,x,err)
131      xoffset=0
132      call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,arrx,xoffset,err)
133      call ceedvectorcreate(ceed,nelem*q,qdata,err)
134
135      call ceedoperatorsetfield(op_setup,'weights',ceed_elemrestriction_none,&
136     & bx,ceed_vector_none,err)
137      call ceedoperatorsetfield(op_setup,'dx',erestrictx,bx,&
138     & ceed_vector_active,err)
139      call ceedoperatorsetfield(op_setup,'qdata',erestrictui,&
140     & ceed_basis_collocated,ceed_vector_active,err)
141      call ceedoperatorsetfield(op_massfine,'qdata',erestrictui,&
142     & ceed_basis_collocated,qdata,err)
143      call ceedoperatorsetfield(op_massfine,'u',erestrictufine,bufine,&
144     & ceed_vector_active,err)
145      call ceedoperatorsetfield(op_massfine,'v',erestrictufine,bufine,&
146     & ceed_vector_active,err)
147
148      call ceedoperatorapply(op_setup,x,qdata,ceed_request_immediate,err)
149
150! Create multigrid level
151      call ceedvectorcreate(ceed,ncomp*nufine,pmultfine,err)
152      val=1.0
153      call ceedvectorsetvalue(pmultfine,val,err)
154      call ceedoperatormultigridlevelcreate(op_massfine,pmultfine,&
155     & erestrictucoarse,bucoarse,op_masscoarse,op_prolong,op_restrict,err)
156
157! Coarse problem
158      call ceedvectorcreate(ceed,ncomp*nucoarse,ucoarse,err)
159      val=1.0
160      call ceedvectorsetvalue(ucoarse,val,err)
161      call ceedvectorcreate(ceed,ncomp*nucoarse,vcoarse,err)
162      call ceedoperatorapply(op_masscoarse,ucoarse,vcoarse,&
163     & ceed_request_immediate,err)
164
165! Check output
166      call ceedvectorgetarrayread(vcoarse,ceed_mem_host,hv,voffset,err)
167      total=0.
168      do i=1,nucoarse*ncomp
169        total=total+hv(voffset+i)
170      enddo
171      if (abs(total-2.)>1.0d-10) then
172! LCOV_EXCL_START
173        write(*,*) 'Computed Area: ',total,' != True Area: 1.0'
174! LCOV_EXCL_STOP
175      endif
176      call ceedvectorrestorearrayread(vcoarse,hv,voffset,err)
177
178! Prolong coarse u
179      call ceedvectorcreate(ceed,ncomp*nufine,ufine,err)
180      call ceedoperatorapply(op_prolong,ucoarse,ufine,&
181     & ceed_request_immediate,err)
182
183! Fine problem
184      call ceedvectorcreate(ceed,ncomp*nufine,vfine,err)
185      call ceedoperatorapply(op_massfine,ufine,vfine,&
186     & ceed_request_immediate,err)
187
188! Check output
189      call ceedvectorgetarrayread(vfine,ceed_mem_host,hv,voffset,err)
190      total=0.
191      do i=1,nufine*ncomp
192        total=total+hv(voffset+i)
193      enddo
194      if (abs(total-2.)>1.0d-10) then
195! LCOV_EXCL_START
196        write(*,*) 'Computed Area Fine Grid: ',total,' != True Area: 1.0'
197! LCOV_EXCL_STOP
198      endif
199      call ceedvectorrestorearrayread(vfine,hv,voffset,err)
200
201! Restrict state to coarse grid
202      call ceedoperatorapply(op_restrict,vfine,vcoarse,&
203     & ceed_request_immediate,err)
204
205! Check output
206      call ceedvectorgetarrayread(vcoarse,ceed_mem_host,hv,voffset,err)
207      total=0.
208      do i=1,nucoarse*ncomp
209        total=total+hv(voffset+i)
210      enddo
211      if (abs(total-2.)>1.0d-10) then
212! LCOV_EXCL_START
213        write(*,*) 'Computed Area Coarse Grid: ',total,' != True Area: 1.0'
214! LCOV_EXCL_STOP
215      endif
216      call ceedvectorrestorearrayread(vcoarse,hv,voffset,err)
217
218      call ceedvectordestroy(qdata,err)
219      call ceedvectordestroy(x,err)
220      call ceedvectordestroy(ucoarse,err)
221      call ceedvectordestroy(ufine,err)
222      call ceedvectordestroy(vcoarse,err)
223      call ceedvectordestroy(vfine,err)
224      call ceedvectordestroy(pmultfine,err)
225      call ceedoperatordestroy(op_masscoarse,err)
226      call ceedoperatordestroy(op_massfine,err)
227      call ceedoperatordestroy(op_prolong,err)
228      call ceedoperatordestroy(op_restrict,err)
229      call ceedoperatordestroy(op_setup,err)
230      call ceedqfunctiondestroy(qf_mass,err)
231      call ceedqfunctiondestroy(qf_setup,err)
232      call ceedbasisdestroy(bucoarse,err)
233      call ceedbasisdestroy(bufine,err)
234      call ceedbasisdestroy(bx,err)
235      call ceedelemrestrictiondestroy(erestrictucoarse,err)
236      call ceedelemrestrictiondestroy(erestrictufine,err)
237      call ceedelemrestrictiondestroy(erestrictx,err)
238      call ceedelemrestrictiondestroy(erestrictui,err)
239      call ceeddestroy(ceed,err)
240      end
241!-----------------------------------------------------------------------
242