• R/O
  • SSH

musubi: 提交

Lattice Boltzmann solver in APES. Now found at https://github.com/apes-suite/musubi


Commit MetaInfo

修訂35d71a54ed5d6c224e06bea9f84ca31000e50900 (tree)
時間2022-08-30 19:25:45
作者Harald Klimach <harald.klimach@dlr....>
CommiterHarald Klimach

Log Message

Fixed compiler warnings.

Change Summary

差異

diff -r b6956abee5af -r 35d71a54ed5d source/bc/mus_bc_fluid_module.fpp
--- a/source/bc/mus_bc_fluid_module.fpp Thu Jun 09 14:16:03 2022 +0200
+++ b/source/bc/mus_bc_fluid_module.fpp Tue Aug 30 12:25:45 2022 +0200
@@ -2417,19 +2417,13 @@
24172417 !> mixture info
24182418 type(mus_mixture_type), intent(in) :: mixture
24192419 ! -------------------------------------------------------------------- !
2420- real(kind=rk) :: rhoDef(globBC%nElems(iLevel)) ! Density on boundary element
2421- real(kind=rk) :: rhoF ! density on fluid element
2422- real(kind=rk) :: uxF(globBC%nElems(iLevel)*3) ! Velocity on fluid element
2423- real(kind=rk) :: uxN(globBC%nElems(iLevel)*3) ! Velocity on neighbor
24242420 ! Equilibrium Plus for bounceback
24252421 real(kind=rk) :: fEqPlus
2426- real(kind=rk) :: inv_rho_phy
24272422 real(kind=rk) :: fTmp( nScalars )
24282423 real(kind=rk) :: fTmpPre( nScalars )
24292424 real(kind=rk) :: fTmpPre_2( nScalars )
2430- real(kind=rk) :: omega ! Relaxation parameter
2431- integer :: iDir, iElem, QQ, invDir, elemPos
2432- integer :: posInBuffer, bcPress_pos
2425+ integer :: iDir, iElem, QQ, elemPos
2426+ integer :: posInBuffer
24332427 real(kind=rk) :: swap
24342428 real(kind=rk) :: L1, L2, L3, L4, L5
24352429 real(kind=rk) :: dudx, dvdx, dwdx, drhodx
@@ -2710,19 +2704,13 @@
27102704 !> mixture info
27112705 type(mus_mixture_type), intent(in) :: mixture
27122706 ! -------------------------------------------------------------------- !
2713- real(kind=rk) :: rhoDef(globBC%nElems(iLevel)) ! Density on boundary element
2714- real(kind=rk) :: rhoF ! density on fluid element
2715- real(kind=rk) :: uxF(globBC%nElems(iLevel)*3) ! Velocity on fluid element
2716- real(kind=rk) :: uxN(globBC%nElems(iLevel)*3) ! Velocity on neighbor
27172707 ! Equilibrium Plus for bounceback
27182708 real(kind=rk) :: fEqPlus
2719- real(kind=rk) :: inv_rho_phy
27202709 real(kind=rk) :: fTmp( nScalars )
27212710 real(kind=rk) :: fTmpPre( nScalars )
27222711 real(kind=rk) :: fTmpPre_2( nScalars )
2723- real(kind=rk) :: omega ! Relaxation parameter
2724- integer :: iDir, iElem, QQ, invDir, elemPos
2725- integer :: posInBuffer, bcPress_pos
2712+ integer :: iDir, iElem, QQ, elemPos
2713+ integer :: posInBuffer
27262714 real(kind=rk) :: swap
27272715 real(kind=rk) :: L1, L2, L3, L4, L5
27282716 real(kind=rk) :: dudx, dvdx, dwdx, drhodx
diff -r b6956abee5af -r 35d71a54ed5d source/compute/mus_compute_d2q9_module.fpp
--- a/source/compute/mus_compute_d2q9_module.fpp Thu Jun 09 14:16:03 2022 +0200
+++ b/source/compute/mus_compute_d2q9_module.fpp Tue Aug 30 12:25:45 2022 +0200
@@ -82,37 +82,37 @@
8282 ! moment(iDir) = sum( PDF(:) * MMtrD2Q9(:,iDir) )
8383 ! end do
8484 ! W S E N SW NW SE NE 0
85- real(kind=rk), dimension(9,9),parameter :: MMtrD2Q9 = &
86- reshape((/ &
87- 1._rk, 1._rk, 1._rk, 1._rk, 1._rk, 1._rk, 1._rk, 1._rk, 1._rk,& !
88- -1._rk, -1._rk, -1._rk, -1._rk, 2._rk, 2._rk, 2._rk, 2._rk, -4._rk,& !mom2
89- -2._rk, -2._rk, -2._rk, -2._rk, 1._rk, 1._rk, 1._rk, 1._rk, 4._rk,& !mom3
90- -1._rk, 0._rk, 1._rk, 0._rk, -1._rk, -1._rk, 1._rk, 1._rk, 0._rk,&
91- 2._rk, 0._rk, -2._rk, 0._rk, -1._rk, -1._rk, 1._rk, 1._rk, 0._rk,& !mom5
92- 0._rk, -1._rk, 0._rk, 1._rk, -1._rk, 1._rk, -1._rk, 1._rk, 0._rk,&
93- 0._rk, 2._rk, 0._rk, -2._rk, -1._rk, 1._rk, -1._rk, 1._rk, 0._rk,& !mom7
94- 1._rk, -1._rk, 1._rk, -1._rk, 0._rk, 0._rk, 0._rk, 0._rk, 0._rk,& !mom8
95- 0._rk, 0._rk, 0._rk, 0._rk, 1._rk, -1._rk, -1._rk, 1._rk, 0._rk & !mom9
96- /),(/9,9/))
85+!! real(kind=rk), dimension(9,9),parameter :: MMtrD2Q9 = &
86+!! reshape((/ &
87+!! 1._rk, 1._rk, 1._rk, 1._rk, 1._rk, 1._rk, 1._rk, 1._rk, 1._rk,& !
88+!! -1._rk, -1._rk, -1._rk, -1._rk, 2._rk, 2._rk, 2._rk, 2._rk, -4._rk,& !mom2
89+!! -2._rk, -2._rk, -2._rk, -2._rk, 1._rk, 1._rk, 1._rk, 1._rk, 4._rk,& !mom3
90+!! -1._rk, 0._rk, 1._rk, 0._rk, -1._rk, -1._rk, 1._rk, 1._rk, 0._rk,&
91+!! 2._rk, 0._rk, -2._rk, 0._rk, -1._rk, -1._rk, 1._rk, 1._rk, 0._rk,& !mom5
92+!! 0._rk, -1._rk, 0._rk, 1._rk, -1._rk, 1._rk, -1._rk, 1._rk, 0._rk,&
93+!! 0._rk, 2._rk, 0._rk, -2._rk, -1._rk, 1._rk, -1._rk, 1._rk, 0._rk,& !mom7
94+!! 1._rk, -1._rk, 1._rk, -1._rk, 0._rk, 0._rk, 0._rk, 0._rk, 0._rk,& !mom8
95+!! 0._rk, 0._rk, 0._rk, 0._rk, 1._rk, -1._rk, -1._rk, 1._rk, 0._rk & !mom9
96+!! /),(/9,9/))
9797
9898 ! D2Q9 MRT moment -> pdf transformation matrix
9999 ! How to use:
100100 ! do iDir = 1, QQ
101101 ! pdf(iDir) = sum( moment(:) * MMIvD2Q9(iDir,:) )
102102 ! end do
103- real(kind=rk), dimension(9,9), parameter :: MMIvD2Q9 = &
104- reshape((/ &
105-! 1 2 3 4 5 6 7 8 9
106- div1_9, -div1_36, -div1_18, -div1_6, div1_6, 0._rk, 0._rk, div1_4, 0._rk, & ! W
107- div1_9, -div1_36, -div1_18, 0._rk, 0._rk, -div1_6, div1_6, -div1_4, 0._rk, & ! S
108- div1_9, -div1_36, -div1_18, div1_6, -div1_6, 0._rk, 0._rk, div1_4, 0._rk, & ! E
109- div1_9, -div1_36, -div1_18, 0._rk, 0._rk, div1_6, -div1_6, -div1_4, 0._rk, & ! N
110- div1_9, div1_18, div1_36, -div1_6, -div1_12, -div1_6, -div1_12, 0._rk, div1_4, & ! SW
111- div1_9, div1_18, div1_36, -div1_6, -div1_12, div1_6, div1_12, 0._rk, -div1_4, & ! NW
112- div1_9, div1_18, div1_36, div1_6, div1_12, -div1_6, -div1_12, 0._rk, -div1_4, & ! SE
113- div1_9, div1_18, div1_36, div1_6, div1_12, div1_6, div1_12, 0._rk, div1_4, & ! NE
114- div1_9, -div1_9, div1_9, 0._rk, 0._rk, 0._rk, 0._rk, 0._rk, 0._rk & ! 0
115- /),(/9,9/), order=(/2,1/))
103+!! real(kind=rk), dimension(9,9), parameter :: MMIvD2Q9 = &
104+!! reshape((/ &
105+!!! 1 2 3 4 5 6 7 8 9
106+!! div1_9, -div1_36, -div1_18, -div1_6, div1_6, 0._rk, 0._rk, div1_4, 0._rk, & ! W
107+!! div1_9, -div1_36, -div1_18, 0._rk, 0._rk, -div1_6, div1_6, -div1_4, 0._rk, & ! S
108+!! div1_9, -div1_36, -div1_18, div1_6, -div1_6, 0._rk, 0._rk, div1_4, 0._rk, & ! E
109+!! div1_9, -div1_36, -div1_18, 0._rk, 0._rk, div1_6, -div1_6, -div1_4, 0._rk, & ! N
110+!! div1_9, div1_18, div1_36, -div1_6, -div1_12, -div1_6, -div1_12, 0._rk, div1_4, & ! SW
111+!! div1_9, div1_18, div1_36, -div1_6, -div1_12, div1_6, div1_12, 0._rk, -div1_4, & ! NW
112+!! div1_9, div1_18, div1_36, div1_6, div1_12, -div1_6, -div1_12, 0._rk, -div1_4, & ! SE
113+!! div1_9, div1_18, div1_36, div1_6, div1_12, div1_6, div1_12, 0._rk, div1_4, & ! NE
114+!! div1_9, -div1_9, div1_9, 0._rk, 0._rk, 0._rk, 0._rk, 0._rk, 0._rk & ! 0
115+!! /),(/9,9/), order=(/2,1/))
116116
117117
118118 contains
diff -r b6956abee5af -r 35d71a54ed5d source/derived/mus_auxFieldVar_module.fpp
--- a/source/derived/mus_auxFieldVar_module.fpp Thu Jun 09 14:16:03 2022 +0200
+++ b/source/derived/mus_auxFieldVar_module.fpp Tue Aug 30 12:25:45 2022 +0200
@@ -1971,7 +1971,7 @@
19711971 real(kind=rk) :: inv_rho_phy, inv_vel_phy
19721972 integer :: iElem, nElems, elemOff
19731973 real(kind=rk) :: spongeField(fun%elemLvl(iLevel)%nElems)
1974- real(kind=rk) :: dens, vel(3), spongeForce(27)
1974+ real(kind=rk) :: dens, vel(3)
19751975 real(kind=rk) :: dens_ref, vel_ref(3), sponDens, sponVel(3)
19761976 ! ------------------------------------------------------------------------ !
19771977 ! position of density and velocity field in auxField
diff -r b6956abee5af -r 35d71a54ed5d source/derived/mus_derQuanIncomp_module.fpp
--- a/source/derived/mus_derQuanIncomp_module.fpp Thu Jun 09 14:16:03 2022 +0200
+++ b/source/derived/mus_derQuanIncomp_module.fpp Tue Aug 30 12:25:45 2022 +0200
@@ -145,6 +145,11 @@
145145 public :: deriveAuxIncomp_fromState
146146 public :: deriveEquilIncomp_fromAux
147147
148+ public :: deriveVelocityincomp
149+ public :: deriveVelocityincomp_fromIndex
150+ public :: mus_deriveVelocityincomp
151+ public :: mus_deriveKeincomp
152+
148153 ! source variables
149154 public :: derive_absorbLayerIncomp
150155 public :: applySrc_absorbLayerIncomp
diff -r b6956abee5af -r 35d71a54ed5d source/derived/mus_derQuanIsothermAcEq_module.fpp
--- a/source/derived/mus_derQuanIsothermAcEq_module.fpp Thu Jun 09 14:16:03 2022 +0200
+++ b/source/derived/mus_derQuanIsothermAcEq_module.fpp Tue Aug 30 12:25:45 2022 +0200
@@ -466,68 +466,6 @@
466466 end subroutine derVelocityIsothermAcEq
467467 ! ****************************************************************************** !
468468
469-! ****************************************************************************** !
470- !> Initiates the calculation of Velocity.
471- !! This routine sets the function Pointer for velocity calcualtion and calls
472- !! the generice get Value of Index routine
473- !!
474- !! The interface has to comply to the abstract interface
475- !! [[tem_varSys_module:tem_varSys_proc_getValOfIndex]].
476- !!
477- recursive subroutine derVelocityIsothermAcEq_fromIndex( fun, varSys, time, &
478- & iLevel, idx, idxLen, nVals, res )
479- !> Description of the method to obtain the variables, here some preset
480- !! values might be stored, like the space time function to use or the
481- !! required variables.
482- class(tem_varSys_op_type), intent(in) :: fun
483-
484- !> The variable system to obtain the variable from.
485- type(tem_varSys_type), intent(in) :: varSys
486-
487- !> Point in time at which to evaluate the variable.
488- type(tem_time_type), intent(in) :: time
489-
490- !> Level on which values are requested
491- integer, intent(in) :: iLevel
492-
493- !> Index of points in the growing array and variable val array to
494- !! return.
495- !! Size: most times nVals, if contiguous arrays are used it depends
496- !! on the number of first indices
497- integer, intent(in) :: idx(:)
498-
499- !> With idx as start index in contiguous memory,
500- !! idxLength defines length of each contiguous memory
501- !! Size: dependes on number of first index for contiguous array,
502- !! but the sum of all idxLen is equal to nVals
503- integer, optional, intent(in) :: idxLen(:)
504-
505- !> Number of values to obtain for this variable (vectorized access).
506- integer, intent(in) :: nVals
507-
508- !> Resulting values for the requested variable.
509- !!
510- !! Dimension: n requested entries x nComponents of this variable
511- !! Access: (iElem-1)*fun%nComponents + iComp
512- real(kind=rk), intent(out) :: res(:)
513- ! ---------------------------------------------------------------------------
514- !> Function pointer to perform specific operation.
515- procedure(mus_derive_fromPDF), pointer :: fnCalcPtr
516- ! ---------------------------------------------------------------------------
517- fnCalcPtr => mus_derVelocityIsothermAcEq
518-
519- call mus_generic_varFromPDF_fromIndex( &
520- & fun = fun, &
521- & varSys = varSys, &
522- & time = time, &
523- & iLevel = iLevel, &
524- & idx = idx, &
525- & nVals = nVals, &
526- & fnCalcPtr = fnCalcPtr, &
527- & res = res )
528-
529- end subroutine derVelocityIsothermAcEq_fromIndex
530-! ****************************************************************************** !
531469
532470 ! ****************************************************************************** !
533471 !> Calculate the equlibrium of given elements with the given input state
diff -r b6956abee5af -r 35d71a54ed5d source/derived/mus_derQuanMSLiquid_module.fpp
--- a/source/derived/mus_derQuanMSLiquid_module.fpp Thu Jun 09 14:16:03 2022 +0200
+++ b/source/derived/mus_derQuanMSLiquid_module.fpp Tue Aug 30 12:25:45 2022 +0200
@@ -3462,353 +3462,6 @@
34623462 ! ************************************************************************ !
34633463
34643464
3465- ! ************************************************************************ !
3466- !> Update state with source variable "electric_field"
3467- !! Simuilar to derive routine but it updates the state whereas derive
3468- !! is used for tracking
3469- !! @todo species electricField
3470- !!
3471- !! This subroutine's interface must match the abstract interface definition
3472- !! [[proc_apply_source]] in derived/[[mus_source_var_module]].f90 in order to
3473- !! be callable via [[mus_source_op_type:applySrc]] function pointer.
3474- subroutine applySrc_electricFieldMSLiquid_2order( fun, inState, outState, &
3475- & neigh, auxField, nPdfSize, &
3476- & iLevel, varSys, time, &
3477- & phyConvFac, derVarPos )
3478- ! -------------------------------------------------------------------- !
3479- !> Description of method to apply source terms
3480- class(mus_source_op_type), intent(in) :: fun
3481-
3482- !> input pdf vector
3483- real(kind=rk), intent(in) :: inState(:)
3484-
3485- !> output pdf vector
3486- real(kind=rk), intent(inout) :: outState(:)
3487-
3488- !> connectivity Array corresponding to state vector
3489- integer,intent(in) :: neigh(:)
3490-
3491- !> auxField array
3492- real(kind=rk), intent(in) :: auxField(:)
3493-
3494- !> number of elements in state Array
3495- integer, intent(in) :: nPdfSize
3496-
3497- !> current level
3498- integer, intent(in) :: iLevel
3499-
3500- !> variable system
3501- type(tem_varSys_type), intent(in) :: varSys
3502-
3503- !> Point in time at which to evaluate the variable.
3504- type(tem_time_type), intent(in) :: time
3505-
3506- !> Physics conversion factor for current level
3507- type(mus_convertFac_type), intent(in) :: phyConvFac
3508-
3509- !> position of derived quantities in varsys
3510- type(mus_derVarPos_type), intent(in) :: derVarPos(:)
3511- ! -------------------------------------------------------------------- !
3512- type(mus_varSys_data_type), pointer :: fPtr
3513- type(mus_scheme_type), pointer :: scheme
3514- real(kind=rk) :: electricField(fun%elemLvl(iLevel)%nElems*3)
3515- real(kind=rk) :: EF_elem(3)
3516- integer :: iElem, nElems, iDir, vPos
3517- integer :: iField, iField_2, nFields, nScalars, QQ, nInputStates
3518- !mass density of nSpecies
3519- real(kind=rk), dimension(varSys%nStateVars) :: mass_dens, num_dens, &
3520- & massFrac, moleFrac, molWeight, phi, chargeNr, chargeTerm
3521- real(kind=rk) :: avgMolWeight, paramB, omega
3522- real(kind=rk) :: gasConst, temp0, minMolWeight
3523- real(kind=rk), dimension(3, varSys%nStateVars ) :: diffForce, electricTerm
3524- real(kind=rk), dimension(varSys%nStateVars, varSys%nStateVars) :: matA, &
3525- & invA, resi_coeff
3526- real(kind=rk), dimension(3, varSys%nStateVars) :: first_moments, velocity, &
3527- & eqVel
3528- real(kind=rk), allocatable :: forceTerm(:), omegaMoments(:,:), &
3529- & omegaMomForce(:,:)
3530- real(kind=rk), allocatable :: fneq(:), fneq_om(:), force_om(:)
3531- real(kind=rk) :: velAvg(3), ucxQuadTerm, totMassDens, totNum_dens
3532- real(kind=rk) :: usqr, ucx, feq, charge_dens
3533- integer :: stateVarMap(varSys%nStateVars)
3534- real(kind=rk) :: pdfTmp( varSys%nScalars ) !< temporary local pdf values
3535- real(kind=rk) :: cxcxTerm(3,3), forceQuadTerm, forceUxTerm, diffForceUxTerm
3536- ! -------------------------------------------------------------------- !
3537-!write(dbgUnit(1),*) 'source variable: ', trim(varSys%varname%val(fun%srcTerm_varPos))
3538- ! convert c pointer to solver type fortran pointer
3539- call c_f_pointer( varSys%method%val( fun%srcTerm_varPos )%method_data, &
3540- & fPtr )
3541- scheme => fPtr%solverData%scheme
3542-
3543- ! Number of elements to apply source terms
3544- nElems = fun%elemLvl(iLevel)%nElems
3545-
3546- ! Get electrical force which is refered in config file either its
3547- ! spacetime variable or operation variable
3548- call varSys%method%val(fun%data_varPos)%get_valOfIndex( &
3549- & varSys = varSys, &
3550- & time = time, &
3551- & iLevel = iLevel, &
3552- & idx = fun%elemLvl(iLevel)%idx(1:nElems), &
3553- & nVals = nElems, &
3554- & res = electricField )
3555-
3556- ! convert physical to lattice
3557- electricField = electricField * fPtr%solverData%physics%coulomb0 &
3558- & / fPtr%solverData%physics%fac(iLevel)%force
3559-
3560- ! number of pdf states this source depends on
3561- ! last input is spacetime function so it is neglected
3562- nInputStates = varSys%method%val(fun%srcTerm_varPos)%nInputs - 1
3563- ! constant parameter
3564- nFields = scheme%nFields
3565- QQ = scheme%layout%fStencil%QQ
3566-
3567- nScalars = varSys%nScalars
3568- stateVarMap = scheme%stateVarMap%varPos%val(:)
3569- ! minimum molecular weight
3570- minMolWeight = minval(scheme%field(:)%fieldProp%species%molWeight)
3571- phi = scheme%field(:)%fieldProp%species%molWeigRatio
3572- molWeight = scheme%field(:)%fieldProp%species%molWeight
3573-!write(dbgUnit(1),*) 'avgMolWeight ', avgMolWeight
3574-!write(dbgUnit(1),*) 'ratio ', molWeight/avgMolWeight
3575- ! species charge
3576- chargeNr(:) = scheme%field(:)%fieldProp%species%chargeNr
3577- ! gas constant
3578- gasConst = scheme%mixture%gasConst_R_LB
3579- ! temperature
3580- temp0 = scheme%mixture%temp0LB
3581- omega = scheme%mixture%relaxLvl(iLevel)%omega_diff
3582-
3583- !resistivities
3584- do iField = 1, nFields
3585- resi_coeff(iField,:) = scheme%field(iField)%fieldProp%species%resi_coeff
3586- enddo
3587-
3588- paramB = scheme%mixture%paramB
3589-
3590- allocate(forceTerm(QQ))
3591- allocate(omegaMoments(QQ,QQ))
3592- allocate(omegaMomForce(QQ,QQ))
3593- allocate(fneq(QQ))
3594- allocate(fneq_om(QQ))
3595- allocate(force_om(QQ))
3596-
3597- nodeloop: do iElem = 1, nElems
3598-!write(dbgUnit(1),*) 'iElem', iElem
3599-
3600- mass_dens = 0._rk
3601- first_moments = 0.0_rk
3602- velocity = 0.0_rk
3603- eqVel = 0.0_rk
3604- pdfTmp = 0._rk
3605-
3606- do iField = 1, nFields
3607- ! compute field density and first moments
3608- do iDir = 1, QQ
3609- vPos = varSys%method%val( stateVarMap(iField) )%state_varPos(iDir)
3610- !store all field pdf in single element to local array pdfTmp
3611- pdfTmp( vPos ) = &
3612-& instate( ?FETCH?( iDir, iField, iElem, QQ, nScalars, nPdfSize,neigh ) )
3613- !field density
3614- mass_dens( iField ) = mass_dens( iField ) + pdfTmp( vPos )
3615-
3616- !field momentum (rho*u)
3617- first_moments( 1, iField ) = first_moments( 1, iField ) &
3618- & + pdfTmp( vPos ) * scheme%layout%fStencil%cxDir( 1, iDir )
3619-
3620- first_moments( 2, iField ) = first_moments( 2, iField ) &
3621- & + pdfTmp( vPos ) * scheme%layout%fStencil%cxDir( 2, iDir )
3622-
3623- first_moments( 3, iField ) = first_moments( 3, iField ) &
3624- & + pdfTmp( vPos ) * scheme%layout%fStencil%cxDir( 3, iDir )
3625-!write(dbgunit,*) 's',s,'iDir', iDir, layout%fStencil%cxDir(3,iDir), pdfTmp( vPos ), gqzsigma(s)
3626- enddo
3627- enddo
3628-
3629- !total mass density
3630- totmassDens = sum(mass_dens)
3631-
3632- !mass fraction
3633- massFrac(:) = mass_dens(:)/totmassDens
3634-
3635- ! solve linear system of equation for actual moments
3636- ! number density of all species
3637- num_dens(:) = mass_dens(:)/molWeight(:)
3638-
3639- !total number density
3640- totNum_dens = sum(num_dens(:))
3641-
3642- !mole fraction
3643- moleFrac(:) = num_dens(:)/totNum_dens
3644-
3645- ! mean molecular weight of mixture
3646- avgMolWeight = totmassDens/totNum_dens
3647- !Electric force per element E * Faraday
3648- EF_elem = electricField((iElem-1)*3+1 : iElem*3)
3649- ! n_k * z_k
3650- do iField = 1, nFields
3651- chargeTerm(iField) = num_dens(iField) * chargeNr(iField) &
3652- & * scheme%mixture%faradayLB
3653- end do
3654- charge_dens = sum(chargeTerm)
3655-!write(dbgUnit(1),*) 'massdens ', mass_dens
3656-!write(dbgUnit(1),*) 'numdens ', num_dens
3657-!write(dbgUnit(1),*) 'chargeTerm', chargeTerm
3658-!write(dbgUnit(1),*) 'EF_elem ', EF_elem
3659-!write(dbgUnit(1),*) 'charge_dens ', charge_dens
3660-!write(dbgUnit(1),*) 'fac ', minMolWeight/(gasConst*temp0)
3661- do iField = 1, nFields
3662- diffForce(:, iField) = minMolWeight / (gasConst*temp0) &
3663- & * ( chargeTerm(iField) &
3664- & - massFrac(iField)*charge_dens ) &
3665- & * EF_elem(:)
3666- electricTerm(:, iField) = chargeTerm(iField) * EF_elem(:)
3667-
3668-!\todo activate
3669- first_moments(:, iField) = first_moments(:, iField) &
3670- & + 0.5_rk*chargeTerm(iField)*EF_elem(:)
3671- end do
3672-!write(dbgUnit(1),*) 'spcForceTerm ', spcForceTerm
3673-!write(dbgUnit(1),*) 'electricTerm ', electricTerm
3674-!write(dbgUnit(1),*) 'firstMoments ', first_moments
3675-
3676- matA = 0.0_rk
3677- !build up matrix to solver LSE for actual velocity
3678- do iField = 1, nFields
3679- !set diagonal part
3680- matA(iField, iField) = 1.0_rk
3681- do iField_2 = 1, nFields
3682- matA(iField, iField) = matA(iField, iField) + omega * 0.5_rk &
3683- & * resi_coeff(iField, iField_2) * phi(iField) &
3684- & * moleFrac(iField_2) / paramB
3685- end do
3686- !set non-diagonal part
3687- do iField_2 = 1, nFields
3688- matA(iField, iField_2) = matA(iField, iField_2) - omega * 0.5_rk &
3689- & * resi_coeff(iField, iField_2) * phi(iField_2)&
3690- & * moleFrac(iField) / paramB
3691- end do
3692- end do
3693-
3694- ! invert matrix
3695- invA = invert_matrix( matA )
3696-
3697- ! actual velocity of all species
3698- velocity(1, :) = matmul( invA, first_moments(1,:) ) / mass_dens(:)
3699- velocity(2, :) = matmul( invA, first_moments(2,:) ) / mass_dens(:)
3700- velocity(3, :) = matmul( invA, first_moments(3,:) ) / mass_dens(:)
3701-
3702- ! compute equilibrim velocity
3703- do iField = 1, nFields
3704- eqVel( :, iField ) = velocity( :, iField )
3705- do iField_2 = 1, nFields
3706- eqVel( :, iField ) = eqVel( :, iField ) &
3707- & + resi_coeff( iField, iField_2 ) * phi(iField) &
3708- & * moleFrac(iField_2) &
3709- & * (velocity(:, iField_2) - velocity(:,iField)) &
3710- & / paramB
3711- end do
3712- end do
3713-
3714- !compute mass averaged mixture velocity
3715- velAvg(1) = dot_product( mass_dens, velocity(1,:) )/totmassDens
3716- velAvg(2) = dot_product( mass_dens, velocity(2,:) )/totmassDens
3717- velAvg(3) = dot_product( mass_dens, velocity(3,:) )/totmassDens
3718-
3719- !compute equilibrium and do collision
3720- do iField = 1, nFields
3721- !omega moments = (M^-1 RelaxMat M)*(I+(M^-1 RelMat M)/2)^-1
3722- omegaMoments = scheme%field(iField)%fieldProp%species%mrt(iLevel)%omegaMoments
3723- omegaMomForce = scheme%field(iField)%fieldProp%species%mrt(iLevel)%omegaMomForce
3724-
3725- usqr = dot_product( velAvg, velAvg )
3726- usqr = usqr * t2cs2inv
3727-
3728- do iDir = 1, QQ
3729- vPos = varSys%method%val( stateVarMap(iField) )%state_varPos(iDir)
3730-
3731- ucx = scheme%layout%fStencil%cxDirRK( 1, iDir )*eqVel(1, iField) &
3732- & + scheme%layout%fStencil%cxDirRK( 2, iDir )*eqVel(2, iField) &
3733- & + scheme%layout%fStencil%cxDirRK( 3, iDir )*eqVel(3, iField)
3734-
3735- ucxQuadTerm = dot_product( &
3736- & scheme%layout%fStencil%cxDirRK(:, iDir), velAvg )
3737-
3738- feq = scheme%layout%weight(iDir) * mass_dens(iField) * ( phi(iField) &
3739- & + ucx * cs2inv + ucxQuadTerm * ucxQuadTerm * t2cs4inv - usqr )
3740-
3741- cxcxTerm = 0.0_rk
3742- if ( iDir == scheme%layout%fStencil%restPosition ) then
3743- ! equilibrium at rest
3744- select case( trim(scheme%header%layout) )
3745- case('d2q9')
3746- feq = scheme%layout%weight( iDir ) * mass_dens(iField) * ( &
3747- & ( 9._rk - 5._rk * phi(iField) )/4._rk - usqr )
3748- cxcxTerm(1,1) = scheme%layout%fStencil%cxcx(1,iDir) - cs2
3749- cxcxTerm(2,2) = scheme%layout%fStencil%cxcx(2,iDir) - cs2
3750- cxcxTerm(1,2) = scheme%layout%fStencil%cxcx(3,iDir)
3751- cxcxTerm(2,1) = scheme%layout%fStencil%cxcx(3,iDir)
3752-
3753- case('d3q19')
3754- feq = scheme%layout%weight( iDir ) * mass_dens(iField) * ( &
3755- & ( 3._rk - 2._rk * phi(iField) ) - usqr )
3756-
3757- cxcxTerm(1,1) = scheme%layout%fStencil%cxcx(1,iDir) - cs2
3758- cxcxTerm(2,2) = scheme%layout%fStencil%cxcx(2,iDir) - cs2
3759- cxcxTerm(3,3) = scheme%layout%fStencil%cxcx(3,iDir) - cs2
3760- cxcxTerm(1,2) = scheme%layout%fStencil%cxcx(4,iDir)
3761- cxcxTerm(2,1) = scheme%layout%fStencil%cxcx(4,iDir)
3762- cxcxTerm(2,3) = scheme%layout%fStencil%cxcx(5,iDir)
3763- cxcxTerm(3,2) = scheme%layout%fStencil%cxcx(5,iDir)
3764- cxcxTerm(1,3) = scheme%layout%fStencil%cxcx(6,iDir)
3765- cxcxTerm(3,1) = scheme%layout%fStencil%cxcx(6,iDir)
3766-
3767- end select
3768- end if
3769- fneq(iDir) = ( feq - pdfTmp( vPos ) )
3770-
3771- forceQuadTerm = dot_product( electricTerm(:, iField), &
3772- & matmul( cxcxTerm, velocity(:,iField) ) )
3773-
3774- forceUxTerm = dot_product( scheme%layout%fStencil%cxDirRK( :, iDir ),&
3775- & electricTerm(:, iField) )
3776-
3777- diffForceUxTerm = dot_product( scheme%layout%fStencil%cxDirRK( :, iDir ),&
3778- & diffForce(:, iField) )
3779-
3780- forceTerm(iDir) = scheme%layout%weight(iDir) &
3781- & * ( cs2inv * forceUxTerm + cs4inv * forceQuadTerm )!+ diffForceUxTerm)
3782-!write(dbgUnit(1),*) 'forceQuadTerm ', forceQuadTerm, ' forceUxTerm ',forceUxTerm
3783-!write(dbgUnit(1),*) 'diffForceUxTerm ', diffForceUxTerm
3784- enddo !iDir
3785-
3786-!write(dbgUnit(1),*) 'forceTerm ', forceTerm
3787- fneq_om = matmul( omegaMoments, fneq )
3788-
3789- force_om = matmul( omegaMomForce, forceTerm )
3790-
3791- do iDir = 1, QQ
3792- vPos = varSys%method%val( stateVarMap(iField) )%state_varPos(iDir)
3793- outstate( &
3794-& ?SAVE?( iDir, iField, iElem, QQ, nScalars, nPdfSize,neigh ) ) = &
3795- !& pdfTmp( vPos ) + fneq_om( iDir ) + forceTerm( iDir )
3796- & pdfTmp( vPos ) + fneq_om( iDir ) + force_om( iDir )
3797- enddo !iDir
3798-
3799- enddo !iField
3800- enddo nodeloop
3801-
3802- deallocate(forceTerm)
3803- deallocate(omegaMoments)
3804- deallocate(omegaMomForce)
3805- deallocate(fneq)
3806- deallocate(fneq_om)
3807- deallocate(force_om)
3808-
3809- end subroutine applySrc_electricFieldMSLiquid_2order
3810- ! ************************************************************************ !
3811-
38123465
38133466 ! ************************************************************************ !
38143467 !> Update state with source variable "force" with 2nd order integration
@@ -6117,34 +5770,5 @@
61175770 end function equilFromMacroWTDF
61185771 ! ************************************************************************ !
61195772
6120- ! ************************************************************************** !
6121- !> This function find the corresponding iField of variables in varPos by
6122- !! comparing input_varPos of the variable against its position in derVarPos
6123- pure function findDepField(derVarPos, nFields, input_varPos, nInputs) &
6124- & result(depField)
6125- ! ------------------------------------------------------------------------ !
6126- !> Position of requested derive variable in varSys for all fields
6127- integer, intent(in) :: derVarPos(:)
6128- !> Number of fields
6129- integer, intent(in) :: nFields
6130- !> Position of requested derive variable in varSys for given nInputs
6131- integer, intent(in) :: input_varPos(:)
6132- !> Number of input variable in input_varPos
6133- integer, intent(in) :: nInputs
6134- !> iField of input variable
6135- integer :: depField(nInputs)
6136- ! ------------------------------------------------------------------------ !
6137- integer :: iIn, iField
6138- ! ------------------------------------------------------------------------ !
6139- do iIn = 1, nInputs
6140- do iField = 1, nFields
6141- if ( input_varPos(iIn) == derVarpos(iField) ) &
6142- & depField(iIn) = iField
6143- end do
6144- end do
6145-
6146- end function findDepField
6147- ! ************************************************************************** !
6148-
61495773 end module mus_derQuanMSLiquid_module
61505774 ! **************************************************************************** !
diff -r b6956abee5af -r 35d71a54ed5d source/derived/mus_derQuanPoisson_module.fpp
--- a/source/derived/mus_derQuanPoisson_module.fpp Thu Jun 09 14:16:03 2022 +0200
+++ b/source/derived/mus_derQuanPoisson_module.fpp Tue Aug 30 12:25:45 2022 +0200
@@ -99,6 +99,8 @@
9999 public :: mus_append_derVar_poisson
100100 public :: deriveAuxPoisson_fromState
101101 public :: deriveEquilPoisson_fromAux
102+ public :: derivePotential_forElement
103+ public :: derivePotential_fromIndex
102104
103105 public :: applySrc_chargeDensity_2ndOrd
104106 public :: applySrc_chargeDensity_1stOrd
diff -r b6956abee5af -r 35d71a54ed5d source/derived/mus_derQuan_module.fpp
--- a/source/derived/mus_derQuan_module.fpp Thu Jun 09 14:16:03 2022 +0200
+++ b/source/derived/mus_derQuan_module.fpp Tue Aug 30 12:25:45 2022 +0200
@@ -121,10 +121,15 @@
121121 private
122122
123123 public :: mus_append_derVar_fluid
124+ public :: mus_deriveVelocity
125+ public :: mus_deriveKE
126+ public :: mus_deriveMomentum
127+ public :: mus_derivePressure
124128
125129 public :: derivePressure, derivePressure_fromIndex
126130 public :: deriveDensity
127131 public :: deriveDensity_fromIndex
132+ public :: deriveVelocity
128133 public :: deriveShearStress
129134 public :: deriveShearMag
130135 public :: deriveWSS3D
@@ -132,6 +137,7 @@
132137 public :: deriveTemp
133138 public :: deriveShearRate
134139 public :: deriveBndForce
140+ public :: deriveMomentumChange
135141 ! equilbrium from macro uses different interface defined in
136142 ! mus_variable_module
137143 public :: deriveEquil_FromMacro
@@ -3655,70 +3661,6 @@
36553661 ! ****************************************************************************** !
36563662
36573663
3658-! ****************************************************************************** !
3659- !> Initiates the calculation of Velocity.
3660- !! This routine sets the function Pointer for velocity calcualtion and calls
3661- !! the generice get Value of Index routine
3662- !!
3663- !! The interface has to comply to the abstract interface
3664- !! [[tem_varSys_module:tem_varSys_proc_getValOfIndex]].
3665- !!
3666- recursive subroutine deriveVelocity_fromIndex( fun, varSys, time, iLevel, &
3667- & idx, idxLen, nVals, res )
3668- !> Description of the method to obtain the variables, here some preset
3669- !! values might be stored, like the space time function to use or the
3670- !! required variables.
3671- class(tem_varSys_op_type), intent(in) :: fun
3672-
3673- !> The variable system to obtain the variable from.
3674- type(tem_varSys_type), intent(in) :: varSys
3675-
3676- !> Point in time at which to evaluate the variable.
3677- type(tem_time_type), intent(in) :: time
3678-
3679- !> Level on which values are requested
3680- integer, intent(in) :: iLevel
3681-
3682- !> Index of points in the growing array and variable val array to
3683- !! return.
3684- !! Size: most times nVals, if contiguous arrays are used it depends
3685- !! on the number of first indices
3686- integer, intent(in) :: idx(:)
3687-
3688- !> With idx as start index in contiguous memory,
3689- !! idxLength defines length of each contiguous memory
3690- !! Size: dependes on number of first index for contiguous array,
3691- !! but the sum of all idxLen is equal to nVals
3692- integer, optional, intent(in) :: idxLen(:)
3693-
3694- !> Number of values to obtain for this variable (vectorized access).
3695- integer, intent(in) :: nVals
3696-
3697- !> Resulting values for the requested variable.
3698- !!
3699- !! Dimension: n requested entries x nComponents of this variable
3700- !! Access: (iElem-1)*fun%nComponents + iComp
3701- real(kind=rk), intent(out) :: res(:)
3702- ! ---------------------------------------------------------------------------
3703- !> Function pointer to perform specific operation.
3704- procedure(mus_derive_fromPDF), pointer :: fnCalcPtr
3705- ! ---------------------------------------------------------------------------
3706-
3707- fnCalcPtr => mus_derivevelocity
3708-
3709- call mus_generic_varFromPDF_fromIndex( &
3710- & fun = fun, &
3711- & varSys = varSys, &
3712- & time = time, &
3713- & iLevel = iLevel, &
3714- & idx = idx, &
3715- & nVals = nVals, &
3716- & fnCalcPtr = fnCalcPtr, &
3717- & res = res )
3718-
3719- end subroutine deriveVelocity_fromIndex
3720-! ****************************************************************************** !
3721-
37223664
37233665 ! ****************************************************************************** !
37243666 !> Initiates the calculation of pressure.
diff -r b6956abee5af -r 35d71a54ed5d source/derived/mus_derVarPos_module.f90
--- a/source/derived/mus_derVarPos_module.f90 Thu Jun 09 14:16:03 2022 +0200
+++ b/source/derived/mus_derVarPos_module.f90 Tue Aug 30 12:25:45 2022 +0200
@@ -43,6 +43,8 @@
4343 public :: mus_derive_FromMacro_dummy
4444 public :: mus_derive_FromState_dummy
4545 public :: mus_derive_FromPreColState_dummy
46+ public :: derive_equilFromAux_dummy
47+ public :: derive_auxfromstate_dummy
4648
4749 !> This type stores the position of each variable in the global sys
4850 type mus_derVarPos_type
diff -r b6956abee5af -r 35d71a54ed5d source/derived/mus_transport_var_module.f90
--- a/source/derived/mus_transport_var_module.f90 Thu Jun 09 14:16:03 2022 +0200
+++ b/source/derived/mus_transport_var_module.f90 Tue Aug 30 12:25:45 2022 +0200
@@ -70,6 +70,7 @@
7070 public :: mus_create_poss_transVar
7171 public :: mus_load_transport_var
7272 public :: mus_init_transport_var
73+ public :: get_transport_params
7374
7475 ! ***************************************************************************!
7576 !> Description contains index to access value using variable function
diff -r b6956abee5af -r 35d71a54ed5d source/intp/mus_interpolate_compact_module.fpp
--- a/source/intp/mus_interpolate_compact_module.fpp Thu Jun 09 14:16:03 2022 +0200
+++ b/source/intp/mus_interpolate_compact_module.fpp Tue Aug 30 12:25:45 2022 +0200
@@ -72,7 +72,6 @@
7272
7373 private
7474
75- public :: unitTest_interpolate_quadraticCompact
7675 public :: fillMyGhostsFromFiner_quadraticCompact2d
7776 public :: fillFinerGhostsFromMe_quadraticCompact2d
7877 public :: fillMyGhostsFromFiner_quadraticCompact3d
@@ -1119,109 +1118,6 @@
11191118
11201119
11211120 ! ****************************************************************************** !
1122- !>
1123- !!
1124- subroutine unitTest_interpolate_quadraticCompact( iCase, passed )
1125- ! ---------------------------------------------------------------------------
1126- !>
1127- integer, intent(in) :: iCase
1128- !>
1129- logical, intent(out) :: passed
1130- ! ---------------------------------------------------------------------------
1131- ! non-equilibrium part of the resulting pdf
1132- real(kind=rk) :: sourceMom( 9, 4 )
1133- real(kind=rk) :: targetMom( 9 ) ! non-equilibrium part of the resulting pdf
1134- real(kind=rk) :: refX(2) ! reference values
1135- real(kind=rk) :: error(4) ! error values
1136- integer :: nSourceElems, SE, SW, NE, NW
1137- ! ---------------------------------------------------------------------------
1138- passed = .false.
1139- error = 0._rk
1140- nSourceElems = 4
1141- write(logUnit(1),*)' Performing unit test for interpolation moment '// &
1142- & 'gradient routine'
1143- write(logUnit(1),*)' Current case for child ', iCase
1144-
1145- select case( iCase )
1146- case( 1, 5)
1147- SW = 4
1148- SE = 3
1149- NW = 2
1150- NE = 1
1151- case( 2, 6)
1152- SW = 2
1153- SE = 4
1154- NW = 1
1155- NE = 3
1156- case( 3, 7)
1157- SW = 3
1158- SE = 1
1159- NW = 4
1160- NE = 2
1161- case( 4, 8)
1162- SW = 1
1163- SE = 2
1164- NW = 3
1165- NE = 4
1166- end select
1167-
1168- ! sourceMom(9,4)
1169- sourceMom = 0._rk
1170- sourceMom( 1, : ) = 1._rk
1171- sourceMom( 3, : ) = 0.1_rk
1172-?? if( PUSH ) then
1173- sourceMom( 6, SW ) = -0.074074074074074_rk ! SW
1174- sourceMom( 6, SE ) = -0.037037037037037_rk ! SE
1175- sourceMom( 6, NW ) = 0.074074074074074_rk ! NW
1176- sourceMom( 6, NE ) = 0.037037037037037_rk ! NE
1177-?? else
1178- sourceMom( 6, SW ) = 0.059259259259259_rk ! SW
1179- sourceMom( 6, SE ) = 0.029629629629630_rk ! SE
1180- sourceMom( 6, NW ) = -0.059259259259259_rk ! NW
1181- sourceMom( 6, NE ) = -0.029629629629630_rk ! NE
1182-?? endif
1183- refX(1) = 0.05625_rk
1184- refX(2) = 0.1_rk
1185- ! targetMom(1:2) = mus_interpolate_quadraticCompact2d_leastSq( &
1186- ! & sourceMom(2:3, 1:nSourceElems), &
1187- ! & (/0.25_rk, 0.25_rk/), sourceMom(4:6, 1:nSourceElems), &
1188- ! & nSourceElems, matrix )
1189- error(1:2) = (refX(1:2)-targetMom(1:2))/refX(1:2)*100._rk
1190-?? if( PUSH ) then
1191-?? else
1192- write(logUnit(6),*)'X reference result ', refX(1), ' intp ', targetMom(1),&
1193- & ' Error [%] ', error(1)
1194- write(logUnit(6),*)'Y reference result ', refX(2), ' intp ', targetMom(2),&
1195- & ' Error [%] ', error(2)
1196-?? endif
1197-
1198- refX(1) = 0.075_rk
1199- refX(2) = 0.1_rk
1200- ! targetMom(1:2) = mus_interpolate_quadraticCompact2d_leastSq( &
1201- ! & sourceMom(2:3, 1:nSourceElems), &
1202- ! & (/0.5_rk, 0.5_rk/), sourceMom(4:6, 1:nSourceElems), &
1203- ! & nSourceElems, matrix )
1204- error(3:4) = (refX(1:2)-targetMom(1:2))/refX(1:2)*100._rk
1205-
1206-?? if( PUSH ) then
1207-?? else
1208- write(logUnit(6),*)'X reference result ', refX(1), ' intp ', targetMom(1),&
1209- & ' Error [%] ', error(3)
1210- write(logUnit(6),*)'Y reference result ', refX(2), ' intp ', targetMom(2),&
1211- & ' Error [%] ', error(4)
1212-?? endif
1213-
1214- if( abs( maxval( error ) ) > 0.00001_rk ) then
1215- passed = .false.
1216- else
1217- passed = .true.
1218- endif
1219-
1220- end subroutine unitTest_interpolate_quadraticCompact
1221-! ****************************************************************************** !
1222-
1223-
1224-! ****************************************************************************** !
12251121 !> Calculate the equilibrium distribution function in all directions
12261122 !!
12271123 !! The equilibrium distribution function is:\n
diff -r b6956abee5af -r 35d71a54ed5d source/intp/mus_interpolate_d3q19_module.fpp
--- a/source/intp/mus_interpolate_d3q19_module.fpp Thu Jun 09 14:16:03 2022 +0200
+++ b/source/intp/mus_interpolate_d3q19_module.fpp Tue Aug 30 12:25:45 2022 +0200
@@ -454,69 +454,6 @@
454454 end subroutine fillCoarser_compactd3q19
455455 ! **************************************************************************** !
456456
457-! **************************************************************************** !
458- pure function getStrain( f, rho, vel, omega, cxDir ) result( S )
459- ! -------------------------------------------------------------------- !
460- !> stencil layout
461- integer, intent(in) :: cxDir(3,QQ)
462- !> pdf array ( post-collision value )
463- real(kind=rk), intent(in) :: f(QQ), rho, vel(3)
464- !> relaxation parameter
465- real(kind=rk), intent(in) :: omega
466- !> output array: strain rate tensor
467- real(kind=rk) :: S(3,3)
468- ! -------------------------------------------------------------------- !
469- real(kind=rk) :: fac
470- ! -------------------------------------------------------------------- !
471-
472- S(1,1) = - f(1) - f(4) - sum(f(11:18)) &
473- & + vel(1) * vel(1) + cs2 * rho
474- S(2,2) = - f(2) - f(5) - sum(f(7:10)) - sum(f(15:18)) &
475- & + vel(2) * vel(2) + cs2 * rho
476- S(3,3) = - f(3) - f(6) - sum(f(7:10)) - sum(f(11:14)) &
477- & + vel(3) * vel(3) + cs2 * rho
478-
479- S(2,1) = - f(15) + f(16) + f(17) - f(18) + vel(2) * vel(1)
480- S(3,2) = - f( 7) + f( 8) + f( 9) - f(10) + vel(3) * vel(2)
481- S(3,1) = - f(11) + f(12) + f(13) - f(14) + vel(3) * vel(1)
482- S(1,2) = S(2,1)
483- S(2,3) = S(3,2)
484- S(1,3) = S(3,1)
485-
486- ! Multiply factor and convert to pre-collision value
487- fac = omega * cs2inv / ?post2pre?( omega ) * div1_2
488- s = s * fac
489-
490- end function getStrain
491-! **************************************************************************** !
492-
493-! **************************************************************************** !
494- pure function getNEq( S, omega, cxDirRK, weight ) result( nEq )
495- ! -------------------------------------------------------------------- !
496- !> Strain rate tensor. It is a symmetric 3x3 matrix
497- real(kind=rk), intent(in) :: S(3,3)
498- real(kind=rk), intent(in) :: cxDirRK(3,QQ)
499- real(kind=rk), intent(in) :: weight(QQ)
500- real(kind=rk), intent(in) :: omega
501- real(kind=rk) :: nEq(QQ)
502- ! -------------------------------------------------------------------- !
503- integer :: iDir
504- real(kind=rk) :: fac
505- ! -------------------------------------------------------------------- !
506-
507- do iDir = 1, QQ
508- nEq(iDir) = s(1,1) * ( cxDirRK(1,iDir) * cxDirRK(1,iDir) - cs2 ) &
509- & + s(2,2) * ( cxDirRK(2,iDir) * cxDirRK(2,iDir) - cs2 ) &
510- & + s(3,3) * ( cxDirRK(3,iDir) * cxDirRK(3,iDir) - cs2 ) &
511- & + s(1,2) * cxDirRK(1,iDir) * cxDirRK(2,iDir) * 2.0_rk &
512- & + s(1,3) * cxDirRK(1,iDir) * cxDirRK(3,iDir) * 2.0_rk &
513- & + s(2,3) * cxDirRK(2,iDir) * cxDirRK(3,iDir) * 2.0_rk
514- end do
515- fac = cs2inv / omega * ?pre2post?( omega )
516- nEq(:) = - weight(:) * nEq(:) * fac
517-
518- end function getNEq
519-! **************************************************************************** !
520457
521458 ! **************************************************************************** !
522459 subroutine calc_moments( state, nScalars, nSize, elementList, nElems, &
diff -r b6956abee5af -r 35d71a54ed5d source/intp/mus_interpolate_debug_module.fpp
--- a/source/intp/mus_interpolate_debug_module.fpp Thu Jun 09 14:16:03 2022 +0200
+++ b/source/intp/mus_interpolate_debug_module.fpp Tue Aug 30 12:25:45 2022 +0200
@@ -90,28 +90,28 @@
9090 public :: do_nothing
9191 public :: do_nothing_arbi
9292
93-?? if ( .not. SOA) then
94- ! public :: fillMyGhostsFromFiner_average
95- ! public :: fillFinerGhostsFromMe_average
93+?? if (.not. SOA) then
94+ public :: fillMyGhostsFromFiner_average
95+ public :: fillFinerGhostsFromMe_average
9696
97- ! public :: fillMyGhostsFromFiner_debug
98- ! public :: fillFinerGhostsFromMe_debug
97+ public :: fillMyGhostsFromFiner_debug
98+ public :: fillFinerGhostsFromMe_debug
9999
100- ! public :: fillMyGhostsFromFiner_copyFirst
101- ! public :: fillFinerGhostsFromMe_copyFirst
100+ public :: fillMyGhostsFromFiner_copyFirst
101+ public :: fillFinerGhostsFromMe_copyFirst
102102
103103 public :: fillMyGhostsFromFiner_TGV2D
104104 public :: fillFinerGhostsFromMe_TGV2D
105105
106- ! public :: fillMyGhostsFromFiner_TGV3D
107- ! public :: fillFinerGhostsFromMe_TGV3D
106+ public :: fillMyGhostsFromFiner_TGV3D
107+ public :: fillFinerGhostsFromMe_TGV3D
108108 ?? endif
109109
110110 public :: TGV_2D
111111
112112 real(kind=rk), parameter :: debugValue = 1.0d0
113113
114- contains
114+contains
115115
116116 ! ****************************************************************************** !
117117 !> Fill GhostFromFiner elements on my level with debug value
diff -r b6956abee5af -r 35d71a54ed5d source/mus_construction_module.fpp
--- a/source/mus_construction_module.fpp Thu Jun 09 14:16:03 2022 +0200
+++ b/source/mus_construction_module.fpp Tue Aug 30 12:25:45 2022 +0200
@@ -166,11 +166,11 @@
166166 ! !!
167167 ! !! # Additional Tasks
168168 ! !!
169- ! !! - receive [[tem_construction_module:tem_buildhorizontaldependencies]]
169+ ! !! - receive [[tem_construction_module:tem_build_horizontaldependencies]]
170170 ! !! "horizontal"
171171 ! !! (within a level for the element updates)
172- ! !! - and [[tem_construction:tem_buildverticaldependencies]] "vertical"
173- ! !! dependencies (between levels for ghost-interpolations).
172+ ! !! - and [[tem_construction_module:tem_build_verticaldependencies]]
173+ ! !! "vertical" dependencies (between levels for ghost-interpolations).
174174 ! !! - The main state vector and the neighbor lists on which the kernel then
175175 ! !! acts is created
176176 ! !! - The MPI buffers are created.
diff -r b6956abee5af -r 35d71a54ed5d source/mus_debug_tools_module.fpp
--- a/source/mus_debug_tools_module.fpp Thu Jun 09 14:16:03 2022 +0200
+++ b/source/mus_debug_tools_module.fpp Tue Aug 30 12:25:45 2022 +0200
@@ -48,6 +48,8 @@
4848 private
4949
5050 public :: dumpPdfAll
51+ public :: dumpPdfSingleElem
52+ public :: dumpPdfSingleElemSAVE
5153
5254 contains
5355
diff -r b6956abee5af -r 35d71a54ed5d source/mus_derivedQuantities_module.fpp
--- a/source/mus_derivedQuantities_module.fpp Thu Jun 09 14:16:03 2022 +0200
+++ b/source/mus_derivedQuantities_module.fpp Tue Aug 30 12:25:45 2022 +0200
@@ -257,70 +257,44 @@
257257 ! ****************************************************************************** !
258258
259259
260-! ****************************************************************************** !
261- !> Calculate the velocity in all 3 directions
262- !!
263- pure function getVelocity_forElemFromState_Force( state, elem, force, &
264- & dtLB, stencil, varPos, nScalars ) result( vel )
265- ! ---------------------------------------------------------------------------
266- real(kind=rk), intent(in) :: state(:) !< state array of a level
267- type(tem_stencilHeader_type), intent(in) :: stencil
268- real(kind=rk), intent(in) :: force(3)
269- real(kind=rk), intent(in) :: dtLB
270- integer, intent(in) :: elem
271- integer, intent(in) :: varPos(:) !< varPos of current field variable
272- integer, intent(in) :: nScalars !< number of scalars in global system
273- real(kind=rk) :: vel(3) !< return value
274- ! ---------------------------------------------------------------------------
275- real(kind=rk) :: dens
276- integer :: iDir
277- ! ---------------------------------------------------------------------------
278- integer :: nElems
279- ! ---------------------------------------------------------------------------
280- nElems = size( state ) / nScalars
281-
282- vel = 0._rk
283- dens = 0._rk
284- do iDir = 1,stencil%QQ
285- vel(:) = vel(:) + state( ?IDX?( varPos(iDir), elem, nScalars, nElems)) &
286- & * stencil%cxDirRK( :, iDir ) &
287- ! Forcing Term
288- & + dtLB*0.5_rk*force(:)
289-
290- dens = dens + state( ?IDX?( varPos(iDir), elem, nScalars, nElems))
291- enddo
292- vel = vel / dens
293-
294- end function getVelocity_forElemFromState_Force
295-! ****************************************************************************** !
260+!!! ****************************************************************************** !
261+!! !> Calculate the velocity in all 3 directions
262+!! !!
263+!! pure function getVelocity_forElemFromState_Force( state, elem, force, &
264+!! & dtLB, stencil, varPos, nScalars ) result( vel )
265+!! ! ---------------------------------------------------------------------------
266+!! real(kind=rk), intent(in) :: state(:) !< state array of a level
267+!! type(tem_stencilHeader_type), intent(in) :: stencil
268+!! real(kind=rk), intent(in) :: force(3)
269+!! real(kind=rk), intent(in) :: dtLB
270+!! integer, intent(in) :: elem
271+!! integer, intent(in) :: varPos(:) !< varPos of current field variable
272+!! integer, intent(in) :: nScalars !< number of scalars in global system
273+!! real(kind=rk) :: vel(3) !< return value
274+!! ! ---------------------------------------------------------------------------
275+!! real(kind=rk) :: dens
276+!! integer :: iDir
277+!! ! ---------------------------------------------------------------------------
278+!! integer :: nElems
279+!! ! ---------------------------------------------------------------------------
280+!! nElems = size( state ) / nScalars
281+!!
282+!! vel = 0._rk
283+!! dens = 0._rk
284+!! do iDir = 1,stencil%QQ
285+!! vel(:) = vel(:) + state( ?IDX?( varPos(iDir), elem, nScalars, nElems)) &
286+!! & * stencil%cxDirRK( :, iDir ) &
287+!! ! Forcing Term
288+!! & + dtLB*0.5_rk*force(:)
289+!!
290+!! dens = dens + state( ?IDX?( varPos(iDir), elem, nScalars, nElems))
291+!! enddo
292+!! vel = vel / dens
293+!!
294+!! end function getVelocity_forElemFromState_Force
295+!!! ****************************************************************************** !
296296
297297
298-! ****************************************************************************** !
299- !> Calculate the momentum in all 3 directions for a subset of pdfs
300- !! the subset is provided and must be of the size of DOFs per element
301- !!
302- pure function getMomentum_forPdfSubset_force( subset, force, dtLB, &
303- & stencil, varPos ) result( mom )
304- ! ---------------------------------------------------------------------------
305- real(kind=rk), intent(in) :: force(3)
306- type(tem_stencilHeader_type), intent(in) :: stencil
307- real(kind=rk), intent(in) :: dtLB
308- real(kind=rk), intent(in) :: subset(:)
309- integer, intent(in) :: varPos(:) !< varPos of current field variable
310- real(kind=rk) :: mom(3)!< return value
311- ! ---------------------------------------------------------------------------
312- integer :: iDir
313- ! ---------------------------------------------------------------------------
314-
315- mom = 0._rk
316- do iDir = 1,stencil%QQ
317- mom(:) = mom(:) + subset(varPos(iDir)) * real(stencil%cxDir(:,iDir), rk) &
318- ! Forcing Term
319- & + dtLB*div1_2*force(:)
320- enddo
321-
322- end function getMomentum_forPdfSubset_force
323-! ****************************************************************************** !
324298
325299
326300 ! ****************************************************************************** !
@@ -396,34 +370,6 @@
396370 ! ****************************************************************************** !
397371
398372
399-! ****************************************************************************** !
400- !> Incompressible formulation of the equilibrium for specific direction.
401- !!
402- !! @warning This function is deprecated. Use the macro EqDir_incomp instead.
403- pure function getEqByDensVelPerDir_incomp( dens, vel, cx, weight ) &
404- & result( equil )
405- ! --------------------------------------------------------------------------
406- real(kind=rk), intent(in) :: dens
407- real(kind=rk), intent(in) :: vel(3)
408- real(kind=rk), intent(in) :: cx(3)
409- real(kind=rk), intent(in) :: weight
410- real(kind=rk) :: equil !< return value
411- ! --------------------------------------------------------------------------
412- real(kind=rk) :: ucx, usq
413- ! --------------------------------------------------------------------------
414- ! square of velocity
415- usq = vel(1)*vel(1) + vel(2)*vel(2) + vel(3)*vel(3)
416-
417- ! velocity times lattice unit velocity
418- ucx = cx(1) * vel(1) + cx(2) * vel(2) + cx(3) * vel(3)
419-
420- ! calculate equilibrium density
421- equil = weight * ( dens + rho0*(ucx*cs2inv &
422- & + ucx*ucx*cs2inv*cs2inv*div1_2 &
423- & - usq*cs2inv*div1_2 ) )
424-
425- end function getEqByDensVelPerDir_incomp
426-! ****************************************************************************** !
427373
428374 ! ****************************************************************************** !
429375 !> Calculate the equilibrium distribution function in all directions
@@ -661,65 +607,65 @@
661607 ! ****************************************************************************** !
662608
663609
664-! ****************************************************************************** !
665- !> Calculate the shear rate tensor (strain rate) by acoustic scaling.
666- !! @todo: This function actually uses diffusive scaling to calculate shear
667- !! rate.
668- !!
669- pure function getShearRateTensor_acoustic_forPdfSubset( subset, omega, &
670- & layout, incompressible ) result( S )
671- ! ---------------------------------------------------------------------------
672- type(mus_scheme_layout_type), intent(in) :: layout !scheme layout
673- real(kind=rk), intent(in) :: subset(:) !< pdf array
674- real(kind=rk), intent(in) :: omega
675- real(kind=rk) :: S(3,3)
676- logical, intent(in), optional :: incompressible
677- ! ---------------------------------------------------------------------------
678- real(kind=rk) :: vel(3) ! velocity components
679- real(kind=rk) :: diagVal(3,3), rho, rho0
680- real(kind=rk) :: fTmp(layout%fStencil%QQ), rhoD
681- integer :: iDir, iVal, jVal, order(3)
682- logical :: incomp
683- ! ---------------------------------------------------------------------------
684- incomp = .false.
685- if ( present( incompressible )) incomp = incompressible
686-
687- Rho = 0._rk
688- Vel = 0._rk
689- do iDir = 1, layout%fStencil%QQ
690- fTmp( iDir ) = subset( iDir )
691- vel(:) = vel(:) + fTmp( iDir ) * dble(layout%fStencil%cxDir(:,iDir))
692- end do
693- rho = sum( fTmp(:) )
694-
695- if ( .not. incomp ) then
696- vel = vel / rho
697- rhoD = rho
698- else
699- rhoD = rho0
700- end if
701-
702- ! Get the second velocity moments of the source element's pdf
703- diagVal(:,:) = 0._rk
704- diagVal(1,1) = cs2*Rho
705- diagVal(2,2) = cs2*Rho
706- diagVal(3,3) = cs2*Rho
707- do jVal = 1, layout%fStencil%nDims
708- do iVal = 1, layout%fStencil%nDims
709- order = 0
710- order(iVal) = order(iVal) +1
711- order(jVal) = order(jVal) +1
712- S(iVal,jVal) = get_moment( layout%fStencil%QQ, layout%fStencil%cxDir, order, fTmp)
713-
714- ! Now compute the S(1) from the second velocity moment
715- S(iVal,jVal) = rhoD * Vel(iVal) * Vel(jVal) + diagVal(iVal,jVal) - S(iVal,jVal)
716- end do
717- end do
718-
719- S(:,:) = S(:,:) * omega * cs2inv * convPrePost(omega) * div1_2
720-
721- end function getShearRateTensor_acoustic_forPdfSubset
722-! ****************************************************************************** !
610+!!! ****************************************************************************** !
611+!! !> Calculate the shear rate tensor (strain rate) by acoustic scaling.
612+!! !! @todo: This function actually uses diffusive scaling to calculate shear
613+!! !! rate.
614+!! !!
615+!! pure function getShearRateTensor_acoustic_forPdfSubset( subset, omega, &
616+!! & layout, incompressible ) result( S )
617+!! ! ---------------------------------------------------------------------------
618+!! type(mus_scheme_layout_type), intent(in) :: layout !scheme layout
619+!! real(kind=rk), intent(in) :: subset(:) !< pdf array
620+!! real(kind=rk), intent(in) :: omega
621+!! real(kind=rk) :: S(3,3)
622+!! logical, intent(in), optional :: incompressible
623+!! ! ---------------------------------------------------------------------------
624+!! real(kind=rk) :: vel(3) ! velocity components
625+!! real(kind=rk) :: diagVal(3,3), rho, rho0
626+!! real(kind=rk) :: fTmp(layout%fStencil%QQ), rhoD
627+!! integer :: iDir, iVal, jVal, order(3)
628+!! logical :: incomp
629+!! ! ---------------------------------------------------------------------------
630+!! incomp = .false.
631+!! if ( present( incompressible )) incomp = incompressible
632+!!
633+!! Rho = 0._rk
634+!! Vel = 0._rk
635+!! do iDir = 1, layout%fStencil%QQ
636+!! fTmp( iDir ) = subset( iDir )
637+!! vel(:) = vel(:) + fTmp( iDir ) * dble(layout%fStencil%cxDir(:,iDir))
638+!! end do
639+!! rho = sum( fTmp(:) )
640+!!
641+!! if ( .not. incomp ) then
642+!! vel = vel / rho
643+!! rhoD = rho
644+!! else
645+!! rhoD = rho0
646+!! end if
647+!!
648+!! ! Get the second velocity moments of the source element's pdf
649+!! diagVal(:,:) = 0._rk
650+!! diagVal(1,1) = cs2*Rho
651+!! diagVal(2,2) = cs2*Rho
652+!! diagVal(3,3) = cs2*Rho
653+!! do jVal = 1, layout%fStencil%nDims
654+!! do iVal = 1, layout%fStencil%nDims
655+!! order = 0
656+!! order(iVal) = order(iVal) +1
657+!! order(jVal) = order(jVal) +1
658+!! S(iVal,jVal) = get_moment( layout%fStencil%QQ, layout%fStencil%cxDir, order, fTmp)
659+!!
660+!! ! Now compute the S(1) from the second velocity moment
661+!! S(iVal,jVal) = rhoD * Vel(iVal) * Vel(jVal) + diagVal(iVal,jVal) - S(iVal,jVal)
662+!! end do
663+!! end do
664+!!
665+!! S(:,:) = S(:,:) * omega * cs2inv * convPrePost(omega) * div1_2
666+!!
667+!! end function getShearRateTensor_acoustic_forPdfSubset
668+!!! ****************************************************************************** !
723669
724670
725671 ! ****************************************************************************** !
@@ -908,69 +854,6 @@
908854 ! ****************************************************************************** !
909855
910856
911-! ****************************************************************************** !
912- !> author: Jiaxing Qi
913- !! Calculate the viscous shear stress (exclude pressure)
914- !!
915- !! The shear stress formula is:\n
916- !! \[
917- !! \tau_{\alpha \beta}=
918- !! -(1-\frac{\omega}{2}) \sum_{i} f^{neq}_{i} c_{i\alpha} c_{i\beta}
919- !! \]\n
920- !! where \( \tau_{\alpha \beta} \) is the stress
921- !! in the \(\beta\)-direction on a face normal to the \(\alpha\)-axis,\n
922- !! \( f^{neq}_i = f_i - f^{eq}_i \) is the non-equilibirium density.\n
923- !! For more information, please refer to:\n
924- !! Krueger T, Varnik F, Raabe D. Shear stress in lattice Boltzmann
925- !! simulations. Physical Review E. 2009;79(4):1-14.
926- !!
927- !! Notice:\n
928- !! Normally shear stress should be calculated by PRE-collision pdf.\n
929- !! But this routine utilize POST-collision pdf (SAVE)!
930- !!
931- pure function getShearStressPost( state, neigh, elem, omega, layout, iField, &
932- & nScalars ) result( tau )
933- ! ---------------------------------------------------------------------------
934- type(mus_scheme_layout_type), intent(in) :: layout !scheme layout
935- real(kind=rk), intent(in) :: state(:) !< pdf array
936- real(kind=rk), intent(in) :: omega
937- integer, intent(in) :: neigh(:) !< connectivity vector
938- integer, intent(in) :: elem !< treeID of target element
939- integer, intent(in) :: iField !< current field
940- integer, intent(in) :: nScalars !< number of scalars in global system
941- real(kind=rk) :: tau(6) !< Shear stress includes xx, yy, zz, xy, yz, xz
942- ! ---------------------------------------------------------------------------
943- ! local variables
944- real(kind=rk) :: equil(layout%fStencil%QQ), nonEq(layout%fStencil%QQ)
945- real(kind=rk) :: rho, vel(3)
946- real(kind=rk) :: fTmp(layout%fStencil%QQ) !< pdf array of current element
947- integer :: iDir
948- ! ---------------------------------------------------------------------------
949- integer :: nElems, QQ
950- ! ---------------------------------------------------------------------------
951- QQ = layout%fStencil%QQ
952- nElems = size( neigh ) / QQ
953-
954- rho = 0._rk
955- vel = 0._rk
956- do iDir = 1, QQ
957- ! use POST-collision pdf
958- fTmp(iDir) = state( ?SAVE?( iDir, iField, elem, QQ, nScalars, nElems, neigh ) )
959- ! compute density and velocity
960- rho = rho + fTmp(iDir)
961- vel(:) = vel(:) + fTmp(iDir) * layout%fStencil%cxDirRK(:,iDir)
962- end do
963-
964- ! compute Eq and nonEq
965- equil(:) = getEquilibrium( rho, vel, layout )
966- nonEq = equil - fTmp
967-
968- ! compute shear stress
969- tau(:) = secondMom( layout%fStencil%cxcx(:,:), nonEq(:), QQ )
970- tau(:) = ( omega * div1_2 - 1._rk ) * tau(:)
971-
972- end function getShearStressPost
973-! ****************************************************************************** !
974857
975858 ! ****************************************************************************** !
976859 !> author: Jiaxing Qi
@@ -1064,36 +947,6 @@
1064947 end function getWSS
1065948 ! ****************************************************************************** !
1066949
1067-! ****************************************************************************** !
1068- !> author: Jiaxing Qi
1069- !! Calculate the wall shear stress (WSS)
1070- !!
1071- !! @todo: add formule and reference!
1072- function getWSS_vonMises( state, neigh, elem, omega, layout, iField, &
1073- & nScalars ) result( wss )
1074- ! ---------------------------------------------------------------------------
1075- type(mus_scheme_layout_type), intent(in) :: layout !scheme layout
1076- real(kind=rk), intent(in) :: state(:)
1077- integer, intent(in) :: neigh(:) !< connectivity vector
1078- real(kind=rk), intent(in) :: omega
1079- integer, intent(in) :: elem
1080- integer, intent(in) :: iField !< current field
1081- integer, intent(in) :: nScalars !< number of scalars in global system
1082- real(kind=rk) :: wss !< return value
1083- ! ---------------------------------------------------------------------------
1084- ! local variables
1085- real(kind=rk) :: tau(6) ! shear stress components
1086- ! ---------------------------------------------------------------------------
1087- wss = 0._rk
1088- tau(:) = getShearstressPost( state, neigh, elem, omega, layout, iField, &
1089- & nScalars )
1090- ! Check wikipedia article on von Mises stress criterion
1091- wss = (tau(1)-tau(2))**2 + (tau(2)-tau(3))**2 + (tau(1)-tau(3))**2 &
1092- & + 6._rk*( tau(4)*tau(4) + tau(5)*tau(5) + tau(6)*tau(6))
1093- wss = sqrt( wss * div1_2 )
1094-
1095- end function getWSS_vonMises
1096-! ****************************************************************************** !
1097950
1098951
1099952 ! ****************************************************************************** !
diff -r b6956abee5af -r 35d71a54ed5d source/mus_flow_module.fpp
--- a/source/mus_flow_module.fpp Thu Jun 09 14:16:03 2022 +0200
+++ b/source/mus_flow_module.fpp Tue Aug 30 12:25:45 2022 +0200
@@ -1357,77 +1357,6 @@
13571357 ! ****************************************************************************** !
13581358
13591359
1360- ! **************************************************************************** !
1361- !> This routine initialize omega values in state vector for nonNewtonian fluid
1362- !!
1363- subroutine mus_init_omega( state, tree, nElems, nSize, varSys, omegaPos, &
1364- & total, ic )
1365- ! ---------------------------------------------------------------------------
1366- real(kind=rk), intent(inout) :: state(:) !< state vector
1367- !> treelmesh
1368- type( treelmesh_type ), intent(in) :: tree
1369- !> initial condition
1370- type( tem_ini_condition_type ) :: ic
1371- !> Number of Elements
1372- integer, intent(in) :: nElems
1373- !> Number of Elements as Size
1374- integer, intent(in) :: nSize
1375- !> variable system
1376- type( tem_varSys_type ), intent(in) :: varSys
1377- !> Position of omega variable in varSys
1378- integer, intent(in) :: omegaPos
1379- !> total list of elements
1380- integer(kind=long_k), intent(in) :: total(:)
1381- ! ---------------------------------------------------------------------------
1382- real(kind=rk), allocatable :: omegas(:) ! omega values from ic
1383- real(kind=rk), allocatable :: coord(:,:) ! coordinates
1384- integer :: iElem, iChunk, nChunks, chunkSize, nChunkElems, elemPos, elemoff
1385- integer :: nScalars, posInState
1386- ! ---------------------------------------------------------------------------
1387-
1388- write(logUnit(1),"(A)") "Initialize omega by omegas from initial condition."
1389-
1390- ! iElem = nSize
1391- nScalars = varSys%nScalars
1392- posInState = varSys%method%val(omegaPos)%state_varPos(1)
1393- write(logUnit(7),"(A,I0)") "omega position in state: ", posInState
1394-
1395- ! find chunksize and number of chunks required for initialzation
1396- chunkSize = io_buffer_size / nScalars
1397- nChunks = ceiling( real(nElems, rk) / real(chunkSize, rk) )
1398-
1399- do iChunk = 1, nChunks
1400- ! Number of elements read so far in previous chunks.
1401- elemOff = ( (iChunk-1)*chunksize )
1402- nChunkElems = min(chunkSize, nElems - elemOff)
1403- allocate( omegas(nChunkElems) )
1404- allocate( coord( nChunkElems, 3 ) )
1405-
1406- do iElem = 1, nChunkElems
1407- ! Calculate the coordinates
1408- elemPos = elemOff + iElem
1409- coord(iElem,1:3) = tem_BaryOfId( tree, total( elemPos ) )
1410- end do
1411-
1412- ! obtain omega values from initial condition
1413- omegas(:) = tem_spatial_for( me = ic%ini_state(11), &
1414- & coord = coord, &
1415- & n = nChunkElems )
1416-
1417- do iElem = 1, nChunkElems
1418- elemPos = elemOff + iElem
1419- state( ?IDX?( posInState, elemPos, nScalars, nSize ) ) = omegas(iElem)
1420- ! write(logUnit(9),"(A,F10.5)") 'omega: ', omegas(iElem)
1421- end do ! iElem
1422-
1423- deallocate(omegas)
1424- deallocate(coord)
1425-
1426- end do ! iChunk
1427- ! call tem_abort()
1428-
1429- end subroutine mus_init_omega
1430- ! **************************************************************************** !
14311360
14321361 ! ****************************************************************************** !
14331362 !> Initialize the isothermal acEq flow from density and velocity\n
diff -r b6956abee5af -r 35d71a54ed5d source/mus_harvesting/mus_hvs_construction_module.f90
--- a/source/mus_harvesting/mus_hvs_construction_module.f90 Thu Jun 09 14:16:03 2022 +0200
+++ b/source/mus_harvesting/mus_hvs_construction_module.f90 Tue Aug 30 12:25:45 2022 +0200
@@ -75,11 +75,11 @@
7575 !!
7676 !! # Additional Tasks
7777 !!
78- !! - receive [[tem_construction_module:tem_buildhorizontaldependencies]]
78+ !! - receive [[tem_construction_module:tem_build_horizontaldependencies]]
7979 !! "horizontal"
8080 !! (within a level for the element updates)
81- !! - and [[tem_construction:tem_buildverticaldependencies]] "vertical"
82- !! dependencies (between levels for ghost-interpolations).
81+ !! - and [[tem_construction_module:tem_build_verticaldependencies]] "vertical"
82+ !! dependencies (between levels for ghost-interpolations).
8383 !! - The main state vector and the neighbor lists on which the kernel then
8484 !! acts is created
8585 !! - The MPI buffers are created.
diff -r b6956abee5af -r 35d71a54ed5d source/mus_mrtRelaxation_module.fpp
--- a/source/mus_mrtRelaxation_module.fpp Thu Jun 09 14:16:03 2022 +0200
+++ b/source/mus_mrtRelaxation_module.fpp Tue Aug 30 12:25:45 2022 +0200
@@ -42,6 +42,7 @@
4242 public :: mus_assign_intp_nonEqScalingFacs_ptr
4343 public :: mus_proc_mrt
4444 public :: mus_proc_intp_nonEqScalingFacs
45+ public :: mus_intp_getNonEqScalingFacs
4546 public :: mrt_d2q9
4647 public :: mrt_d3q19
4748 public :: mrt_d3q27
diff -r b6956abee5af -r 35d71a54ed5d source/turbulence/mus_WALE_module.f90
--- a/source/turbulence/mus_WALE_module.f90 Thu Jun 09 14:16:03 2022 +0200
+++ b/source/turbulence/mus_WALE_module.f90 Tue Aug 30 12:25:45 2022 +0200
@@ -82,8 +82,6 @@
8282 real(kind=rk):: gradU(3,3,vlen)
8383 real(kind=rk):: gradU_sqr(3,3,vlen)
8484 integer :: ndims
85- integer :: leftngh(3), rightngh(3)
86- real(kind=rk) :: leftvel(3,3), rightvel(3,3)
8785 integer :: nChunks, iChunks, nChunkElems, low_bound, elempos
8886 ! --------------------------------------------------------------------------
8987 ! viscosity coeff
@@ -187,14 +185,13 @@
187185 real(kind=rk) :: gradU_sqr(2,2,vlen)
188186 real(kind=rk) :: SR(3), Sd(3), onehalf_trSd, Sd_sqr, SR_sqr, OP1, OP2
189187 integer :: ndims
190- integer :: leftngh(2), rightngh(2)
191- real(kind=rk) :: leftvel(2,2), rightvel(2,2)
192188 integer :: nChunks, iChunks, nChunkElems, low_bound, elempos
193189 real(kind=rk) :: visc_coeff
194190 !> gradient of velocity
195191 ! --------------------------------------------------------------------------
196192 visc_coeff = (turbConfig%coeff%C_w * dxL )**2.0_rk
197193 ndims = 2
194+ nChunks = ceiling(real(nSolve, kind=rk) / real(vlen, kind=rk))
198195
199196 do iChunks = 1, nChunks
200197 ! calculate the end number of iElem loop
@@ -213,7 +210,7 @@
213210 ! square of velocity gradient. gradU . gradU
214211 gradU_sqr(:,:,iElem) = matmul(gradU(:,:,iElem), gradU(:,:,iElem))
215212 end do !iElem
216-
213+
217214 do iElem = 1, nChunkElems
218215 elempos = low_bound + iElem
219216 ! traceless symmetric part of the square of the velocity gradient tensor
diff -r b6956abee5af -r 35d71a54ed5d wscript
--- a/wscript Thu Jun 09 14:16:03 2022 +0200
+++ b/wscript Tue Aug 30 12:25:45 2022 +0200
@@ -35,40 +35,53 @@
3535 def subconf(conf):
3636 conf.setenv('')
3737 if not conf.options.no_harvesting:
38- conf.env.build_hvs = True
38+ conf.env.build_hvs = True
3939 else:
40- conf.env.build_hvs = False
40+ conf.env.build_hvs = False
4141
4242 conf.env['with_ext_tdf'] = False
4343 if conf.options.with_ext_tdf:
44- conf.env['with_ext_tdf'] = True
44+ conf.env['with_ext_tdf'] = True
4545
4646 if not conf.env.LIB_PRECICE and conf.env['with_ext_tdf']:
47- conf.setenv('')
48- conf.setenv('cenv')
49- # load c++ compiler for external thermodynamicFactor
50- conf.load('compiler_cxx')
51- # KM:todo set debug CXX flags only for debug variant
52- #conf.env['CXXFLAGS'] = ['-g']
53- conf.check(lib='stdc++', uselib_store='STDCXX', mandatory=True)
54- # Need to have the libraries in the Fortran environment for linking
55- conf.all_envs[''].LIB_STDCXX = conf.env.LIB_STDCXX
56- conf.all_envs[''].LIBPATH_STDCXX = conf.env.LIBPATH_STDCXX
47+ conf.setenv('')
48+ conf.setenv('cenv')
49+ # load c++ compiler for external thermodynamicFactor
50+ conf.load('compiler_cxx')
51+ # KM:todo set debug CXX flags only for debug variant
52+ #conf.env['CXXFLAGS'] = ['-g']
53+ conf.check(lib='stdc++', uselib_store='STDCXX', mandatory=True)
54+ # Need to have the libraries in the Fortran environment for linking
55+ conf.all_envs[''].LIB_STDCXX = conf.env.LIB_STDCXX
56+ conf.all_envs[''].LIBPATH_STDCXX = conf.env.LIBPATH_STDCXX
5757
5858 if conf.options.libdmapp:
59- conf.setenv('')
60- conf.env['libdmapp'] = True
61- conf.check_fc(lib='dmapp', uselib_store='DMAPP', mandatory=False)
59+ conf.setenv('')
60+ conf.env['libdmapp'] = True
61+ conf.check_fc(lib='dmapp', uselib_store='DMAPP', mandatory=False)
6262
6363 def configure(conf):
64+ from waflib import Logs
6465 conf.setenv('')
6566
6667 # Compile flags are set in treelm!
6768 # Check in treelm/wscript for the flags that are set, for the various
6869 # available targets.
70+
6971 conf.recurse('treelm')
7072 subconf(conf)
7173
74+ # Avoid some warnings in gfortran:
75+ if not conf.options.nowarn:
76+ for key, fenv in conf.all_envs.items():
77+ if fenv.FC_NAME == 'GFORTRAN':
78+ fenv.FCFLAGS.append('-Wno-unused-dummy-argument')
79+ fenv.FCFLAGS.append('-fmax-stack-var-size=131072')
80+
81+ Logs.warn('Musubi modified flags:')
82+ Logs.warn('Default flags: {0}'.format(' '.join(conf.all_envs[''].FCFLAGS)))
83+ Logs.warn('Debug flags: {0}'.format(' '.join(conf.all_envs['debug'].FCFLAGS)))
84+
7285 conf.setenv('ford', conf.env)
7386 conf.env.ford_mainpage = 'mus_mainpage.md'
7487
@@ -78,8 +91,8 @@
7891 bld.add_group()
7992
8093 if bld.cmd != 'gendoxy':
81- # Build Treelm and aotus recursive
82- bld.recurse('treelm')
94+ # Build Treelm and aotus recursive
95+ bld.recurse('treelm')
8396
8497 # Build the Musubi objects
8598 build_mus_objs(bld)
@@ -115,24 +128,24 @@
115128 # source files to compute thermodynamic factors
116129 extLib_objs=[]
117130 if bld.env.with_ext_tdf:
118- ext_tf_source = bld.path.ant_glob('external/thermFactor/Density.cpp')
119- ext_tf_source += bld.path.ant_glob('external/thermFactor/eNRTL.cpp')
120- ext_tf_source += bld.path.ant_glob('external/thermFactor/initialize.cpp')
121- ext_tf_source += bld.path.ant_glob('external/thermFactor/Matrixoperation.cpp')
122- mus_sources += bld.path.ant_glob('source/mus_eNRTL_module.f90')
131+ ext_tf_source = bld.path.ant_glob('external/thermFactor/Density.cpp')
132+ ext_tf_source += bld.path.ant_glob('external/thermFactor/eNRTL.cpp')
133+ ext_tf_source += bld.path.ant_glob('external/thermFactor/initialize.cpp')
134+ ext_tf_source += bld.path.ant_glob('external/thermFactor/Matrixoperation.cpp')
135+ mus_sources += bld.path.ant_glob('source/mus_eNRTL_module.f90')
123136
124- bld(
125- features = 'cxx',
126- source = ext_tf_source,
127- target = 'ext_tdf')
128- extLib_objs.append('ext_tdf')
129- extLib_objs.append('STDCXX')
137+ bld(
138+ features = 'cxx',
139+ source = ext_tf_source,
140+ target = 'ext_tdf')
141+ extLib_objs.append('ext_tdf')
142+ extLib_objs.append('STDCXX')
130143 else:
131- mus_sources += bld.path.ant_glob('source/mus_eNRTL_dummy.f90')
132-
144+ mus_sources += bld.path.ant_glob('source/mus_eNRTL_dummy.f90')
145+
133146 if bld.env.libdmapp:
134- extLib_objs.append('DMAPP')
135-
147+ extLib_objs.append('DMAPP')
148+
136149 bld(
137150 features = 'coco fc',
138151 source = mus_sources,
Show on old repository browser