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