ETS-Core  version:0.0.4-46-ge2d8
Core actors for the ETS-6
 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  USE ids_schemas
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 (IDS_REAL) :: rho(solver%nrho) !radii
44  !REAL (IDS_REAL) :: RHOMAX
45 
46  REAL (IDS_REAL) :: amix !fraction of new sollution mixed
47 
48  REAL (IDS_REAL) :: y(solver%nrho), ym(solver%nrho) !function at the current amd previous time steps
49  REAL (IDS_REAL) :: dy(solver%nrho) !derivative of function
50 
51  REAL (IDS_REAL) :: yy(solver%nrho) !new function
52 
53  REAL (IDS_REAL) :: a(solver%nrho), b(solver%nrho) !coefficients
54  REAL (IDS_REAL) :: c(solver%nrho), d(solver%nrho) !coefficients
55  REAL (IDS_REAL) :: e(solver%nrho), f(solver%nrho) !coefficients
56  REAL (IDS_REAL) :: g(solver%nrho), h !coefficients
57  REAL (IDS_REAL) :: dd(solver%nrho), de(solver%nrho) !derivatives of coefficients
58 
59  REAL (IDS_REAL) :: v(2), u(2), w(2) !boundary conditions
60 
61 ! +++ Coefficients used by internal numerical solver:
62  !REAL (IDS_REAL) :: X(SOLVER%NRHO) !normalized radii
63  !REAL (IDS_REAL) :: 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 
108  CALL deriv_holob(nrho,rho,d,dd)
109  CALL deriv_holob(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_solver(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_ids_real-amix) + yy*amix
129 
130  CALL deriv_holob(nrho,rho,y,dy)
131 
132  solver%Y(idim,1:solver%NRHO) = y
133  solver%DY(idim,1:solver%NRHO) = dy
134 
135 20 CONTINUE
136 
137  END DO equation_loop1
138 
139 !#endif
140 
141  RETURN
142 
143 END SUBROUTINE solutionfem
144 
145 
146 !#ifdef WANTCOS
147 
148 ! + + + + + + + + + NUMERICAL PART + + + + + + + + + + +
149 
150 !********************************************************************************
151 !*********************************************************************************
152 SUBROUTINE fem_pde_solver (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
159  IMPLICIT NONE
160 
161  INTEGER, INTENT(in) :: n
162 
163  REAL (IDS_REAL), INTENT(in) :: x(n), dt_in
164  REAL (IDS_REAL), INTENT(in) :: a(n),b(n),c(n),d(n),e(n),f(n)
165  REAL (IDS_REAL), INTENT(in) :: g(n),dd(n),de(n),ym(n)
166 
167  REAL (IDS_REAL), INTENT(in) :: u(2), v(2), w(2)
168 
169  REAL (IDS_REAL), INTENT(inOUT) :: yy(n)
170  !
171 
172  !Global matrices
173  REAL (IDS_REAL), DIMENSION(:,:), ALLOCATABLE :: mat,det,kof,det1,det2,det3
174 
175  REAL (IDS_REAL), 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 (IDS_REAL) :: ksi(4), weight(4)
184  REAL (IDS_REAL) :: n1(4), n2(4)
185  REAL (IDS_REAL) :: nn1(2,2), nn2(2,2), nn3(2,2), nn4(2,2)
186  REAL (IDS_REAL) :: dnn1(2,2), dnn2(2,2), dnn3(2,2), dnn4(2,2)
187  REAL (IDS_REAL) :: dndn(2,2)
188 
189  !counters in loops, auxiliary index
190  INTEGER :: k, i,j, i1, i2
191 
192  !variables inside loop
193  REAL (IDS_REAL) :: delx
194  REAL (IDS_REAL) :: a_in(4), b_in(4), c_in(4), d_in(4)
195  REAL (IDS_REAL) :: e_in(4), f_in(4), g_in(4)
196  REAL (IDS_REAL) :: de_in(4), dd_in(4)
197  REAL (IDS_REAL) :: af(4), bf(4), cf(4), df(4)
198 
199  !local matrices
200  REAL (IDS_REAL) :: matlok(2,2), koflok(1,2)
201  REAL (IDS_REAL) :: det1lok(2,2),det2lok(2,2), det3lok(2,2)
202 
203  !auxiliary
204  REAL (IDS_REAL) :: pom1(1,2), pom2(2,1)
205  REAL (IDS_REAL) :: tmp1(4,2), tmp22(4), tmp2(1,4)
206 
207 
208  !allocation
209  ALLOCATE(mat(n,n))
210  ALLOCATE(det(n,n))
211  ALLOCATE(det1(n,n))
212  ALLOCATE(det2(n,n))
213  ALLOCATE(det3(n,n))
214  ALLOCATE(kof(1,n))
215 
216  ALLOCATE(rhs(1,n))
217  ALLOCATE(lhs(n,n))
218  ALLOCATE(rhs_dirichle(1,n-1))
219  ALLOCATE(lhs_dirichle(n-1,n-1))
220 
221  ALLOCATE(pivot_dirichle(n-1))
222  ALLOCATE(pivot(n))
223 
224  !initialization
225  mat(:,:) = 0.0e0_ids_real
226  det(:,:) = 0.0e0_ids_real
227  det1(:,:)= 0.0e0_ids_real
228  det2(:,:)= 0.0e0_ids_real
229  det3(:,:)= 0.0e0_ids_real
230  kof(1,:) = 0.0e0_ids_real
231 
232  rhs(1,:) =0.0e0_ids_real
233  lhs(:,:) =0.0e0_ids_real
234 
235  rhs_dirichle(1,:)=0.0e0_ids_real
236  lhs_dirichle(:,:)=0.0e0_ids_real
237 
238 
239  ksi = (/-0.86113631e0_ids_real,-0.33998104e0_ids_real,&
240  0.33998104e0_ids_real,0.86113631e0_ids_real/)
241  weight = (/ 0.34785485e0_ids_real,0.65214515e0_ids_real,&
242  0.65214515e0_ids_real,0.34785485e0_ids_real/)
243 
244  n1 = (1-ksi)/2
245  n2 = (1+ksi)/2
246 
247  ! NN1, NN2, NN3, NN4
248  pom1(1,1) = n1(1)
249  pom1(1,2) = n2(1)
250  pom2(1,1) = n1(1)
251  pom2(2,1) = n2(1)
252  nn1 = matmul(pom2,pom1)
253 
254  pom1(1,1) = n1(2)
255  pom1(1,2) = n2(2)
256  pom2(1,1) = n1(2)
257  pom2(2,1) = n2(2)
258  nn2 = matmul(pom2,pom1)
259 
260  pom1(1,1) = n1(3)
261  pom1(1,2) = n2(3)
262  pom2(1,1) = n1(3)
263  pom2(2,1) = n2(3)
264  nn3 = matmul(pom2,pom1)
265 
266  pom1(1,1) = n1(4)
267  pom1(1,2) = n2(4)
268  pom2(1,1) = n1(4)
269  pom2(2,1) = n2(4)
270  nn4 = matmul(pom2,pom1)
271 
272  dndn = reshape((/1, -1, -1, 1/), shape(dndn))
273 
274 
275 
276  solveeqloop: DO k =1, n-1
277 
278  delx = x(k+1) - x(k)
279 
280  ! dNN1, dNN2, dNN3, dNN4
281  pom2(1,1) = -1/delx
282  pom2(2,1) = 1/delx
283 
284  pom1(1,1) = n1(1)
285  pom1(1,2) = n2(1)
286  dnn1 = matmul(pom2,pom1)
287 
288  pom1(1,1) = n1(2)
289  pom1(1,2) = n2(2)
290  dnn2 = matmul(pom2,pom1)
291 
292  pom1(1,1) = n1(3)
293  pom1(1,2) = n2(3)
294  dnn3 = matmul(pom2,pom1)
295 
296  pom1(1,1) = n1(4)
297  pom1(1,2) = n2(4)
298  dnn4 = matmul(pom2,pom1)
299 
300 
301  a_in = a(k) * n1 + a(k+1) *n2
302  b_in = b(k) * n1 + b(k+1) *n2
303  c_in = c(k) * n1 + c(k+1) *n2
304  d_in = d(k) * n1 + d(k+1) *n2
305  e_in = e(k) * n1 + e(k+1) *n2
306  f_in = f(k) * n1 + f(k+1) *n2
307  g_in = g(k) * n1 + g(k+1) *n2
308  de_in = de(k) * n1 + de(k+1) *n2
309  dd_in = dd(k) * n1 + dd(k+1) *n2
310 
311 
312  matlok = (delx/2) *&
313  (( c_in(1)*a_in(1) + c_in(1)*g_in(1)*dt_in + de_in(1)*dt_in )*nn1*weight(1) +&
314  ( c_in(2)*a_in(2) + c_in(2)*g_in(2)*dt_in + de_in(2)*dt_in )*nn2*weight(2) +&
315  ( c_in(3)*a_in(3) + c_in(3)*g_in(3)*dt_in + de_in(3)*dt_in )*nn3*weight(3) +&
316  ( c_in(4)*a_in(4) + c_in(4)*g_in(4)*dt_in + de_in(4)*dt_in )*nn4*weight(4) )
317 
318  det1lok = (delx/2) * (dt_in*e_in(1)*dnn1*weight(1) + dt_in*e_in(2)*dnn2*weight(2) +&
319  dt_in*e_in(3)*dnn3*weight(3) + dt_in*e_in(4)*dnn4*weight(4))
320 
321 
322  det2lok = 1/(2*delx) *(dt_in*d_in(1)*weight(1) + dt_in*d_in(2)*weight(2) +&
323  dt_in*d_in(3)*weight(3) + dt_in*d_in(4)*weight(4)) * dndn
324 
325  det3lok = (delx/2) * (c_in(1)*b_in(1)*nn1*weight(1) + c_in(2)*b_in(2)*nn2*weight(2) +&
326  c_in(3)*b_in(3)*nn3*weight(3) + c_in(4)*b_in(4)*nn4*weight(4))
327 
328 
329 
330  tmp1 = reshape((/n1,n2/),shape(tmp1))
331  tmp22 = (c_in*f_in)*weight*dt_in
332  tmp2 = reshape((/tmp22/),shape(tmp2))
333  koflok = matmul(tmp2,tmp1)
334  koflok = (delx/2)*koflok
335 
336  i1 = k
337  i2 = k+1
338 
339  mat(i1:i2,i1:i2)= mat(i1:i2,i1:i2)+matlok
340  det1(i1:i2,i1:i2)=det1(i1:i2,i1:i2)+det1lok
341  det2(i1:i2,i1:i2)=det2(i1:i2,i1:i2)+det2lok
342  det3(i1:i2,i1:i2)=det3(i1:i2,i1:i2)+det3lok
343  kof(:,i1:i2)=kof(:,i1:i2)+koflok
344 
345  END DO solveeqloop
346 
347  !left-hand side
348  lhs = mat + det1 + det2
349 
350  !right-hand side
351  rhs = reshape( (/ym/), shape(rhs))
352  rhs = matmul(rhs,det3) + kof
353 
354 
355 
356 ! General form of boundary condition is split in two parts:
357 
358 ! Robin's BC, (special case would be Neuman)
359 IF (v(2).NE.0) THEN
360 
361  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
362  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)
363 
364  lhs = transpose(lhs)
365  CALL dgesv(n, 1, lhs, n, pivot, rhs, n, ok)
366  ! the solution is in RHS
367 
368 ! Dirichle's boundary condition
369 ELSEIF (v(2).EQ.0) THEN
370 
371  ! Dirichle's boundary condition requires the following procedure:
372 
373  ! last row and last column thrown out
374  !MAT_Dirichle(:,:) = LHS(1:N-1,1:N-1)
375  do i=1,n-1
376  do j=1,n-1
377  lhs_dirichle(i,j) = lhs(i,j)
378  enddo
379  enddo
380 
381  ! adding the last column multipied with w(2)/u(2) on the RHS
382  !KOF_dirichle(1,:) = RHS(1,1:N-1) - W(2)/U(2) * LHS(1:N-1,N) !LHS(N-1,1:N)
383  do i=1,n-1
384  rhs_dirichle(1,i) = rhs(1,i) - w(2)/u(2) * lhs(n,i)
385  enddo
386  !
387  lhs_dirichle = transpose(lhs_dirichle)
388 
389  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);
390 
391  rhs(1,:n-1) = rhs_dirichle(1,:)
392  rhs(1,n) = w(2)/u(2)
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(det)) DEALLOCATE(det)
407  IF(ALLOCATED(kof)) DEALLOCATE(kof)
408  !DEALLOCATE(MAT,DET,KOF)
409 
410  IF(ALLOCATED(det1)) DEALLOCATE(det1)
411  IF(ALLOCATED(det2)) DEALLOCATE(det2)
412  IF(ALLOCATED(det3)) DEALLOCATE(det3)
413  !DEALLOCATE(DET1,DET2,DET3)
414 
415  IF(ALLOCATED(rhs)) DEALLOCATE(rhs)
416  IF(ALLOCATED(lhs)) DEALLOCATE(lhs)
417  !DEALLOCATE(RHS,LHS)
418 
419  IF(ALLOCATED(rhs_dirichle)) DEALLOCATE(rhs_dirichle)
420  IF(ALLOCATED(lhs_dirichle)) DEALLOCATE(lhs_dirichle)
421 
422 END SUBROUTINE fem_pde_solver
423 
424 
425 
426 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
433 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
434 SUBROUTINE deriv_holob(N,X1,F,FC)
435 
436  USE ids_schemas
437 
438  IMPLICIT NONE
439 
440  INTEGER, INTENT(in) :: n
441  REAL (IDS_REAL), INTENT(in) :: x1(n), f(n)
442  REAL (IDS_REAL), INTENT (out) :: fc(n)
443 
444  fc(1)= (-1.0)*(-1.0/2.0)*((f(1)-f(2))/(x1(1)-x1(2)))&
445  &+(-2.0)*(-1.0) *((f(1)-f(3))/(x1(1)-x1(3)))&
446  &+(-3.0)*(1.0/2.0) *((f(1)-f(4))/(x1(1)-x1(4)))
447 
448  fc(2)= (-1.0)*(-1.0/2.0)*((f(2)-f(3))/(x1(2)-x1(3)))&
449  &+(-2.0)*(-1.0) *((f(2)-f(4))/(x1(2)-x1(4)))&
450  &+(-3.0)*(1.0/2.0) *((f(2)-f(5))/(x1(2)-x1(5)))
451 
452  fc(3)= ((1.0/4.0)*((f(4)-f(2))/(x1(4)-x1(2)))*(2.0*1.0)) &
453  &+((1.0/8.0)*((f(5)-f(1))/(x1(5)-x1(1)))*(2.0*2.0))
454 
455  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)&
456  &+(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&
457  &+(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
458 
459  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))&
460  &+((1.0/8.0)*((f(n-1+1) -f(n-1+1-4))/(x1(n-1+1) -x1(n-1+1-4)))*(2*2)))
461 
462  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)))&
463  &+(-2.0)*(-1.0) *((f(n-1+1-1)-f(n-1+1-3))/(x1(n-1+1-1)-x1(n-1+1-3)))&
464  &+(-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)))
465 
466  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)))&
467  &+(-2.0)*(-1.0) *((f(n-1+1)-f(n-1+1-2))/(x1(n-1+1)-x1(n-1+1-2)))&
468  &+(-3.0)*(1.0/2.0) *((f(n-1+1)-f(n-1+1-3))/(x1(n-1+1)-x1(n-1+1-3)))
469 
470  RETURN
471 END SUBROUTINE deriv_holob
472 
473 !#endif
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
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:152
deriv_holob
subroutine deriv_holob(N, X1, F, FC)
These subroutines calculate derivative, DY of function Y respect to argument X.
Definition: SOLVERFEM.f90:434
type_transport_solver_numerics
The module declares types of variables used by numerical solvers.
Definition: type_transport_solver_numerics.f90:12
Generated by   doxygen 1.8.5