ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
SOLVERFEM.f90
Go to the documentation of this file.
1 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
8 ! 19.01.2017.
9 ! + + + + + + + + + NUMERICAL SOLUTION + + + + + + + + +
10 
11 SUBROUTINE solutionfem (SOLVER, ifail)
12 
13 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
14 ! +
15 ! This subroutine is prepared to solve single +
16 ! transport equation in standardised form +
17 ! adopted by the ETS. +
18 ! +
19 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
20 ! Source: FESB Team +
21 ! Developers: Anna Šušnjara, Vicko Dorić +
22 ! +
23 ! Comments: finite element method based solver +
24 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
25 
26 
27  USE type_solver
28  USE itm_types
29 
30  IMPLICIT NONE
31 
32 ! +++ Input/Output with ETS:
33  INTEGER, INTENT (INOUT) :: ifail
34  TYPE (numerics), INTENT (INOUT) :: solver !contains all I/O quantities to numerics part
35 
36 !#ifdef WANTCOS
37 
38 ! +++ Internal input/output parameters:
39  INTEGER :: idim, ndim !equation index and total number of equations to be solved
40  INTEGER :: irho, nrho !radius index, number of radial points
41 
42  REAL (R8) :: rho(solver%nrho) !radii
43  !REAL (R8) :: RHOMAX
44 
45  REAL (R8) :: amix !fraction of new sollution mixed
46 
47  REAL (R8) :: y(solver%nrho), ym(solver%nrho) !function at the current amd previous time steps
48  REAL (R8) :: dy(solver%nrho) !derivative of function
49 
50  REAL (R8) :: yy(solver%nrho) !new function
51 
52  REAL (R8) :: a(solver%nrho), b(solver%nrho) !coefficients
53  REAL (R8) :: c(solver%nrho), d(solver%nrho) !coefficients
54  REAL (R8) :: e(solver%nrho), f(solver%nrho) !coefficients
55  REAL (R8) :: g(solver%nrho), h !coefficients
56  REAL (R8) :: dd(solver%nrho), de(solver%nrho) !derivatives of coefficients
57 
58  REAL (R8) :: v(2), u(2), w(2) !boundary conditions
59 
60 ! +++ Coefficients used by internal numerical solver:
61  !REAL (R8) :: X(SOLVER%NRHO) !normalized radii
62  !REAL (R8) :: RUB2, RUB1 !radii of inner and outer boundaries
63 
64 
65 
66 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
67 ! +++ Set up local variables (input, obtained from ETS):
68 ! Control parameters:
69  ndim = solver%NDIM
70  nrho = solver%NRHO
71  amix = solver%AMIX
72  h = solver%H
73 
74 
75 ! +++ Solution of equations starting from 1 to NDIM:
76  equation_loop1: DO idim = 1, ndim
77  IF (solver%EQ_FLAG(idim).EQ.0) goto 20 !equation is not solved
78 
79 
80 
81 ! +++ Numerical coefficients obtained from physics part in form:
82 ! (A*Y-B*Y(t-1))/H + 1/C * (-D*Y' + E*Y)' = F - G*Y
83 
84  ! RHOMAX = SOLVER%RHO(NRHO)
85 
86  rho_loop1: DO irho=1,nrho
87  rho(irho) = solver%RHO(irho)
88 
89  y(irho) = solver%Y(idim,irho)
90  dy(irho) = solver%DY(idim,irho)
91  ym(irho) = solver%YM(idim,irho)
92 
93  a(irho) = solver%A(idim,irho)
94  b(irho) = solver%B(idim,irho)
95  c(irho) = solver%C(idim,irho)
96  d(irho) = solver%D(idim,irho)
97  e(irho) = solver%E(idim,irho)
98  f(irho) = solver%F(idim,irho)
99  g(irho) = solver%G(idim,irho)
100  END DO rho_loop1
101 
102  !X = RHO/RHOMAX
103 
104  !CALL DERIV_HOLOB(NRHO,X,D,DD)
105  !CALL DERIV_HOLOB(NRHO,X,E,DE)
106 
107  CALL deriv_holob(nrho,rho,d,dd)
108  CALL deriv_holob(nrho,rho,e,de)
109 
110 
111 ! +++ Boundary conditions in form:
112 ! V*Y' + U*Y = W
113 
114  u = solver%U(idim,:)
115  v = solver%V(idim,:)
116  w = solver%W(idim,:)
117 
118 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
119 ! +++ Solution of transport equation:
120 
121  CALL fem_pde_solver(a,b,c,d,dd,e,de,f,g,h, &
122  ym,yy,rho,nrho,u, v, w)
123  !CALL fem_pde_solver (A,B,C,D,DD,E,DE,F,G,H, &
124  ! YM,YY,RUB1,RUB2,X,NRHO&
125  ! U, V, W)
126 
127  y = y*(1.e0_r8-amix) + yy*amix
128 
129  CALL deriv_holob(nrho,rho,y,dy)
130 
131  solver%Y(idim,:) = y
132  solver%DY(idim,:) = dy
133 
134 20 CONTINUE
135 
136  END DO equation_loop1
137 
138 !#endif
139 
140  RETURN
141 
142 END SUBROUTINE solutionfem
143 
144 
145 !#ifdef WANTCOS
146 
147 ! + + + + + + + + + NUMERICAL PART + + + + + + + + + + +
148 
149 !********************************************************************************
150 !*********************************************************************************
151 SUBROUTINE fem_pde_solver (A,B,C,D,DD,E,DE,F,G,dt_in,&
152  ym,yy,x,n,&
153  u, v, w)
154 !SUBROUTINE fem_pde_solver (A,B,C,D,DD,E,DE,F,G,dt_in,&
155  !YM,YY,RUB1,RUB2,X,N&
156  !U, V, W)
157  USE itm_types
158  IMPLICIT NONE
159 
160  INTEGER, INTENT(in) :: n
161 
162  REAL (R8), INTENT(in) :: x(n), dt_in
163  REAL (R8), INTENT(in) :: a(n),b(n),c(n),d(n),e(n),f(n)
164  REAL (R8), INTENT(in) :: g(n),dd(n),de(n),ym(n)
165 
166  REAL (R8), INTENT(in) :: u(2), v(2), w(2)
167 
168  REAL (R8), INTENT(inOUT) :: yy(n)
169  !
170 
171  !Global matrices
172  REAL (R8), DIMENSION(:,:), ALLOCATABLE :: mat,det,kof,det1,det2,det3
173 
174  REAL (R8), DIMENSION(:,:), ALLOCATABLE :: rhs,rhs_dirichle,&
175  lhs,lhs_dirichle
176 
177  ! pivot and ok are arguments for SGESV()
178  INTEGER, DIMENSION(:), ALLOCATABLE :: pivot, pivot_dirichle
179  INTEGER :: ok
180 
181  !shape functions and their derivatives, gauss points and weights
182  REAL (R8) :: ksi(4), weight(4)
183  REAL (R8) :: n1(4), n2(4)
184  REAL (R8) :: nn1(2,2), nn2(2,2), nn3(2,2), nn4(2,2)
185  REAL (R8) :: dnn1(2,2), dnn2(2,2), dnn3(2,2), dnn4(2,2)
186  REAL (R8) :: dndn(2,2)
187 
188  !counters in loops, auxiliary index
189  INTEGER :: k, i,j, i1, i2
190 
191  !variables inside loop
192  REAL (R8) :: delx
193  REAL (R8) :: a_in(4), b_in(4), c_in(4), d_in(4)
194  REAL (R8) :: e_in(4), f_in(4), g_in(4)
195  REAL (R8) :: de_in(4), dd_in(4)
196  REAL (R8) :: af(4), bf(4), cf(4), df(4)
197 
198  !local matrices
199  REAL (R8) :: matlok(2,2), koflok(1,2)
200  REAL (R8) :: det1lok(2,2),det2lok(2,2), det3lok(2,2)
201 
202  !auxiliary
203  REAL (R8) :: pom1(1,2), pom2(2,1)
204  REAL (R8) :: tmp1(4,2), tmp22(4), tmp2(1,4)
205 
206 
207  !allocation
208  ALLOCATE(mat(n,n))
209  ALLOCATE(det(n,n))
210  ALLOCATE(det1(n,n))
211  ALLOCATE(det2(n,n))
212  ALLOCATE(det3(n,n))
213  ALLOCATE(kof(1,n))
214 
215  ALLOCATE(rhs(1,n))
216  ALLOCATE(lhs(n,n))
217  ALLOCATE(rhs_dirichle(1,n-1))
218  ALLOCATE(lhs_dirichle(n-1,n-1))
219 
220  ALLOCATE(pivot_dirichle(n-1))
221  ALLOCATE(pivot(n))
222 
223  !initialization
224  mat(:,:) = 0.0e0_r8
225  det(:,:) = 0.0e0_r8
226  det1(:,:)= 0.0e0_r8
227  det2(:,:)= 0.0e0_r8
228  det3(:,:)= 0.0e0_r8
229  kof(1,:) = 0.0e0_r8
230 
231  rhs(1,:) =0.0e0_r8
232  lhs(:,:) =0.0e0_r8
233 
234  rhs_dirichle(1,:)=0.0e0_r8
235  lhs_dirichle(:,:)=0.0e0_r8
236 
237 
238  ksi = (/-0.86113631e0_r8,-0.33998104e0_r8,&
239  0.33998104e0_r8,0.86113631e0_r8/)
240  weight = (/ 0.34785485e0_r8,0.65214515e0_r8,&
241  0.65214515e0_r8,0.34785485e0_r8/)
242 
243  n1 = (1-ksi)/2
244  n2 = (1+ksi)/2
245 
246  ! NN1, NN2, NN3, NN4
247  pom1(1,1) = n1(1)
248  pom1(1,2) = n2(1)
249  pom2(1,1) = n1(1)
250  pom2(2,1) = n2(1)
251  nn1 = matmul(pom2,pom1)
252 
253  pom1(1,1) = n1(2)
254  pom1(1,2) = n2(2)
255  pom2(1,1) = n1(2)
256  pom2(2,1) = n2(2)
257  nn2 = matmul(pom2,pom1)
258 
259  pom1(1,1) = n1(3)
260  pom1(1,2) = n2(3)
261  pom2(1,1) = n1(3)
262  pom2(2,1) = n2(3)
263  nn3 = matmul(pom2,pom1)
264 
265  pom1(1,1) = n1(4)
266  pom1(1,2) = n2(4)
267  pom2(1,1) = n1(4)
268  pom2(2,1) = n2(4)
269  nn4 = matmul(pom2,pom1)
270 
271  dndn = reshape((/1, -1, -1, 1/), shape(dndn))
272 
273 
274 
275  solveeqloop: DO k =1, n-1
276 
277  delx = x(k+1) - x(k)
278 
279  ! dNN1, dNN2, dNN3, dNN4
280  pom2(1,1) = -1/delx
281  pom2(2,1) = 1/delx
282 
283  pom1(1,1) = n1(1)
284  pom1(1,2) = n2(1)
285  dnn1 = matmul(pom2,pom1)
286 
287  pom1(1,1) = n1(2)
288  pom1(1,2) = n2(2)
289  dnn2 = matmul(pom2,pom1)
290 
291  pom1(1,1) = n1(3)
292  pom1(1,2) = n2(3)
293  dnn3 = matmul(pom2,pom1)
294 
295  pom1(1,1) = n1(4)
296  pom1(1,2) = n2(4)
297  dnn4 = matmul(pom2,pom1)
298 
299 
300  a_in = a(k) * n1 + a(k+1) *n2
301  b_in = b(k) * n1 + b(k+1) *n2
302  c_in = c(k) * n1 + c(k+1) *n2
303  d_in = d(k) * n1 + d(k+1) *n2
304  e_in = e(k) * n1 + e(k+1) *n2
305  f_in = f(k) * n1 + f(k+1) *n2
306  g_in = g(k) * n1 + g(k+1) *n2
307  de_in = de(k) * n1 + de(k+1) *n2
308  dd_in = dd(k) * n1 + dd(k+1) *n2
309 
310 
311  matlok = (delx/2) *&
312  (( c_in(1)*a_in(1) + c_in(1)*g_in(1)*dt_in + de_in(1)*dt_in )*nn1*weight(1) +&
313  ( c_in(2)*a_in(2) + c_in(2)*g_in(2)*dt_in + de_in(2)*dt_in )*nn2*weight(2) +&
314  ( c_in(3)*a_in(3) + c_in(3)*g_in(3)*dt_in + de_in(3)*dt_in )*nn3*weight(3) +&
315  ( c_in(4)*a_in(4) + c_in(4)*g_in(4)*dt_in + de_in(4)*dt_in )*nn4*weight(4) )
316 
317  det1lok = (delx/2) * (dt_in*e_in(1)*dnn1*weight(1) + dt_in*e_in(2)*dnn2*weight(2) +&
318  dt_in*e_in(3)*dnn3*weight(3) + dt_in*e_in(4)*dnn4*weight(4))
319 
320 
321  det2lok = 1/(2*delx) *(dt_in*d_in(1)*weight(1) + dt_in*d_in(2)*weight(2) +&
322  dt_in*d_in(3)*weight(3) + dt_in*d_in(4)*weight(4)) * dndn
323 
324  det3lok = (delx/2) * (c_in(1)*b_in(1)*nn1*weight(1) + c_in(2)*b_in(2)*nn2*weight(2) +&
325  c_in(3)*b_in(3)*nn3*weight(3) + c_in(4)*b_in(4)*nn4*weight(4))
326 
327 
328 
329  tmp1 = reshape((/n1,n2/),shape(tmp1))
330  tmp22 = (c_in*f_in)*weight*dt_in
331  tmp2 = reshape((/tmp22/),shape(tmp2))
332  koflok = matmul(tmp2,tmp1)
333  koflok = (delx/2)*koflok
334 
335  i1 = k
336  i2 = k+1
337 
338  mat(i1:i2,i1:i2)= mat(i1:i2,i1:i2)+matlok
339  det1(i1:i2,i1:i2)=det1(i1:i2,i1:i2)+det1lok
340  det2(i1:i2,i1:i2)=det2(i1:i2,i1:i2)+det2lok
341  det3(i1:i2,i1:i2)=det3(i1:i2,i1:i2)+det3lok
342  kof(:,i1:i2)=kof(:,i1:i2)+koflok
343 
344  END DO solveeqloop
345 
346  !left-hand side
347  lhs = mat + det1 + det2
348 
349  !right-hand side
350  rhs = reshape( (/ym/), shape(rhs))
351  rhs = matmul(rhs,det3) + kof
352 
353 
354 
355 ! General form of boundary condition is split in two parts:
356 
357 ! Robin's BC, (special case would be Neuman)
358 IF (v(2).NE.0) THEN
359 
360  rhs(1,n)=rhs(1,n)+w(2)/v(2)*d(n)*dt_in ! On RHS w(2)/u(2) should be added (for u(2) = 0 this is basically Neumman's boundary condition
361  lhs(n,n) = lhs(n,n) + u(2)/v(2)*d(n)*dt_in ! On the LHS u(2)/v(2) should be added for the last element in matrix --- this multiplies the psi(N)
362 
363  lhs = transpose(lhs)
364  CALL dgesv(n, 1, lhs, n, pivot, rhs, n, ok)
365  ! the solution is in RHS
366 
367 ! Dirichle's boundary condition
368 ELSEIF (v(2).EQ.0) THEN
369 
370  ! Dirichle's boundary condition requires the following procedure:
371 
372  ! last row and last column thrown out
373  !MAT_Dirichle(:,:) = LHS(1:N-1,1:N-1)
374  do i=1,n-1
375  do j=1,n-1
376  lhs_dirichle(i,j) = lhs(i,j)
377  enddo
378  enddo
379 
380  ! adding the last column multipied with w(2)/u(2) on the RHS
381  !KOF_dirichle(1,:) = RHS(1,1:N-1) - W(2)/U(2) * LHS(1:N-1,N) !LHS(N-1,1:N)
382  do i=1,n-1
383  rhs_dirichle(1,i) = rhs(1,i) - w(2)/u(2) * lhs(n,i)
384  enddo
385  !
386  lhs_dirichle = transpose(lhs_dirichle)
387 
388  call dgesv(n-1, 1, lhs_dirichle, n-1, pivot_dirichle, rhs_dirichle, n-1, ok) !PSI=(M-dt*D_D)\(M*psi.'+dt*K);
389 
390  rhs(1,:n-1) = rhs_dirichle(1,:)
391  rhs(1,n) = w(2)/u(2)
392 
393 ENDIF
394 
395  !return solution
396  rho_loop2: DO i=1, n
397  yy(i) = rhs(1,i)
398  END DO rho_loop2
399 
400  !deallocation
401  IF(ALLOCATED(pivot))DEALLOCATE(pivot)
402  IF(ALLOCATED(pivot_dirichle)) DEALLOCATE(pivot_dirichle)
403 
404  IF(ALLOCATED(mat)) DEALLOCATE(mat)
405  IF(ALLOCATED(det)) DEALLOCATE(det)
406  IF(ALLOCATED(kof)) DEALLOCATE(kof)
407  !DEALLOCATE(MAT,DET,KOF)
408 
409  IF(ALLOCATED(det1)) DEALLOCATE(det1)
410  IF(ALLOCATED(det2)) DEALLOCATE(det2)
411  IF(ALLOCATED(det3)) DEALLOCATE(det3)
412  !DEALLOCATE(DET1,DET2,DET3)
413 
414  IF(ALLOCATED(rhs)) DEALLOCATE(rhs)
415  IF(ALLOCATED(lhs)) DEALLOCATE(lhs)
416  !DEALLOCATE(RHS,LHS)
417 
418  IF(ALLOCATED(rhs_dirichle)) DEALLOCATE(rhs_dirichle)
419  IF(ALLOCATED(lhs_dirichle)) DEALLOCATE(lhs_dirichle)
420 
421 END SUBROUTINE fem_pde_solver
422 
423 
424 
425 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
432 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
433 SUBROUTINE deriv_holob(N,X1,F,FC)
434 
435  USE itm_types
436 
437  IMPLICIT NONE
438 
439  INTEGER, INTENT(in) :: n
440  REAL (r8), INTENT(in) :: x1(n), f(n)
441  REAL (r8), INTENT (out) :: fc(n)
442 
443  fc(1)= (-1.0)*(-1.0/2.0)*((f(1)-f(2))/(x1(1)-x1(2)))&
444  &+(-2.0)*(-1.0) *((f(1)-f(3))/(x1(1)-x1(3)))&
445  &+(-3.0)*(1.0/2.0) *((f(1)-f(4))/(x1(1)-x1(4)))
446 
447  fc(2)= (-1.0)*(-1.0/2.0)*((f(2)-f(3))/(x1(2)-x1(3)))&
448  &+(-2.0)*(-1.0) *((f(2)-f(4))/(x1(2)-x1(4)))&
449  &+(-3.0)*(1.0/2.0) *((f(2)-f(5))/(x1(2)-x1(5)))
450 
451  fc(3)= ((1.0/4.0)*((f(4)-f(2))/(x1(4)-x1(2)))*(2.0*1.0)) &
452  &+((1.0/8.0)*((f(5)-f(1))/(x1(5)-x1(1)))*(2.0*2.0))
453 
454  fc(4:n-1-2)= ((5.0/32.0)*((f(5:n-1-1)-f(3:n-1-3))/(x1(5:n-1-1)-x1(3:n-1-3)))*2.0*1.0)&
455  &+(1.0/8.00)*((f(6:n-1) -f(2:n-1-4))/(x1(6:n-1) -x1(2:n-1-4)))*2.0*2.0&
456  &+(1.0/32.0)*((f(7:n-1+1)-f(1:n-1-5))/(x1(7:n-1+1)-x1(1:n-1-5)))*2.0*3.0
457 
458  fc(n-1-1)= (((1.0/4.0)*((f(n-1+1-1)-f(n-1+1-2))/(x1(n-1+1-1)-x1(n-1+1-2)))*(2*1))&
459  &+((1.0/8.0)*((f(n-1+1) -f(n-1+1-4))/(x1(n-1+1) -x1(n-1+1-4)))*(2*2)))
460 
461  fc(n-1)= (-1.0)*(-1.0/2.0)*((f(n-1+1-1)-f(n-1+1-2))/(x1(n-1+1-1)-x1(n-1+1-2)))&
462  &+(-2.0)*(-1.0) *((f(n-1+1-1)-f(n-1+1-3))/(x1(n-1+1-1)-x1(n-1+1-3)))&
463  &+(-3.0)*(1.0/2.0) *((f(n-1+1-1)-f(n-1+1-4))/(x1(n-1+1-1)-x1(n-1+1-4)))
464 
465  fc(n-1+1)= (-1.0)*(-1.0/2.0)*((f(n-1+1)-f(n-1+1-1))/(x1(n-1+1)-x1(n-1+1-1)))&
466  &+(-2.0)*(-1.0) *((f(n-1+1)-f(n-1+1-2))/(x1(n-1+1)-x1(n-1+1-2)))&
467  &+(-3.0)*(1.0/2.0) *((f(n-1+1)-f(n-1+1-3))/(x1(n-1+1)-x1(n-1+1-3)))
468 
469  RETURN
470 END SUBROUTINE deriv_holob
471 
472 !#endif
solver
Definition: solver.f90:1
type_solver
Definition: type_solver.f90:1
solutionfem
subroutine solutionfem(SOLVER, ifail)
This subroutine is prepared to solve single transport equation in standardised form adopted by the ET...
Definition: SOLVERFEM.f90:11
fem_pde_solver
subroutine fem_pde_solver(A, B, C, D, DD, E, DE, F, G, dt_in, YM, YY, X, N, U, V, W)
Definition: SOLVERFEM.f90:151
deriv_holob
subroutine deriv_holob(N, X1, F, FC)
These subroutines calculate derivative, DY of function Y respect to argument X.
Definition: SOLVERFEM.f90:433
type_solver::numerics
Definition: type_solver.f90:22