SFEMaNS  version 4.1 (work in progress)
Reference documentation for SFEMaNS
 All Classes Files Functions Variables Groups Pages
sub_mass.f90
Go to the documentation of this file.
2  USE my_util
3  USE boundary
4 
6  PRIVATE
7 CONTAINS
8 
9  SUBROUTINE three_level_mass(comm_one_d, time, level_set_LA_P1, level_set_LA_P2, list_mode, &
10  mesh_p1, mesh_p2, chmp_vit_p2, max_vel, level_set_per, density_m2, density_m1, density, &
11  level_set_m1, level_set, visc_entro_level)
12  USE def_type_mesh
13  USE input_data
15  USE my_util
16  IMPLICIT NONE
17  REAL(KIND=8) :: time
18  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
19  type(petsc_csr_la) :: level_set_la_p1
20  type(petsc_csr_la) :: level_set_la_p2
21  TYPE(periodic_type), INTENT(IN) :: level_set_per
22  TYPE(mesh_type), INTENT(IN) :: mesh_p1
23  TYPE(mesh_type), INTENT(IN) :: mesh_p2
24  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: chmp_vit_p2
25  REAL(KIND=8), DIMENSION(:,:,:), INTENT(INOUT) :: density, density_m1, density_m2
26  REAL(KIND=8), DIMENSION(:,:,:,:), INTENT(INOUT) :: level_set, level_set_m1
27  REAL(KIND=8), INTENT(INOUT) :: max_vel
28  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: visc_entro_level
29  REAL(KIND=8), DIMENSION(mesh_P1%np,6,SIZE(list_mode)) :: chmp_vit_p1
30  INTEGER :: n, i, k
31 #include "petsc/finclude/petsc.h"
32  mpi_comm, DIMENSION(:), POINTER :: comm_one_d
33 
34  IF (inputs%if_level_set) THEN
35 
36  IF (inputs%if_level_set_P2) THEN
37  DO n = 1, inputs%nb_fluid-1
38  CALL three_level_level_set(comm_one_d, time, level_set_la_p2, inputs%dt, list_mode, &
39  mesh_p2, level_set_m1(n,:,:,:), level_set(n,:,:,:), chmp_vit_p2, max_vel, &
40  inputs%my_par_level_set, inputs%level_set_list_dirichlet_sides, level_set_per, n, &
41  visc_entro_level)
42  END DO
43  ELSE
44  DO i = 1, SIZE(list_mode)
45  DO k = 1, 6
46  CALL project_p2_p1(mesh_p2%jj, mesh_p1%jj, chmp_vit_p2(:,k,i), chmp_vit_p1(:,k,i))
47  END DO
48  END DO
49 
50  DO n = 1, inputs%nb_fluid-1
51  CALL three_level_level_set(comm_one_d, time, level_set_la_p1, inputs%dt, list_mode, &
52  mesh_p1, level_set_m1(n,:,:,:), level_set(n,:,:,:), chmp_vit_p1, max_vel, &
53  inputs%my_par_level_set, inputs%level_set_list_dirichlet_sides, level_set_per, n, &
54  visc_entro_level)
55  END DO
56  END IF
57 
58  !===Update densities
59  density_m2 = density_m1
60  density_m1 = density
61  CALL reconstruct_variable(comm_one_d, list_mode, mesh_p1, mesh_p2, level_set, &
62  inputs%density_fluid, density)
63  RETURN
64  ELSE
65  RETURN
66  END IF
67 
68  END SUBROUTINE three_level_mass
69 
70  SUBROUTINE reconstruct_variable(comm_one_d, list_mode, mesh_P1, mesh_P2, level_set, values, variable)
71  !==============================
72  USE def_type_mesh
73  USE sft_parallele
74  USE input_data
75  IMPLICIT NONE
76  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
77  TYPE(mesh_type), INTENT(IN) :: mesh_p1
78  TYPE(mesh_type), INTENT(IN) :: mesh_p2
79  REAL(KIND=8), DIMENSION(:,:,:,:), INTENT(IN) :: level_set
80  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: values
81  REAL(KIND=8), DIMENSION(:,:,:), INTENT(INOUT) :: variable
82  LOGICAL, SAVE :: once = .true.
83  INTEGER, SAVE :: m_max_c
84  INTEGER, SAVE :: m_max_pad
85  INTEGER, SAVE :: bloc_size
86  INTEGER, SAVE :: nb_procs
87  INTEGER :: i, code, k, nb_inter
88  REAL(KIND=8), DIMENSION(mesh_P2%np,2,SIZE(list_mode)) :: rho_phi
89  REAL(KIND=8), DIMENSION(inputs%nb_fluid-1,mesh_P2%np,2,SIZE(list_mode)) :: level_set_p2
90  !Communicators for Petsc, in space and Fourier------------------------------
91 #include "petsc/finclude/petsc.h"
92  mpi_comm, DIMENSION(:), POINTER :: comm_one_d
93  !------------------------------END OF DECLARATION--------------------------------------
94 
95  IF (once) THEN
96  once = .false.
97  !-------------DIMENSIONS-------------------------------------------------------
98  m_max_c = SIZE(list_mode)
99  !-------------FFT VARIABLES----------------------------------------------------
100  CALL mpi_comm_size(comm_one_d(2), nb_procs, code)
101  bloc_size = mesh_p2%np/nb_procs+1
102  m_max_pad = 3*m_max_c*nb_procs/2
103  END IF
104 
105  IF (.NOT.inputs%if_level_set) THEN
106  ! we don't reconstruct density, viscosity if no level set
107  RETURN
108  ELSE
109  IF (inputs%if_level_set_P2) THEN
110  level_set_p2=level_set
111  ELSE
112  DO nb_inter = 1, inputs%nb_fluid-1
113  DO i = 1, SIZE(list_mode)
114  DO k = 1, 2
115  CALL inject_p1_p2(mesh_p1%jj, mesh_p2%jj, level_set(nb_inter,:,k,i), level_set_p2(nb_inter,:,k,i))
116  END DO
117  END DO
118  END DO
119  END IF
120  IF (maxval(abs(values(1)-values(:))) .LE. 1.d-6) THEN
121  variable = 0.d0
122  DO i = 1, m_max_c
123  IF (list_mode(i)==0) THEN
124  variable(:,1,i) = values(1)
125  END IF
126  END DO
127  ELSE IF (inputs%level_set_reconstruction_type == 'lin') THEN
128  IF (inputs%if_kill_overshoot) THEN
129  IF (nb_procs==1.AND.SIZE(list_mode)==1.AND.list_mode(1)==0) THEN !case axisym
130  level_set_p2 = min(1.d0, level_set_p2)
131  level_set_p2 = max(0.d0, level_set_p2)
132  ELSE !level set depends of theta
133  DO k = 1, inputs%nb_fluid-1
134  CALL fft_no_overshoot_level_set(comm_one_d(2), level_set_p2(k,:,:,:), &
135  nb_procs, bloc_size, m_max_pad)
136  END DO
137  END IF
138  END IF
139  variable = 0.d0
140  DO i = 1, m_max_c
141  IF (list_mode(i)==0) THEN
142  variable(:,1,i) = values(1)
143  END IF
144  END DO
145  variable = variable + (values(2)-values(1))*level_set_p2(1,:,:,:)
146  IF (inputs%nb_fluid.GE.3) THEN
147  DO i = 1, inputs%nb_fluid-2
148  CALL fft_par_prod_dcl(comm_one_d(2), variable, level_set_p2(i+1,:,:,:), rho_phi, &
149  nb_procs, bloc_size, m_max_pad)
150  variable = variable -rho_phi + values(i+2)*level_set_p2(i+1,:,:,:)
151  END DO
152  END IF
153  ELSE
154  CALL fft_heaviside_dcl(comm_one_d(2), level_set_p2, values, variable, &
155  nb_procs, bloc_size, m_max_pad, inputs%level_set_tanh_coeff_reconstruction)
156  DO i = 1, m_max_c
157  IF (list_mode(i)==0) THEN
158  variable(:,2,i) = 0.d0
159  END IF
160  END DO
161  END IF
162  END IF
163  END SUBROUTINE reconstruct_variable
164 
165  SUBROUTINE total_mass(comm_one_d, list_mode, mass_mesh, level_set, mass_tot)
166  !==============================
167  USE def_type_mesh
168  USE sft_parallele
169  USE input_data
170  IMPLICIT NONE
171  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
172  TYPE(mesh_type), INTENT(IN) :: mass_mesh
173  REAL(KIND=8), DIMENSION(:,:,:,:), INTENT(IN) :: level_set
174  REAL(KIND=8), INTENT(OUT) :: mass_tot
175  REAL(KIND=8), DIMENSION(SIZE(level_set,2),SIZE(level_set,3), & SIZE(level_set,4)) :: density_loc
176  INTEGER :: m_max_pad, bloc_size, nb_procs
177  INTEGER :: i, code, my_petscworld_rank, m, l
178  REAL(KIND=8) :: mass_loc, mass_f, ray
179  REAL(KIND=8), DIMENSION(mass_mesh%np,2,SIZE(list_mode)) :: rho_phi
180  INTEGER, DIMENSION(mass_mesh%gauss%n_w) :: j_loc
181  REAL(KIND=8) :: pi= 3.14159265358979323846d0
182  !Communicators for Petsc, in space and Fourier------------------------------
183 #include "petsc/finclude/petsc.h"
184  mpi_comm, DIMENSION(:), POINTER :: comm_one_d
185  !------------------------------END OF DECLARATION--------------------------------------
186 
187  CALL mpi_comm_rank(petsc_comm_world,my_petscworld_rank,code)
188  CALL mpi_comm_size(comm_one_d(2), nb_procs, code)
189 
190  IF(.NOT.inputs%if_level_set) THEN
191  CALL error_petsc('BUG in sub_mass : you should not compute any mass')
192  ELSE
193  IF (inputs%level_set_reconstruction_type == 'lin') THEN
194  density_loc = 0.d0
195  DO i = 1, SIZE(list_mode)
196  IF (list_mode(i)==0) THEN
197  density_loc(:,1,i) = inputs%density_fluid(1)
198  END IF
199  END DO
200  density_loc = density_loc + (inputs%density_fluid(2)-inputs%density_fluid(1))*level_set(i,:,:,:)
201 
202  bloc_size = SIZE(level_set,2)/nb_procs+1
203  m_max_pad = 3*SIZE(list_mode)*nb_procs/2
204  IF (inputs%nb_fluid.GE.3) THEN
205  DO i = 1, inputs%nb_fluid-2
206  CALL fft_par_prod_dcl(comm_one_d(2), density_loc, level_set(i+1,:,:,:), rho_phi, &
207  nb_procs, bloc_size, m_max_pad)
208  density_loc = density_loc -rho_phi + inputs%density_fluid(i+2)*level_set(i+1,:,:,:)
209  END DO
210  END IF
211  ELSE
212  bloc_size = SIZE(level_set,2)/nb_procs+1
213  m_max_pad = 3*SIZE(list_mode)*nb_procs/2
214  CALL fft_heaviside_dcl(comm_one_d(2), level_set, inputs%density_fluid, &
215  density_loc, nb_procs, bloc_size, m_max_pad, inputs%level_set_tanh_coeff_reconstruction)
216  END IF
217 
218  mass_loc = 0.d0
219  DO i = 1, SIZE(list_mode)
220  IF (list_mode(i)==0) THEN
221  DO m = 1, mass_mesh%me
222  j_loc = mass_mesh%jj(:,m)
223  DO l = 1, mass_mesh%gauss%l_G
224  !===Compute radius of Gauss point
225  ray = sum(mass_mesh%rr(1,j_loc)*mass_mesh%gauss%ww(:,l))
226  mass_loc = mass_loc + sum(density_loc(j_loc,1,i)* &
227  mass_mesh%gauss%ww(:,l))*ray*mass_mesh%gauss%rj(l,m)
228  END DO
229  END DO
230  END IF
231  END DO
232  mass_loc = mass_loc*2*pi
233  CALL mpi_allreduce(mass_loc, mass_f, 1, mpi_double_precision, mpi_sum, &
234  comm_one_d(2), code)
235  CALL mpi_allreduce(mass_f, mass_tot, 1, mpi_double_precision, mpi_sum, &
236  comm_one_d(1), code)
237  END IF
238 
239  END SUBROUTINE total_mass
240 
241  SUBROUTINE inject_p1_p2(jj_c, jj_f, pp_c, pp_f)
243  IMPLICIT NONE
244  INTEGER, DIMENSION(:,:), INTENT(IN) :: jj_c, jj_f
245  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: pp_c
246  REAL(KIND=8), DIMENSION(:), INTENT(OUT) :: pp_f
247  REAL(KIND=8) :: half = 0.5
248  INTEGER:: m
249  IF (SIZE(jj_c,1)==3) THEN
250  DO m = 1, SIZE(jj_f,2)
251  pp_f(jj_f(1:3,m)) = pp_c(jj_c(:,m))
252  pp_f(jj_f(4,m)) = (pp_c(jj_c(2,m)) + pp_c(jj_c(3,m)))*half
253  pp_f(jj_f(5,m)) = (pp_c(jj_c(3,m)) + pp_c(jj_c(1,m)))*half
254  pp_f(jj_f(6,m)) = (pp_c(jj_c(1,m)) + pp_c(jj_c(2,m)))*half
255  END DO
256  ELSE
257  CALL error_petsc('BUG in inject_P1_P2: finite element not yet programmed')
258  END IF
259 
260  END SUBROUTINE inject_p1_p2
261 
262  SUBROUTINE project_p2_p1(jj_P2, jj_P1, pp_P2, pp_P1)
264  IMPLICIT NONE
265  INTEGER, DIMENSION(:,:), INTENT(IN) :: jj_p2, jj_p1
266  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: pp_p2
267  REAL(KIND=8), DIMENSION(:), INTENT(OUT) :: pp_p1
268  INTEGER:: m
269 
270  IF (SIZE(jj_p1,1)==3) THEN
271  DO m = 1, SIZE(jj_p1,2)
272  pp_p1(jj_p1(:,m)) = pp_p2(jj_p2(1:3,m))
273  END DO
274  ELSE
275  CALL error_petsc('BUG in inject_P2_P1: finite element not yet programmed')
276  END IF
277 
278  END SUBROUTINE project_p2_p1
279 
280 END MODULE subroutine_mass
281 
subroutine, public reconstruct_variable(comm_one_d, list_mode, mesh_P1, mesh_P2, level_set, values, variable)
Definition: sub_mass.f90:70
subroutine, public total_mass(comm_one_d, list_mode, mass_mesh, level_set, mass_tot)
Definition: sub_mass.f90:165
subroutine, public fft_par_prod_dcl(communicator, c1_in, c2_in, c_out, nb_procs, bloc_size, m_max_pad, temps)
subroutine, public three_level_level_set(comm_one_d, time, cc_1_LA, dt, list_mode, cc_mesh, cn_m1, cn, chmp_vit, max_vel, my_par_cc, cc_list_dirichlet_sides, cc_per, nb_inter, visc_entro_level)
subroutine, public inject_p1_p2(jj_c, jj_f, pp_c, pp_f)
Definition: sub_mass.f90:242
subroutine, public fft_no_overshoot_level_set(communicator, c1_inout, nb_procs, bloc_size, m_max_pad, temps)
subroutine, public fft_heaviside_dcl(communicator, V1_in, values, V_out, nb_procs, bloc_size, m_max_pad, coeff_tanh, temps)
subroutine mass_tot(communicator, mesh, tempn, RESLT)
subroutine, public project_p2_p1(jj_P2, jj_P1, pp_P2, pp_P1)
Definition: sub_mass.f90:263
subroutine error_petsc(string)
Definition: my_util.f90:15
subroutine, public three_level_mass(comm_one_d, time, level_set_LA_P1, level_set_LA_P2, list_mode, mesh_P1, mesh_P2, chmp_vit_P2, max_vel, level_set_per, density_m2, density_m1, density, level_set_m1, level_set, visc_entro_level)
Definition: sub_mass.f90:9