sources_viscosity.f90
Go to the documentation of this file.
1!#############################################################################
2!# #
3!# fosite - 3D hydrodynamical simulation program #
4!# module: sources_viscosity.f90 #
5!# #
6!# Copyright (C) 2008-2024 #
7!# Bjoern Sperling <sperling@astrophysik.uni-kiel.de> #
8!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de> #
9!# #
10!# This program is free software; you can redistribute it and/or modify #
11!# it under the terms of the GNU General Public License as published by #
12!# the Free Software Foundation; either version 2 of the License, or (at #
13!# your option) any later version. #
14!# #
15!# This program is distributed in the hope that it will be useful, but #
16!# WITHOUT ANY WARRANTY; without even the implied warranty of #
17!# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or #
18!# NON INFRINGEMENT. See the GNU General Public License for more #
19!# details. #
20!# #
21!# You should have received a copy of the GNU General Public License #
22!# along with Mesh program; if not, write to the Free Software #
23!# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. #
24!# #
25!#############################################################################
37!----------------------------------------------------------------------------!
60!----------------------------------------------------------------------------!
61MODULE sources_viscosity_mod
62 USE sources_base_mod
63 USE physics_base_mod
64 USE fluxes_base_mod
65 USE mesh_base_mod
66 USE marray_base_mod
67 USE marray_compound_mod
68 USE logging_base_mod
69 USE common_dict
70 IMPLICIT NONE
71 !--------------------------------------------------------------------------!
72 PRIVATE
73 CHARACTER(LEN=32), PARAMETER :: source_name = "viscosity of Newtonian fluid"
74 ! flags for viscosity model
75 INTEGER,PARAMETER :: molecular = 1 ! constant viscosity
76 INTEGER,PARAMETER :: alpha = 2 ! Shakura-Sunyaev prescription
77 INTEGER,PARAMETER :: beta = 3 ! Duschl prescription
78 INTEGER,PARAMETER :: powerlaw = 4 ! power law depending on spec. angular momentum
79 INTEGER,PARAMETER :: alpha_alt = 5 ! alternative Shakura-Sunyaev
80 CHARACTER(LEN=32),PARAMETER, DIMENSION(5) :: viscosity_name = (/ &
81 "constant viscosity ", &
82 "turbulent Shakura-Sunyaev ", &
83 "turbulent Duschl ", &
84 "power law viscosity ", &
85 "alternative Shakura-Sunyaev "/)
86
87 TYPE, EXTENDS(sources_base) :: sources_viscosity
88 CLASS(logging_base), ALLOCATABLE :: viscosity
89 CLASS(marray_base), ALLOCATABLE :: ephir, & !< azimuthal unit vector */ distance to axis
90 ephir_tmp, &
91 dynvis, & !< dynamic viscosity
92 kinvis, & !< kinematic viscosity
93 bulkvis
94 REAL :: dynconst,bulkconst,power
95 REAL, DIMENSION(:,:,:), POINTER :: & !< comp. of stress tensor
96 btxx,btyy,btzz,btxy,btxz,btyz,tmp,tmp2,tmp3
97 REAL, DIMENSION(:,:), POINTER :: &
98 sxx,syy,szz,sxy,sxz,syz
99 CONTAINS
100 PROCEDURE :: initsources
101 PROCEDURE :: infosources
102 PROCEDURE :: setoutput
103 PROCEDURE :: updateviscosity
104 PROCEDURE :: externalsources
105 PROCEDURE :: calctimestep
106 final :: finalize
107 END TYPE
108 !--------------------------------------------------------------------------!
109 PUBLIC :: &
110 ! types
111 sources_viscosity, &
112 ! constants
113 molecular,alpha,beta,powerlaw,alpha_alt,viscosity_name
114 !--------------------------------------------------------------------------!
115
116CONTAINS
117
118 SUBROUTINE initsources(this,Mesh,Physics,Fluxes,config,IO)
119 USE geometry_base_mod
120 USE geometry_cylindrical_mod, ONLY: geometry_cylindrical
121 USE geometry_logcylindrical_mod, ONLY: geometry_logcylindrical
122 USE physics_eulerisotherm_mod, ONLY : physics_eulerisotherm
123 IMPLICIT NONE
124 !------------------------------------------------------------------------!
125 CLASS(sources_viscosity), INTENT(INOUT) :: this
126 CLASS(mesh_base), INTENT(IN) :: Mesh
127 CLASS(physics_base), INTENT(IN) :: Physics
128 CLASS(fluxes_base), INTENT(IN) :: Fluxes
129 TYPE(dict_typ), POINTER :: config, IO
130 !------------------------------------------------------------------------!
131 INTEGER :: stype
132 INTEGER :: err, viscosity_number ,k,l
133 !------------------------------------------------------------------------!
134 CALL getattr(config, "stype", stype)
135 CALL this%InitSources_base(stype,source_name)
136
137 SELECT TYPE(phys => physics)
138 CLASS IS(physics_eulerisotherm)
139 IF (phys%VDIM.EQ.1) &
140 CALL this%Error("InitSources_viscosity","viscosity is currently not supported in 1D physics")
141 CLASS DEFAULT
142 CALL this%Error("InitSources_viscosity", &
143 "viscosity is currently only supported in eulerisotherm/euler physics")
144 END SELECT
145
146 ! viscosity model
147 CALL getattr(config, "vismodel", viscosity_number)
148 ALLOCATE(logging_base::this%viscosity)
149 CALL this%viscosity%InitLogging(viscosity_number,viscosity_name(viscosity_number))
150
151 ! dynamic viscosity constant
152 CALL getattr(config, "dynconst", this%dynconst, 0.1)
153 ! bulk viscosity constant
154 ! At the moment this only affects the molecular viscosity. The default
155 ! value refers to a trace free viscosity tensor.
156 ! All other viscosity models should always be trace free. Therefore this
157 ! is hard coded and cannot be changed.
158 CALL getattr(config, "bulkconst", this%bulkconst, -2./3.*this%dynconst)
159
160 ! power law exponent used for POWERLAW viscosity;
161 ! defaults to constant kinematic viscosity, i.e. nu ~ l^0 = const
162 CALL getattr(config, "exponent", this%power, 0.0)
163
164 ! set viscous courant number
165 CALL getattr(config, "cvis", this%cvis, 0.5)
166
167 ALLOCATE(this%dynvis,this%kinvis,this%bulkvis, &
168 this%btxx(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX), &
169 this%btyy(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX), &
170 this%btzz(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX), &
171 this%btxy(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX), &
172 this%btxz(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX), &
173 this%btyz(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX), &
174 stat=err)
175 IF (err.NE.0) &
176 CALL this%Error("InitSources_viscosity","Memory allocation failed.")
177
178 ! initialize mesh arrays
179 this%dynvis = marray_base()
180 this%kinvis = marray_base()
181 this%bulkvis = marray_base()
182 this%dynvis%data1d(:) = this%dynconst
183 this%kinvis%data1d(:) = 0.0
184 this%bulkvis%data1d(:) = this%bulkconst
185
186 ! do the model specific initialization
187 SELECT CASE(this%viscosity%GetType())
188 CASE(molecular)
189 ! do nothing
190 CASE(alpha,beta,powerlaw)
191 ! compute the projector to determine the local angular velocity
192 ! ATTENTION: it is always assumed that the fluid rotates around the z-axis
193 SELECT TYPE(phys => physics)
194 CLASS IS(physics_eulerisotherm)
195 ALLOCATE(this%ephir,this%ephir_tmp,stat=err)
196 IF (err.NE.0) &
197 CALL this%Error("InitSources_viscosity","Memory allocation failed.")
198 this%ephir_tmp = marray_base(3)
199 this%ephir = marray_base(3)
200 ! compute ephi/r (r: distance to z-axis, ephi: azimuthal unit vector)
201 ! set cartesian vector components using cartesian coordinates
202 this%ephir%data4d(:,:,:,1) = -mesh%cart%bcenter(:,:,:,2) / mesh%radius%bcenter(:,:,:)**2
203 this%ephir%data4d(:,:,:,2) = mesh%cart%bcenter(:,:,:,1) / mesh%radius%bcenter(:,:,:)**2
204 this%ephir%data4d(:,:,:,3) = 0.0
205 ! convert to curvilinear vector components
206 CALL mesh%geometry%Convert2Curvilinear(mesh%curv%bcenter,this%ephir%data4d,this%ephir_tmp%data4d)
207 ! reduce vector components to match Physics%VDIM
208 IF (phys%VDIM.LT.3) THEN
209 DEALLOCATE(this%ephir,stat=err)
210 IF (err.EQ.0) THEN
211 ALLOCATE(this%ephir,stat=err)
212 IF (err.EQ.0) this%ephir = marray_base(phys%VDIM)
213 END IF
214 IF (err.NE.0) &
215 CALL this%Error("InitSources_viscosity","Memory reallocation failed.")
216 END IF
217 ! copy vector components actually used
218 l = 1
219 DO k=1,3
220 IF (btest(mesh%VECTOR_COMPONENTS,k-1)) THEN
221 this%ephir%data2d(:,l) = this%ephir_tmp%data2d(:,k)
222 l = l + 1
223 END IF
224 END DO
225 ! destroy temporary marray
226 DEALLOCATE(this%ephir_tmp)
227 CLASS DEFAULT
228 CALL this%Error("InitSources_viscosity",&
229 "Physics not supported for alpha-/beta-viscosity")
230 END SELECT
231 CASE(alpha_alt)
233 CALL this%Error("InitSources_viscosity",&
234 "alternative alpha-viscosity currently not supported")
235! IF (Physics%VDIM.NE.2) &
236! CALL this%Error("InitSources_viscosity",&
237! "alternative alpha-viscosity works only for flat disks")
238 ! IF (.NOT.Timedisc%always_update_bccsound) &
239 ! CALL Error(this,"InitSources_viscosity","always_update_bccsound must be enabled in timedisc")
240 CASE DEFAULT
241 CALL this%Error("InitSources_viscosity",&
242 "Viscosity prescription not supported.")
243 END SELECT
244 ! set initial time to negative value;
245 ! this guarantees update of viscosity for time t=0.0
246 this%time = -1.0
247 ! set stress tensor arrays to zero
248 this%btxx(:,:,:) = 0.0
249 this%btyy(:,:,:) = 0.0
250 this%btzz(:,:,:) = 0.0
251 this%btxy(:,:,:) = 0.0
252 this%btxz(:,:,:) = 0.0
253 this%btyz(:,:,:) = 0.0
254
255 CALL this%SetOutput(mesh,config,io)
256
257 CALL this%InfoSources()
258 END SUBROUTINE initsources
259
260 SUBROUTINE setoutput(this,Mesh,config,IO)
261 IMPLICIT NONE
262 !------------------------------------------------------------------------!
263 CLASS(sources_viscosity),INTENT(IN) :: this
264 CLASS(mesh_base) ,INTENT(IN) :: Mesh
265 TYPE(dict_typ),POINTER :: config,IO
266 !------------------------------------------------------------------------!
267 INTEGER :: valwrite
268 !------------------------------------------------------------------------!
269 ! Hier wird der Stresstensor per default in 3D ausgegeben! Vorher nur in 2D plus abhängig
270 ! von der Physik auch die Komponenten der dritten Dimension
271 CALL getattr(config, "output/stress", valwrite, 0)
272 IF (valwrite .EQ. 1) THEN
273 CALL setattr(io, "stress:xx", this%btxx(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
274 CALL setattr(io, "stress:xy", this%btxy(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
275 CALL setattr(io, "stress:xz", this%btxz(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
276 CALL setattr(io, "stress:yy", this%btyy(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
277 CALL setattr(io, "stress:yz", this%btyz(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
278 CALL setattr(io, "stress:zz", this%btzz(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
279 END IF
280
281 CALL getattr(config, "output/dynvis", valwrite, 0)
282 IF (valwrite .EQ. 1) &
283 CALL setattr(io,"dynvis", this%dynvis%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
284
285 CALL getattr(config, "output/kinvis", valwrite, 0)
286 IF (valwrite .EQ. 1) &
287 CALL setattr(io, "kinvis", this%kinvis%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
288
289 CALL getattr(config, "output/bulkvis", valwrite, 0)
290 IF (valwrite .EQ. 1) &
291 CALL setattr(io, "bulkvis", this%bulkvis%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
292
293 END SUBROUTINE setoutput
294
295
296 SUBROUTINE infosources(this)
297 IMPLICIT NONE
298 !------------------------------------------------------------------------!
299 CLASS(sources_viscosity),INTENT(IN) :: this
300 !------------------------------------------------------------------------!
301 CHARACTER(LEN=32) :: dynconst_str,bulkconst_str
302 !------------------------------------------------------------------------!
303 WRITE (dynconst_str,'(ES9.2)') this%dynconst
304 CALL this%Info(" viscosity model: " // this%viscosity%GetName())
305 SELECT CASE(this%viscosity%GetType())
306 CASE(molecular)
307 WRITE (bulkconst_str,'(ES9.2)') this%bulkconst
308 CALL this%Info(" dynamic viscosity: " // trim(dynconst_str))
309 CALL this%Info(" bulk viscosity: " // trim(bulkconst_str))
310 CASE(alpha,alpha_alt)
311 CALL this%Info(" alpha: " // trim(dynconst_str))
312 CASE(beta)
313 CALL this%Info(" beta: " // trim(dynconst_str))
314 CASE(powerlaw)
315 WRITE (bulkconst_str,'(ES9.2)') this%power
316 CALL this%Info(" coupling constant: " // trim(dynconst_str))
317 CALL this%Info(" exponent: " // trim(bulkconst_str))
318 END SELECT
319 END SUBROUTINE infosources
320
321
356 SUBROUTINE updateviscosity(this,Mesh,Physics,Fluxes,time,pvar,cvar)
357 USE physics_eulerisotherm_mod, ONLY : physics_eulerisotherm, statevector_eulerisotherm
358 USE physics_euler_mod, ONLY : physics_euler, statevector_euler
359 IMPLICIT NONE
360 !------------------------------------------------------------------------!
361 CLASS(sources_viscosity),INTENT(INOUT) :: this
362 CLASS(mesh_base),INTENT(IN) :: Mesh
363 CLASS(physics_base),INTENT(INOUT) :: Physics
364 CLASS(fluxes_base),INTENT(IN) :: Fluxes
365 REAL,INTENT(IN) :: time
366 CLASS(marray_compound),INTENT(INOUT) :: pvar,cvar
367 !------------------------------------------------------------------------!
368 INTEGER :: l
369 !------------------------------------------------------------------------!
370 ! only update viscosity if time has changed
371 IF ((time.NE.this%time) .OR. (time .EQ. 0.)) THEN
372 ! Set the dynamic viscosity
373 SELECT CASE(this%viscosity%GetType())
374 CASE(molecular)
375 ! do nothing, since it has already been initialized
376 ! in InitSources_viscosity
377 CASE(alpha)
378 ! compute Shakura-Sunyaev type alpha viscosity
379 ! account for non-inertial frames by adding the frame angular velocity
380 SELECT TYPE(p => pvar)
381 CLASS IS(statevector_eulerisotherm)
382 SELECT TYPE(phys => physics)
383 CLASS IS(physics_eulerisotherm)
384 this%dynvis%data1d(:) = etafkt_alpha(this%dynconst,p%density%data1d(:),&
385 phys%bccsound%data1d(:),omega(this%ephir,p%velocity)+mesh%Omega)
386 END SELECT
387 END SELECT
388 CASE(beta)
389 ! compute Duschl type beta viscosity
390 ! account for non-inertial frames by adding the frame angular velocity
391 SELECT TYPE(p => pvar)
392 CLASS IS(statevector_eulerisotherm)
393 this%dynvis%data1d(:) = etafkt_beta(this%dynconst,mesh%radius%data2d(:,2), &
394 p%density%data1d(:),omega(this%ephir,p%velocity)+mesh%Omega)
395 END SELECT
396 CASE(powerlaw)
397 ! kinematic viscosity scales with power of specific angular momentum
398 ! account for non-inertial frames by adding the frame angular velocity
399 SELECT TYPE(p => pvar)
400 CLASS IS(statevector_eulerisotherm)
401 this%dynvis%data1d(:) = etafkt_powerlaw(this%dynconst,this%power, &
402 mesh%radius%data2d(:,2),p%density%data1d(:),omega(this%ephir,p%velocity)+mesh%Omega)
403 END SELECT
404! TODO: height of the disk is calculated in another source module, which is not translated yet. Therefore, ALPHA_ALT is not
405! working
406! CASE(ALPHA_ALT)
407! ! eta = nu*rho = alpha*cs*h * rho
408! ! compute dynamic viscosity
409! this%dynvis(:,:,:) = this%dynconst * Physics%bccsound%data3d(:,:,:) &
410! * Physics%sources%height(:,:,:) * pvar(:,:,:,Physics%DENSITY)
411 END SELECT
412
413 ! Set the bulk viscosity
414 SELECT CASE(this%viscosity%GetType())
415 CASE(molecular)
416 ! do nothing, since it has already been initialized
417 ! in InitSources_viscosity
418 CASE(alpha,beta,powerlaw,alpha_alt)
419 ! All available parametrized viscosities should be trace free.
420 ! This can be achieved by setting the following:
421 this%bulkvis%data1d(:) = -2./3.*this%dynvis%data1d(:)
422 END SELECT
423
424 ! update time
425 this%time = time
426 END IF
427 CONTAINS
428 ! some elemental functions for computation of the viscosity
429
430 ! compute angular velocity
431 PURE FUNCTION omega(ephir,velocity)
432 IMPLICIT NONE
433 !----------------------------------------------------------------------!
434 CLASS(marray_base), INTENT(IN) :: ephir,velocity
435 REAL, DIMENSION(SIZE(ephir%data2d,DIM=1)) :: omega
436 !----------------------------------------------------------------------!
437 INTEGER :: l
438 !----------------------------------------------------------------------!
439 ! project velocity on ephi/r -> local angular velocity
440 omega(:) = ephir%data2d(:,1) * velocity%data2d(:,1)
441 DO l=2,ephir%DIMS(1)
442 omega(:) = omega(:) + ephir%data2d(:,l) * velocity%data2d(:,l)
443 END DO
444 END FUNCTION omega
445
446 ! Alpha viscosity
447 ELEMENTAL FUNCTION etafkt_alpha(alpha,density,cs,omega) RESULT(eta)
448 IMPLICIT NONE
449 !----------------------------------------------------------------------!
450 REAL, INTENT(IN) :: alpha,density,cs,omega
451 REAL :: eta
452 !----------------------------------------------------------------------!
453 eta = alpha * density * cs**2 / abs(omega)
454 END FUNCTION etafkt_alpha
455
456 ! Beta viscosity
457 ELEMENTAL FUNCTION etafkt_beta(beta,r,density,omega) RESULT(eta)
458 IMPLICIT NONE
459 !----------------------------------------------------------------------!
460 REAL, INTENT(IN) :: beta,r,density,omega
461 REAL :: eta
462 !----------------------------------------------------------------------!
463 ! nu = beta * r * v_phi = beta * r^2 * omega
464 ! eta = rho * nu
465 eta = beta * density * r**2 * abs(omega)
466 END FUNCTION etafkt_beta
467
468 ! power law viscosity
469 ELEMENTAL FUNCTION etafkt_powerlaw(beta,q,r,density,omega) RESULT(eta)
470 IMPLICIT NONE
471 !----------------------------------------------------------------------!
472 REAL, INTENT(IN) :: beta,q,r,density,omega
473 REAL :: eta
474 !----------------------------------------------------------------------!
475 eta = beta * density * (r**2 * abs(omega))**q
476 END FUNCTION etafkt_powerlaw
477
478 END SUBROUTINE updateviscosity
479
480
481 SUBROUTINE externalsources(this,Mesh,Physics,Fluxes,Sources,time,dt,pvar,cvar,sterm)
482 USE physics_eulerisotherm_mod, ONLY : physics_eulerisotherm
483 USE physics_euler_mod, ONLY : physics_euler
484 IMPLICIT NONE
485 !------------------------------------------------------------------------!
486 CLASS(sources_viscosity),INTENT(INOUT) :: this
487 CLASS(mesh_base),INTENT(IN) :: Mesh
488 CLASS(physics_base),INTENT(INOUT) :: Physics
489 CLASS(fluxes_base),INTENT(IN) :: Fluxes
490 CLASS(sources_base), INTENT(INOUT) :: Sources
491 REAL,INTENT(IN) :: time, dt
492 CLASS(marray_compound),INTENT(INOUT) :: pvar,cvar,sterm
493 !------------------------------------------------------------------------!
494 CALL this%UpdateViscosity(mesh,physics,fluxes,time,pvar,cvar)
495 SELECT TYPE(phys => physics)
496 CLASS IS(physics_eulerisotherm)
497 CALL phys%CalcStresses(mesh,pvar,this%dynvis,this%bulkvis, &
498 this%btxx,this%btxy,this%btxz,this%btyy,this%btyz,this%btzz)
499 CALL phys%ViscositySources(mesh,pvar,this%btxx,this%btxy,this%btxz, &
500 this%btyy,this%btyz,this%btzz,sterm)
501 END SELECT
502 END SUBROUTINE externalsources
503
504
511 SUBROUTINE calctimestep(this,Mesh,Physics,Fluxes,pvar,cvar,time,dt,dtcause)
512 USE physics_eulerisotherm_mod, ONLY : physics_eulerisotherm, statevector_eulerisotherm
513 USE physics_euler_mod, ONLY : physics_euler, statevector_euler
514 IMPLICIT NONE
515 !------------------------------------------------------------------------!
516 CLASS(sources_viscosity),INTENT(INOUT) :: this
517 CLASS(mesh_base),INTENT(IN) :: Mesh
518 CLASS(physics_base),INTENT(INOUT) :: Physics
519 CLASS(fluxes_base),INTENT(IN) :: Fluxes
520 CLASS(marray_compound), INTENT(INOUT) :: pvar,cvar
521 REAL, INTENT(IN) :: time
522 REAL, INTENT(INOUT) :: dt
523 INTEGER, INTENT(OUT) :: dtcause
524 !------------------------------------------------------------------------!
525 REAL :: invdt
526 !------------------------------------------------------------------------!
527 ! update dynamic and bulk viscosity
528 CALL this%UpdateViscosity(mesh,physics,fluxes,time,pvar,cvar)
529
530 ! compute kinematic viscosity
531 SELECT TYPE(p => pvar)
532 CLASS IS(statevector_eulerisotherm)
533 this%kinvis%data1d(:) = this%dynvis%data1d(:) / p%density%data1d(:)
534 END SELECT
535
536 ! compute maximum of inverse time step for combined advection-diffusion transport
537 IF ((mesh%JNUM.EQ.1).AND.(mesh%KNUM.EQ.1)) THEN
538 ! 1D, only x-direction
539 invdt = maxval((max(fluxes%maxwav%data2d(:,1),-fluxes%minwav%data2d(:,1)) &
540 + 2*this%kinvis%data1d(:) / mesh%dlx%data1d(:)) / mesh%dlx%data1d(:))
541 ELSE IF ((mesh%INUM.EQ.1).AND.(mesh%KNUM.EQ.1)) THEN
542 ! 1D, only y-direction
543 invdt = maxval((max(fluxes%maxwav%data2d(:,1),-fluxes%minwav%data2d(:,1)) &
544 + 2*this%kinvis%data1d(:) / mesh%dly%data1d(:)) / mesh%dly%data1d(:))
545 ELSE IF ((mesh%INUM.EQ.1).AND.(mesh%JNUM.EQ.1)) THEN
546 ! 1D, only z-direction
547 invdt = maxval((max(fluxes%maxwav%data2d(:,1),-fluxes%minwav%data2d(:,1)) &
548 + 2*this%kinvis%data1d(:) / mesh%dlz%data1d(:)) / mesh%dlz%data1d(:))
549 ELSE IF ((mesh%INUM.GT.1).AND.(mesh%JNUM.GT.1).AND.(mesh%KNUM.EQ.1)) THEN
550 ! 2D, x-y-plane
551 invdt = maxval((max(fluxes%maxwav%data2d(:,1),-fluxes%minwav%data2d(:,1)) &
552 + 2*this%kinvis%data1d(:) / mesh%dlx%data1d(:)) / mesh%dlx%data1d(:) &
553 + (max(fluxes%maxwav%data2d(:,2),-fluxes%minwav%data2d(:,2)) &
554 + 2*this%kinvis%data1d(:) / mesh%dly%data1d(:)) / mesh%dly%data1d(:))
555 ELSE IF ((mesh%INUM.GT.1).AND.(mesh%KNUM.GT.1).AND.(mesh%JNUM.EQ.1)) THEN
556 ! 2D, x-z-plane
557 invdt = maxval((max(fluxes%maxwav%data2d(:,1),-fluxes%minwav%data2d(:,1)) &
558 + 2*this%kinvis%data1d(:) / mesh%dlx%data1d(:)) / mesh%dlx%data1d(:) &
559 + (max(fluxes%maxwav%data2d(:,2),-fluxes%minwav%data2d(:,2)) &
560 + 2*this%kinvis%data1d(:) / mesh%dlz%data1d(:)) / mesh%dlz%data1d(:))
561 ELSE IF ((mesh%JNUM.GT.1).AND.(mesh%KNUM.GT.1).AND.(mesh%INUM.EQ.1)) THEN
562 ! 2D, y-z-plane
563 invdt = maxval((max(fluxes%maxwav%data2d(:,1),-fluxes%minwav%data2d(:,1)) &
564 + 2*this%kinvis%data1d(:) / mesh%dly%data1d(:)) / mesh%dly%data1d(:) &
565 + (max(fluxes%maxwav%data2d(:,2),-fluxes%minwav%data2d(:,2)) &
566 + 2*this%kinvis%data1d(:) / mesh%dlz%data1d(:)) / mesh%dlz%data1d(:))
567 ELSE
568 ! full 3D
569 invdt = maxval((max(fluxes%maxwav%data2d(:,1),-fluxes%minwav%data2d(:,1)) &
570 + 2*this%kinvis%data1d(:) / mesh%dlx%data1d(:)) / mesh%dlx%data1d(:) &
571 + (max(fluxes%maxwav%data2d(:,2),-fluxes%minwav%data2d(:,2)) &
572 + 2*this%kinvis%data1d(:) / mesh%dly%data1d(:)) / mesh%dly%data1d(:) &
573 + (max(fluxes%maxwav%data2d(:,3),-fluxes%minwav%data2d(:,3)) &
574 + 2*this%kinvis%data1d(:) / mesh%dlz%data1d(:)) / mesh%dlz%data1d(:))
575 END IF
576 IF (invdt.GT.tiny(invdt)) THEN
577 dt = this%cvis / invdt
578 ELSE
579 dt = huge(dt)
580 END IF
581 END SUBROUTINE calctimestep
582
584 SUBROUTINE finalize(this)
585 IMPLICIT NONE
586 !------------------------------------------------------------------------!
587 TYPE(sources_viscosity),INTENT(INOUT) :: this
588 !------------------------------------------------------------------------!
589 IF (ALLOCATED(this%ephir)) DEALLOCATE(this%ephir)
590 DEALLOCATE(this%dynvis,this%kinvis,this%bulkvis, &
591 this%btxx,this%btyy,this%btzz,this%btxy,this%btxz,this%btyz)
592 END SUBROUTINE finalize
593
594END MODULE sources_viscosity_mod
common_dict
Dictionary for generic data types.
Definition: common_dict.f90:61
fluxes_base_mod
base module for numerical flux functions
Definition: fluxes_base.f90:39
geometry_base_mod
base class for geometrical properties
Definition: geometry_base.f90:36
geometry_cylindrical_mod
defines properties of a 3D cylindrical mesh
Definition: geometry_cylindrical.f90:36
geometry_logcylindrical_mod
defines properties of a 3D logcylindrical mesh
Definition: geometry_logcylindrical.f90:33
logging_base_mod
Basic fosite module.
Definition: logging_base.f90:36
logging_base_mod::finalize
subroutine finalize(this)
Destructor of logging_base class.
Definition: logging_base.f90:160
marray_base_mod
base class for mesh arrays
Definition: marray_base.f90:36
marray_compound_mod
derived class for compound of mesh arrays
Definition: marray_compound.f90:35
mesh_base_mod
basic mesh module
Definition: mesh_base.f90:72
mesh_base_mod::setoutput
subroutine setoutput(this, config, IO)
Setup mesh fields for i/o.
Definition: mesh_base.f90:828
physics_base_mod
Basic physics module.
Definition: physics_base.f90:55
physics_euler_mod
physics module for 1D,2D and 3D non-isothermal Euler equations
Definition: physics_euler.f90:41
physics_eulerisotherm_mod
physics module for 1D,2D and 3D isothermal Euler equations
Definition: physics_eulerisotherm.f90:44
sources_base_mod
generic source terms module providing functionaly common to all source terms
Definition: sources_base.f90:42
sources_viscosity_mod
computes momentum and energy sources due to shear stresses
Definition: sources_viscosity.f90:61
sources_viscosity_mod::infosources
subroutine infosources(this)
Definition: sources_viscosity.f90:297
sources_viscosity_mod::alpha_alt
integer, parameter, public alpha_alt
Definition: sources_viscosity.f90:79
sources_viscosity_mod::powerlaw
integer, parameter, public powerlaw
Definition: sources_viscosity.f90:78
sources_viscosity_mod::molecular
integer, parameter, public molecular
Definition: sources_viscosity.f90:75
sources_viscosity_mod::beta
integer, parameter, public beta
Definition: sources_viscosity.f90:77
sources_viscosity_mod::source_name
character(len=32), parameter source_name
Definition: sources_viscosity.f90:73
sources_viscosity_mod::initsources
subroutine initsources(this, Mesh, Physics, Fluxes, config, IO)
Definition: sources_viscosity.f90:119
sources_viscosity_mod::calctimestep
subroutine calctimestep(this, Mesh, Physics, Fluxes, pvar, cvar, time, dt, dtcause)
compute time step limit for advection-diffusion problems
Definition: sources_viscosity.f90:512
sources_viscosity_mod::externalsources
subroutine externalsources(this, Mesh, Physics, Fluxes, Sources, time, dt, pvar, cvar, sterm)
Definition: sources_viscosity.f90:482
sources_viscosity_mod::viscosity_name
character(len=32), dimension(5), parameter, public viscosity_name
Definition: sources_viscosity.f90:80
sources_viscosity_mod::alpha
integer, parameter, public alpha
Definition: sources_viscosity.f90:76
sources_viscosity_mod::updateviscosity
subroutine updateviscosity(this, Mesh, Physics, Fluxes, time, pvar, cvar)
updates dynamic and bulk viscosity
Definition: sources_viscosity.f90:357
etafkt_beta
elemental real function etafkt_beta(beta, r, density, omega)
Definition: sources_viscosity.f90:458
etafkt_alpha
elemental real function etafkt_alpha(alpha, density, cs, omega)
Definition: sources_viscosity.f90:448
omega
pure real function, dimension(size(ephir%data2d, dim=1)) omega(ephir, velocity)
Definition: sources_viscosity.f90:432
etafkt_powerlaw
elemental real function etafkt_powerlaw(beta, q, r, density, omega)
Definition: sources_viscosity.f90:470
common_dict::dict_typ
Definition: common_dict.f90:123
fluxes_base_mod::fluxes_base
Definition: fluxes_base.f90:61
geometry_cylindrical_mod::geometry_cylindrical
Definition: geometry_cylindrical.f90:41
geometry_logcylindrical_mod::geometry_logcylindrical
Definition: geometry_logcylindrical.f90:39
logging_base_mod::logging_base
common data structure
Definition: logging_base.f90:56
marray_base_mod::marray_base
basic mesh array class
Definition: marray_base.f90:69
marray_compound_mod::marray_compound
data types and methods
Definition: marray_compound.f90:46
mesh_base_mod::mesh_base
mesh data structure
Definition: mesh_base.f90:122
physics_base_mod::physics_base
Definition: physics_base.f90:69
physics_euler_mod::physics_euler
Definition: physics_euler.f90:54
physics_euler_mod::statevector_euler
Definition: physics_euler.f90:105
physics_eulerisotherm_mod::physics_eulerisotherm
Definition: physics_eulerisotherm.f90:56
physics_eulerisotherm_mod::statevector_eulerisotherm
Definition: physics_eulerisotherm.f90:122
sources_base_mod::sources_base
Definition: sources_base.f90:57
sources_viscosity_mod::sources_viscosity
Definition: sources_viscosity.f90:87
Generated by   doxygen 1.9.4