xref: /libCEED/tests/t505-operator-f.f90 (revision d9b786505a4dfcb66b2fcd9e3b61dd507168515d)
1!-----------------------------------------------------------------------
2!
3! Header with QFunctions
4!
5      include 't500-operator-f.h'
6!-----------------------------------------------------------------------
7      program test
8      implicit none
9      include 'ceed/fortran.h'
10
11      integer ceed,err,i,j
12      integer stridesu(3)
13      integer erestrictx,erestrictu,erestrictui
14      integer bx,bu
15      integer qf_setup,qf_mass
16      integer op_setup,op_mass
17      integer qdata,x,u,v
18      integer nelem,p,q
19      parameter(nelem=15)
20      parameter(p=5)
21      parameter(q=8)
22      integer nx,nu
23      parameter(nx=nelem+1)
24      parameter(nu=nelem*(p-1)+1)
25      integer indx(nelem*2)
26      integer indu(nelem*p)
27      real*8 arrx(nx)
28      integer*8 voffset,xoffset
29
30      real*8 hv(nu)
31      real*8 total
32
33      character arg*32
34
35      external setup,mass
36
37      call getarg(1,arg)
38      call ceedinit(trim(arg)//char(0),ceed,err)
39
40      do i=0,nx-1
41        arrx(i+1)=i/(nx-1.d0)
42      enddo
43      do i=0,nelem-1
44        indx(2*i+1)=i
45        indx(2*i+2)=i+1
46      enddo
47
48      call ceedelemrestrictioncreate(ceed,nelem,2,1,1,nx,ceed_mem_host,&
49     & ceed_use_pointer,indx,erestrictx,err)
50
51      do i=0,nelem-1
52        do j=0,p-1
53          indu(p*i+j+1)=i*(p-1)+j
54        enddo
55      enddo
56
57      call ceedelemrestrictioncreate(ceed,nelem,p,1,1,nu,ceed_mem_host,&
58     & ceed_use_pointer,indu,erestrictu,err)
59      stridesu=[1,q,q]
60      call ceedelemrestrictioncreatestrided(ceed,nelem,q,1,q*nelem,stridesu,&
61     & erestrictui,err)
62
63      call ceedbasiscreatetensorh1lagrange(ceed,1,1,2,q,ceed_gauss,bx,err)
64      call ceedbasiscreatetensorh1lagrange(ceed,1,1,p,q,ceed_gauss,bu,err)
65
66      call ceedqfunctioncreateinterior(ceed,1,setup,&
67     &SOURCE_DIR&
68     &//'t500-operator.h:setup'//char(0),qf_setup,err)
69      call ceedqfunctionaddinput(qf_setup,'weight',1,ceed_eval_weight,err)
70      call ceedqfunctionaddinput(qf_setup,'dx',1,ceed_eval_grad,err)
71      call ceedqfunctionaddoutput(qf_setup,'rho',1,ceed_eval_none,err)
72
73      call ceedqfunctioncreateinterior(ceed,1,mass,&
74     &SOURCE_DIR&
75     &//'t500-operator.h:mass'//char(0),qf_mass,err)
76      call ceedqfunctionaddinput(qf_mass,'rho',1,ceed_eval_none,err)
77      call ceedqfunctionaddinput(qf_mass,'u',1,ceed_eval_interp,err)
78      call ceedqfunctionaddoutput(qf_mass,'v',1,ceed_eval_interp,err)
79
80      call ceedoperatorcreate(ceed,qf_setup,ceed_qfunction_none,&
81     & ceed_qfunction_none,op_setup,err)
82      call ceedoperatorcreate(ceed,qf_mass,ceed_qfunction_none,&
83     & ceed_qfunction_none,op_mass,err)
84
85      call ceedvectorcreate(ceed,nx,x,err)
86      xoffset=0
87      call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,arrx,xoffset,err)
88      call ceedvectorcreate(ceed,nelem*q,qdata,err)
89
90      call ceedoperatorsetfield(op_setup,'weight',ceed_elemrestriction_none,&
91     & bx,ceed_vector_none,err)
92      call ceedoperatorsetfield(op_setup,'dx',erestrictx,bx,&
93     & ceed_vector_active,err)
94      call ceedoperatorsetfield(op_setup,'rho',erestrictui,&
95     ceed_basis_none,ceed_vector_active,err)
96      call ceedoperatorsetfield(op_mass,'rho',erestrictui,&
97     ceed_basis_none,qdata,err)
98      call ceedoperatorsetfield(op_mass,'u',erestrictu,bu,&
99     & ceed_vector_active,err)
100      call ceedoperatorsetfield(op_mass,'v',erestrictu,bu,&
101     & ceed_vector_active,err)
102
103      call ceedoperatorapply(op_setup,x,qdata,ceed_request_immediate,err)
104
105      call ceedvectorcreate(ceed,nu,u,err)
106      call ceedvectorsetvalue(u,1.d0,err)
107      call ceedvectorcreate(ceed,nu,v,err)
108      call ceedvectorsetvalue(v,0.d0,err)
109
110      call ceedoperatorapplyadd(op_mass,u,v,ceed_request_immediate,err)
111
112      call ceedvectorgetarrayread(v,ceed_mem_host,hv,voffset,err)
113      total=0.
114      do i=1,nu
115        total=total+hv(voffset+i)
116      enddo
117      if (abs(total-1.)>1.0d-10) then
118! LCOV_EXCL_START
119        write(*,*) 'Computed Area: ',total,' != True Area: 1.0'
120! LCOV_EXCL_STOP
121      endif
122      call ceedvectorrestorearrayread(v,hv,voffset,err)
123
124      call ceedvectorsetvalue(v,1.d0,err)
125      call ceedoperatorapplyadd(op_mass,u,v,ceed_request_immediate,err)
126
127      call ceedvectorgetarrayread(v,ceed_mem_host,hv,voffset,err)
128      total=-nu
129      do i=1,nu
130        total=total+hv(voffset+i)
131      enddo
132      if (abs(total-1.)>1.0d-10) then
133! LCOV_EXCL_START
134        write(*,*) 'Computed Area: ',total,' != True Area: 1.0'
135! LCOV_EXCL_STOP
136      endif
137      call ceedvectorrestorearrayread(v,hv,voffset,err)
138
139      call ceedvectordestroy(qdata,err)
140      call ceedvectordestroy(x,err)
141      call ceedvectordestroy(u,err)
142      call ceedvectordestroy(v,err)
143      call ceedoperatordestroy(op_mass,err)
144      call ceedoperatordestroy(op_setup,err)
145      call ceedqfunctiondestroy(qf_mass,err)
146      call ceedqfunctiondestroy(qf_setup,err)
147      call ceedbasisdestroy(bu,err)
148      call ceedbasisdestroy(bx,err)
149      call ceedelemrestrictiondestroy(erestrictu,err)
150      call ceedelemrestrictiondestroy(erestrictx,err)
151      call ceedelemrestrictiondestroy(erestrictui,err)
152      call ceeddestroy(ceed,err)
153      end
154!-----------------------------------------------------------------------
155