Conversation
* Add D3Q6 in grad data selection * Additional parameter for order of equilibrium * Add a second-order Flekkoy BC
haraldkl
left a comment
There was a problem hiding this comment.
Please also add an example to check this stuff in the recheck, as you mentioned earlier. This pull request is an open branch and in draft mode, you can add stuff and revise the changes as much as you want.
| ! write(*, *) "All populations at Elem 1." | ||
| ! iElem = 1 | ||
| ! elemPos = globBC%elemLvl(iLevel)%elem%val( iElem ) | ||
| ! do iDir = 1, QQ | ||
| ! write(*, *) "Dir ", iDir, ": ", fTmp( (iElem-1)*nScalars + iDir ), "state value: ", & | ||
| ! & state(?FETCH?( iDir, iField, elemPos, QQ, nScalars,nSize, neigh )) | ||
| ! if (globBC%elemLvl(iLevel)%bitmask%val(iDir, iElem)) then | ||
| ! write(*, *) "need to be updated." | ||
| ! invDir = layout%fStencil%cxDirInv(iDir) | ||
| ! write(*, *) "pop at opposite dir: ", - fTmp( (iElem-1)*nScalars + invDir ) | ||
| ! write(*, *) "weight at oppo dir: ", 2._rk*layout%weight( invDir ) | ||
| ! write(*, *) "rho at elem: ", rhoDef(iElem) | ||
| ! else | ||
| ! write(*, *) "don't need to be updated." | ||
| ! end if | ||
| ! end do |
There was a problem hiding this comment.
Please remove commented code before merging into main.
| ! ----------- DEBUG output ------------------------------------------- | ||
| ! write(dbgUnit(6), *) ' fEqPlus: ', fEqPlus | ||
| ! write(dbgUnit(6), *) ' fEqPlusFluid: ', fEqPlusFluid | ||
| ! write(dbgUnit(6), *) ' fPlusFluid: ', fPlusFluid | ||
| ! write(dbgUnit(6), *) ' updated pdf:' | ||
| ! write( dbgUnit(6), *) 'iDir', iDir, 'invDir', invDir, state( & | ||
| ! & ?FETCH?( iDir, iField, globBC%elemLvl(iLevel)%elem%val( iElem ), QQ, nScalars, nSize,neigh )) | ||
| ! -------------------------------------------------------------------------- |
There was a problem hiding this comment.
Please remove commented code before merging into main.
| ! iElem = 1 | ||
| ! elemPos = globBC%elemLvl(iLevel)%elem%val( iElem ) | ||
| ! do iDir = 1, QQ | ||
| ! write(*, *) "Dir ", iDir, ": ", fTmp( (iElem-1)*nScalars + iDir ), "state value: ", & | ||
| ! & state(?FETCH?( iDir, iField, elemPos, QQ, nScalars,nSize, neigh )) | ||
| ! if (globBC%elemLvl(iLevel)%bitmask%val(iDir, iElem)) then | ||
| ! write(*, *) "need to be updated." | ||
| ! else | ||
| ! write(*, *) "don't need to be updated." | ||
| ! end if | ||
| ! end do |
There was a problem hiding this comment.
Please remove commented code before merging into main.
| call aot_get_val( L = conf, & | ||
| & thandle = spc_handle, & | ||
| & key = 'lambda', & | ||
| & val = me%lambda, & | ||
| & ErrCode = iError, & | ||
| & default = 0.25_rk ) |
There was a problem hiding this comment.
Please align the & at the end of the lines.
There was a problem hiding this comment.
BUG: Molecular weight and trt relaxation parameter are both stored in me%lambda
| call aot_get_val( L = conf, & | ||
| & thandle = thandle, & | ||
| & key = 'order', & | ||
| & val = me%order, & | ||
| & default = 'second', & | ||
| & ErrCode = iError ) |
There was a problem hiding this comment.
Please align the & at the end of the lines.
haraldkl
left a comment
There was a problem hiding this comment.
Currently the passive scalar implementation is not working like in the mercurial repository, due to an overhaul of the derived quantity mechanism. The introduced data-type to hold respective function pointers only covers fluid and fluid_incompressible, for the passive scalar we run into an abort at that point. We need to decide on what to properly do here, whether we completely ignore the quantities data structure for that kernel, and just add an empty branch for those types to avoid the abort, or if we remove the aborts.
| getQuantities%rho0Inv_ptr => get_rho0Inv_incompressible | ||
| else | ||
| else | ||
| !! TODO: add pointers to obtain derived quantities for other fluid types |
There was a problem hiding this comment.
Suggestion:
- Replace the abort with a warning that "derived function pointer is not implemented for the scheme kind"
- Rename label_fluid to scheme kind
- Replace if else with select case on scheme kind.
| call tem_abort("get_pdfEq not set for fluid type") | ||
| end if | ||
| case default | ||
| !! TODO: Potentially add pointers to obtain pdfs for other stencils (like d3q6) |
There was a problem hiding this comment.
Thank you for this nice work. It is all good to be merged overall from my point of view, I think. There are just some minor leftovers that should be taken care of. And of course if you do not consider it complete yet, please feel free to work on this branch for as long as you want. Once it is done and ready to be merged from your point of view, remove the draft status from the pull request.
Please be aware of the changes that we are planning for the kernel names in #9 and #10 and the parameters there. This would affect the introduction of the order setting.
| ! write(*, *) "Examine FETCH: " | ||
| ! iElem = 1 | ||
| ! do iDir = 1, varSys%nScalars | ||
| ! pos = ?FETCH?( iDir, 1, iElem, QQ, varSys%nScalars, nElems, neigh) | ||
| ! write(*, *) "iDir ", iDir, ": ", pos | ||
| ! end do | ||
|
|
||
| ! write(*, *) "Examine SAVE: " | ||
| ! iElem = 1 | ||
| ! do iDir = 1, varSys%nScalars | ||
| ! pos = ?SAVE?( iDir, 1, iElem, QQ, varSys%nScalars, nElems, neigh) | ||
| ! write(*, *) iDir, ": ", pos | ||
| ! end do | ||
|
|
||
| ! write(*, *) "Print population before collision" | ||
| ! do iDir = 1, varSys%nScalars | ||
| ! pos = ?FETCH?( iDir, 1, iElem, QQ, varSys%nScalars, nElems, neigh) | ||
| ! write(*, *) "iDir ", iDir, ": ", inState(pos) | ||
| ! end do | ||
|
|
There was a problem hiding this comment.
Please remove commented debug statements before merging.
| ! iElem = 1 | ||
| ! write(*, *) "Print population after collision" | ||
| ! do iDir = 1, varSys%nScalars | ||
| ! pos = ?SAVE?( iDir, 1, iElem, QQ, varSys%nScalars, nElems, neigh) | ||
| ! write(*, *) "iDir ", iDir, ": ", outState(pos) | ||
| ! end do |
There was a problem hiding this comment.
Please remove commented debug statements before merging.
| real(kind=rk) :: ombulkLvl(globalMaxLevels) | ||
| !> relaxation paramete for Nernst-Planck equation | ||
| real(kind=rk) :: omega | ||
| !> relaxation paramete for trt scheme |
There was a problem hiding this comment.
typo: paramete -> parameter
There was a problem hiding this comment.
Note @MikeW097: you can directly read environment variables in Lua.
Like this:
tau = os.getenv('tau')
u_field = os.getenv('u')
order = os.getenv('order')
| bounding_cube = { origin = bc_origin, | ||
| length = length_bnd } | ||
|
|
||
| ebug = {debugMode = true, debugFiles = false, debugMesh='debug/' } |
There was a problem hiding this comment.
I know this is done in several places to "comment out" settings. But I think it is clearer to prepend a NO in front of the renamed variable, like NOdebug.
|
This passes the recheck: We should also include at least one of the example testcases for the regression check in the |
|
Yes but lets note that we dont check passive scalar itself in the regression check. We should add it in future. |
Yeah, might got lost. I wrote that below the recheck results:
So, this is a todo before we merge it in my opinion. |
KannanMasilamani
left a comment
There was a problem hiding this comment.
Nice work, especially the documentation of test case.
I have made few suggestions.
Regarding the order in identify table, I will introduce the variant option as descriped in #10 in next merge. So this branch should be merged after that.
| kind = 'passive_scalar', | ||
| relaxation='trt', | ||
| layout='d2q9', | ||
| order = 'second' |
There was a problem hiding this comment.
``We have observed that there is need to extend this identify table to support several relaxations options for a scheme kind. So, we have fomulated a set of rules see #10.
Here is the suggestion for this case
identify = {
label = 'species',
kind = 'passive_scalar',
relaxation = {
name = 'trt',
variant = 'second_order',
},
layout = 'd2q9'
}
Also see #10 for how to name the compute routine.
| case('pressure_antiBounceBack_pasScal') | ||
| bc( iBnd )%fnct => pressure_antiBounceBack_pasScal |
There was a problem hiding this comment.
No need to introduce 'pressure_antiBounceBack_pasScal' as boundary kind. Use existing 'pressure_antiBounceBack' and select passive scalar routine using scheme kind.
Suggestion:
select case(trim(schemeHeader%kind))
case('fluid')
bc( iBnd )%fnct => pressure_antiBounceBack
case('passive_scalar')
bc( iBnd )%fnct => pressure_antiBounceBack_pasScal
case default
call tem_abort('Unknown scheme kind for pressure_antiBounceBack')
end select
| case('pressure_antiBounceBack_pasScal') | ||
| write(logUnit(1),*) ' Pressure anti bounce back for passive scalar ' & | ||
| & // trim(me( myBCID )%label) | ||
| me( myBCID )%requireNeighBufPost = .true. | ||
| me( myBCID )%nNeighs = 1 | ||
| call tem_load_bc_state( bc = me(myBCID)%BC_states & | ||
| & %pressure, & | ||
| & state_name = 'pressure', & | ||
| & nComp = 1, & | ||
| & conf = conf, & | ||
| & bc_handle = sub_handle, & | ||
| & varDict = me(myBCID)%varDict, & | ||
| & varSys = varSys ) | ||
|
|
There was a problem hiding this comment.
This code duplication is not required if above suggestion is used. So boundary kind is pressure_antiBounceBack and depending on scheme kind current routine is assigned in init_boundary_single.
|
|
||
| public :: mus_append_derVar_lbmPS | ||
| public :: deriveEquilPS_FromMacro | ||
| public :: deriveEquilPS2_FromMacro |
There was a problem hiding this comment.
Suggestion to rename routine deriveEquilPS2ndOrder_FromMacro to be self explaining.
| & applySrc_absorbLayerIncomp | ||
| use mus_derQuanPS_module, only: mus_append_derVar_lbmPS, & | ||
| & deriveEquilPS_FromMacro, & | ||
| & deriveEquilPS2_FromMacro,& |
There was a problem hiding this comment.
See previous comment for renaming this routine.
| case default | ||
| write(logUnit(1),*) 'stencil label = "', trim(label_stencil), '"' | ||
| write(logUnit(1),*) 'fluid type = "', trim(scheme_kind), '"' | ||
| write(logUnit(1),*) "Warning: get_pdfEq not set for fluid type" | ||
| end select |
There was a problem hiding this comment.
Use formating for output. See previous comment.
| write(logUnit(1),*) 'stencil label = "', trim(label_stencil), '"' | ||
| write(logUnit(1),*) 'fluid type = "', trim(scheme_kind), '"' | ||
| write(logUnit(1),*) "Warning: get_pdfEq not set for fluid type" |
There was a problem hiding this comment.
Use formating for output. See previous comment.
| write(logUnit(1),*) 'stencil label = "', trim(label_stencil), '"' | ||
| write(logUnit(1),*) 'fluid type = "', trim(scheme_kind), '"' | ||
| write(logUnit(1),*) "Warning: get_pdfEq not set for fluid type" |
There was a problem hiding this comment.
Use formating for output. See previous comment.
| !> equilibrium order Ex: first, second | ||
| character(len=labelLen) :: order |
There was a problem hiding this comment.
Should be replaced with relaxation variant. This option is to be introduced in change merge.
| ! get order of equilibrium | ||
| call aot_get_val( L = conf, & | ||
| & thandle = thandle, & | ||
| & key = 'order', & | ||
| & val = me%order, & | ||
| & default = 'second', & | ||
| & ErrCode = iError ) | ||
|
|
There was a problem hiding this comment.
It is not required after relaxation variant options is introduced.
|
Closing this PR, as it has been superseded by #25. |
This pull request provides for some updates and improvements to the passive scalar implementation.