firemodels
Non matching velocities in refinement interface
In debugging the global solver for cases with grid refinement I came across this inconsistency, which seems independent on how much I tighten the velocity tolerance for FFT. Add the two lines to write out average velocities from fine and coarse sides in main.f90, line 1029 (end of corrector):
IF(MY_RANK==0) WRITE(0,) ICYC,'MESH 1=',0.5_EB(MESHES(1)%U(20,1,3)+MESHES(1)%U(20,1,4)) IF(MY_RANK==1) WRITE(0,*) ICYC,'MESH 2=',MESHES(2)%U(0,1,2)
compile and run the case, 2 procs:
&HEAD CHID='test', TITLE='Test refinement region'/
&MESH IJK=20,1,10, XB=0.000,0.200,-0.005,0.005,0.0,0.1 /
&MESH IJK=10,1,5, XB=0.200,0.400,-0.005,0.005,0.0,0.1 /
&TIME T_END= 5. /
&MISC SIMULATION_MODE='DNS', STRATIFICATION=.FALSE., NOISE=.FALSE. /
&PRES SOLVER='FFT', MAX_PRESSURE_ITERATIONS=1000, VELOCITY_TOLERANCE=0.00000001/
&SPEC ID='LJ AIR', VISCOSITY=1.2E-5, BACKGROUND=T /
&SURF ID='WALL', FREE_SLIP=.TRUE., DEFAULT=.TRUE. /
&SURF ID='BLOW', VEL=-5., TAU_V=0. /
&OBST XB=0.0,0.4,-0.005,0.005,0.0,0.0255 /
&VENT PBX=0.0, SURF_ID='BLOW' /
&VENT PBX=0.4, SURF_ID='OPEN' /
&SLCF PBY=0.,QUANTITY='VELOCITY',VECTOR=.TRUE., CELL_CENTERED=T /
&SLCF PBY=0.,QUANTITY='TEMPERATURE', CELL_CENTERED=T /
&SLCF PBY=0.,QUANTITY='CARTESIAN VELOCITY DIVERGENCE', CELL_CENTERED=T /
&TAIL /
The write out added in main.f90 focuses in the interface u-velocity in the coarse face where one of the faces in the fine side is blocked.
The coarse side velocity drops towards zero when tightening the VELOCITY_TOLERANCE, yet the fine side average for that coarse face ranges between 0.15 and 1.5. Not clear what is going on.
I'll dig more on this. We need to understand why we are not matching velocities no matter how much we tighten velocity tolerance for this case. Note this only happens for the interface where the coarse cell has a covered fine face on the other side. If you have an idea what could be going on let me know.
A I recall, we don't check the tolerance at these non-matching boundaries.
I see the fine gas cell next to the coarse cell has a external east face of type INTERPOLATED_BOUNDARY. The coarse cell next to it has a west face of type SOLID_BOUNDARY. It doesn't look like they are connected, unless the normal velocity information from the fine side is being placed on the coarse WALL_CELL u_normal. Regardless, I need to figure out a workaround in my logic here.
Kevin have a look at the case when you have a moment. I expected to see an INTERPOLATED_BOUNDARY WALL cell on the coarse side for the external wall cell of mesh 2, IIG,JJG,KKG=1,1,2 (west face), as it is connected to some fine cells in the gas phase (a single cell in mesh 1 I,J,K = 20,1,4).
During set-up, mesh 2 sees a solid cell in mesh 1 and therefore does not assign the BOUNDARY_TYPE to INTERPOLATED_BOUNDARY. It is possible that if the first fine cell was gas and the second was solid, it might then assume that it is an interpolated boundary. There is no logic to look at all of the cells in the other mesh to decide on the boundary type. This has always been a bit of a gray area. We probably shouldn't allow this, but that ship has sailed.
Ideally other than checking if the other mesh cells are gas phase, we want to see what the BC type of the wall cells on the other side is to make them consistent. This way we can handle thin obstructions in the interface.
For the global matrix building having these be inconsistent poses a problem. The graph of the matrix that states which unknowns are connected to which others is non symmetric, and the matrix comes out wrong.
I'm going to push the changes I have. Other that this issue the refinement treatment I have so far seems fine. I'll throw an error if I pickup a gas cell on the other side of a WALL_CELL type SOLID where EWC%NOM>0.
What FDS does is that it averages the normal component of velocity over the adjacent fine grid cells. But I do not know how you would design the matrix for the pressure solve this way.
Well, I have the divergence operator in the Laplacian is casted in as a sum of fluxes dH/Dxn * Af for each cell faces f. The values of H used here for each flux definition are H in the coarse cell and H in the corresponding fine cell, plus center to center distance (0.5(dxc+dxf) in each cartesian direction. So in all gas-gas cell interfaces the matrix coefficients are tied to the fine side faces and Areas. These are computed straightforwardly at setup time. Once the system is solved the value of H is exchanged in to OMESHES and from then interpolated (extrapolated) to ghost cells both in fine and coarse side such that the average gradient of the solution H is preserved. The gradient of H in the coarse side is the average of the gradients of H in the fine faces. The normal velocities can be projected and will match (coarse velocity average of fine values) by construction. So essentially this does what would be done using the local solvers up to zero velocity error on fine coarse interfaces. In the local solvers we interpolate the value from previous pressure iteration H^k-1 to the faces using coarse/fine cell normal distances, and once we solve we extrapolate the solution H^k to ghost cell. This is what I do once the global system is solved (interpolate to face and extrapolate to ghost cell).
In this special case the gradient of H will be zero in solid fine faces, and dH/Dxn=(H_coarse-H_finei)/(0.5(dxc+dxf)) in interpolated fine faces. It "should" check out same as before.
During set-up, mesh 2 sees a solid cell in mesh 1 and therefore does not assign the
BOUNDARY_TYPEtoINTERPOLATED_BOUNDARY. It is possible that if the first fine cell was gas and the second was solid, it might then assume that it is an interpolated boundary. There is no logic to look at all of the cells in the other mesh to decide on the boundary type. This has always been a bit of a gray area. We probably shouldn't allow this, but that ship has sailed.
This is true. In this case:
&HEAD CHID='test4', TITLE='Test refinement region'/
&MESH IJK=20,1,20, XB=0.000,0.200,-0.005,0.005,-0.1,0.1 /
&MESH IJK=10,1,10, XB=0.200,0.400,-0.005,0.005,-0.1,0.1 /
&TIME T_END= 5. /
&MISC SIMULATION_MODE='DNS', STRATIFICATION=.FALSE., NOISE=.FALSE. /
&PRES SOLVER='FFT', MAX_PRESSURE_ITERATIONS=1000, VELOCITY_TOLERANCE=0.00000001/
&SPEC ID='LJ AIR', VISCOSITY=1.2E-5, BACKGROUND=T /
&SURF ID='WALL', FREE_SLIP=.TRUE., DEFAULT=.TRUE. /
&SURF ID='BLOW', VEL=-5., TAU_V=0. /
&OBST XB=0.0,0.4,-0.005,0.005,-0.0255,0.0255 /
&VENT PBX=0.0, SURF_ID='BLOW' /
&VENT PBX=0.4, SURF_ID='OPEN' /
&SLCF PBY=0.,QUANTITY='VELOCITY',VECTOR=.TRUE., CELL_CENTERED=T /
&SLCF PBY=0.,QUANTITY='H', CELL_CENTERED=T /
&SLCF PBY=0.,QUANTITY='PRESSURE', CELL_CENTERED=T /
&SLCF PBY=0.,QUANTITY='TEMPERATURE', CELL_CENTERED=T /
&SLCF PBY=0.,QUANTITY='CARTESIAN VELOCITY DIVERGENCE', CELL_CENTERED=T /
&TAIL
The lower coarse cell in mesh 2 (1,1,4) has a west external wall cell of type INTERPOLATED, whereas the upper coarse cell next to the refinement in mesh 2 (1,1,7) is of type SOLID. I added a framework to exchange other mesh external wall cell types within the global solver. I'm looking into swapping back the coarse SOLID wall cells to INTERPOLATED, s.t. we have external wall cells from both sides with the same wall type. I tried using the SOLID coarse wall cell and adjust the projection when we have an interpolated wall cell in the OMESH, but this requires tests and changes in many places, making the code overly complicated.
Marcos, I do no think we want to bend over backwards to handle this case in the global solver. We should identify it and throw an Error. Personally, I don't think this should be allowed at all, but if it somehow works now than surely people will complain if we remove it. However, for the global solver we can punt on this.
If this is a snapping issue, then the snapping algorithm could detect this and remove the solid from the interpolated cell.
If this is a snapping issue, then the snapping algorithm could detect this and remove the solid from the interpolated cell.
That should also work. At least we should have the wall cell types match between gas-cells on each side, mass fluxes on that face don't match no matter how many pressure iterations are applied.
For the global solver we exchange the wall cell type from the other mesh in OMESH%EWC_TYPE and use it to set both coarse and fine side mesh WALL INTERPOLATED boundary types to SOLID in cases where the fine wall cells are a mix of INTERPOLATED and SOLID or NULL. This way no matrix graph connections are made between fine and coarse cells on these special cases, and the solver runs with no instability. I do this within the GLOBMAT_SOLVER module. For the local solver the wall boundary types don't match, but there are no instabilities. The velocity error is independent of pressure iterations. I'm fine with closing unless you want to take a stab at making the wall types consistent between coarse and fine side in this special case.