diff --git a/example/data/1024x256_initial_hz.png b/example/data/1024x256_initial_hz.png new file mode 100644 index 000000000..646dbf009 Binary files /dev/null and b/example/data/1024x256_initial_hz.png differ diff --git a/example/data/128x32.png b/example/data/128x32.png new file mode 100644 index 000000000..6814d3e02 Binary files /dev/null and b/example/data/128x32.png differ diff --git a/example/data/32x32.png b/example/data/32x32.png new file mode 100644 index 000000000..f4db923df Binary files /dev/null and b/example/data/32x32.png differ diff --git a/example/flow/2d/cave.xml b/example/flow/2d/cave.xml new file mode 100644 index 000000000..5f9889d7d --- /dev/null +++ b/example/flow/2d/cave.xml @@ -0,0 +1,45 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + Solver$Fields$Height[] = rnorm(64*32, mean = 0.02048, sd=0.000002048) + + Solver$Actions$InitFromFields() + Solver$Actions$InitFromFields() + + + + + + + + + + + + + + + + diff --git a/example/flow/2d/fr.xml b/example/flow/2d/fr.xml new file mode 100644 index 000000000..6f55b0d85 --- /dev/null +++ b/example/flow/2d/fr.xml @@ -0,0 +1,76 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + +<<<<<<< HEAD + + + + + + + + + + + + + + + + + + + + +======= + + + + + + + + + + + + + + + + + + + + + + +>>>>>>> soluteMRT + + + diff --git a/example/flow/2d/fracture.xml b/example/flow/2d/fracture.xml new file mode 100644 index 000000000..23d14402f --- /dev/null +++ b/example/flow/2d/fracture.xml @@ -0,0 +1,73 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/example/flow/2d/fracture_dissolution.xml b/example/flow/2d/fracture_dissolution.xml new file mode 100644 index 000000000..e8b04dea4 --- /dev/null +++ b/example/flow/2d/fracture_dissolution.xml @@ -0,0 +1,92 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + Solver$Fields$Height[] = rnorm(100, mean = 0.02, sd=0.01) + + Solver$Actions$InitFromFields() + Solver$Actions$InitFromFields() + Solver$Actions$InitFromFields() + + + + + + + + + + + + + + + + + + + + + + diff --git a/models/flow/d2q9_fracture/Description.md b/models/flow/d2q9_fracture/Description.md new file mode 100644 index 000000000..e9bf8e383 --- /dev/null +++ b/models/flow/d2q9_fracture/Description.md @@ -0,0 +1,8 @@ +The `d2q9` model is a basic 2D Lattice Boltzmann Method model for flow +simulation. + +It is a implementation of a simple Multiple-Relaxation-Time collision +operator. + +The model has one [option](/basics/options/): BC. This option turns on +custom fields for setting non-standard profiles on inlets and other things. diff --git a/models/flow/d2q9_fracture/Dynamics.R b/models/flow/d2q9_fracture/Dynamics.R new file mode 100644 index 000000000..8b7ba4fc2 --- /dev/null +++ b/models/flow/d2q9_fracture/Dynamics.R @@ -0,0 +1,125 @@ +# Density - table of variables of LB Node to stream +# name - variable name to stream +# dx,dy,dz - direction of streaming +# comment - additional comment + +AddDensity(name = "f[0]", dx = 0, dy = 0, group = "f") +AddDensity(name = "f[1]", dx = 1, dy = 0, group = "f") +AddDensity(name = "f[2]", dx = 0, dy = 1, group = "f") +AddDensity(name = "f[3]", dx = -1, dy = 0, group = "f") +AddDensity(name = "f[4]", dx = 0, dy = -1, group = "f") +AddDensity(name = "f[5]", dx = 1, dy = 1, group = "f") +AddDensity(name = "f[6]", dx = -1, dy = 1, group = "f") +AddDensity(name = "f[7]", dx = -1, dy = -1, group = "f") +AddDensity(name = "f[8]", dx = 1, dy = -1, group = "f") + +AddDensity(name = "h", dx = 0, dy = 0, group = "HZ") + +AddDensity(name = "c[0]", dx = 0, dy = 0, group = "c") +AddDensity(name = "c[1]", dx = 1, dy = 0, group = "c") +AddDensity(name = "c[2]", dx = 0, dy = 1, group = "c") +AddDensity(name = "c[3]", dx = -1, dy = 0, group = "c") +AddDensity(name = "c[4]", dx = 0, dy = -1, group = "c") +AddDensity(name = "c[5]", dx = 1, dy = 1, group = "c") +AddDensity(name = "c[6]", dx = -1, dy = 1, group = "c") +AddDensity(name = "c[7]", dx = -1, dy = -1, group = "c") +AddDensity(name = "c[8]", dx = 1, dy = -1, group = "c") + +# AddField(name="f[1]", dx=1); + +# THIS QUANTITIES ARE NEEDED FOR PYTHON INTEGRATION EXAMPLE +# COMMENT OUT FOR PERFORMANCE +# If present thei are used: +# As VelocityX/Y for Boundary conditions +# As mass force (+ GravitationX/Y) in fluid +if (Options$bc) { + AddDensity(name = "BC[0]", group = "BC", parameter = TRUE) + AddDensity(name = "BC[1]", group = "BC", parameter = TRUE) +} + +# Quantities - table of fields that can be exported from the LB lattice (like density, velocity etc) +# name - name of the field +# type - C type of the field, "real_t" - for single/double float, and "vector_t" for 3D vector single/double float +# Every field must correspond to a function in "Dynamics.c". +# If one have filed [something] with type [type], one have to define a function: +# [type] get[something]() { return ...; } + +AddQuantity(name = "Rho", unit = "kg/m3") +AddQuantity(name = "P", unit = "Pa") +AddQuantity(name = "U", unit = "m/s", vector = T) +AddQuantity(name = "H", unit = "m") +AddQuantity(name = "C", unit = "kg/m3") + +# Stages +AddDensity(name = "Height_parameter", group = "init", parameter = TRUE) +AddStage(name = "BaseIteration", main = "Run", load.densities = TRUE, save.fields = TRUE) +AddStage(name = "BaseInitFromFields", main = "InitFromFields", load.densities = TRUE, save.fields = TRUE) +AddStage(name = "BasehEvolution", main = "hEvolution", load.densities = TRUE, save.fields = TRUE) +AddAction(name = "Iteration", "BaseIteration") +AddAction(name = "InitFromFields", "BaseInitFromFields") +AddAction(name = "hEvolution", "BasehEvolution") + +# Settings - table of settings (constants) that are taken from a .xml file +# name - name of the constant variable +# comment - additional comment +# You can state that another setting is 'derived' from this one stating for example: RelaxationRate='1.0/(3*Viscosity + 0.5)' + +AddSetting(name = "RelaxationRate", S2 = "1-RelaxationRate", comment = "one over relaxation time") +AddSetting(name = "Viscosity", RelaxationRate = "1.0/(3*Viscosity + 0.5)", default = 0.16666666, comment = "viscosity") +AddSetting(name = "VelocityX", default = 0, comment = "inlet/outlet/init velocity", zonal = T) +AddSetting(name = "VelocityY", default = 0, comment = "inlet/outlet/init velocity", zonal = T) +AddSetting(name = "Pressure", default = 0, comment = "inlet/outlet/init density", zonal = T) +AddSetting(name = "GravitationX", default = 0, comment = "body/external acceleration", zonal = T) +AddSetting(name = "GravitationY", default = 0, comment = "body/external acceleration", zonal = T) + +AddSetting(name = "S2", default = "0", comment = "MRT Sx") +AddSetting(name = "S3", default = "0", comment = "MRT Sx") +AddSetting(name = "S4", default = "0", comment = "MRT Sx") +AddSetting(name = "MagicNumber", default = 0.25, comment = "TRT magic number") + +AddSetting(name = "Height", default = 1, zonal = TRUE) + +AddSetting(name = "ConcentrationRelaxationRate", SC2 = "1-ConcentrationRelaxationRate", default = 0, comment = "one over concentration relaxation time") +AddSetting(name = "Diffusivity", ConcentrationRelaxationRate = "1.0/(3*Diffusivity+0.5)", default = 0.166666667, comment = "diffusivity") +AddSetting(name = "Concentration", default = 0, comment = "inlet/outlet/init concentration", zonal = T) +AddSetting(name = "Saturation", default = 1, comment = "Saturation concentration") +AddSetting(name = "ReactionConstant", default = 1, comment = "Reaction speed coefficient") +AddSetting(name = "SolidConcentration", default = 1, comment = "Solid concentration") +AddSetting(name = "Sherwood", default = 8, comment = "Sherwood number") +AddSetting(name = "hTimeStep", default = 1, comment = "h evolution time step in pseudo-steady state approach") + + +AddSetting(name = "SC2", default = "0", comment = "MRT concentration SCx") +AddSetting(name = "SC3", default = "0", comment = "MRT concentration SCx") +AddSetting(name = "SC4", default = "0", comment = "MRT concentration SCx") +AddSetting(name = "MagicNumber_D", default = 0.25, comment = "TRT concetration magic number") + +# Globals - table of global integrals that can be monitored and optimized +AddGlobal(name = "PressureLoss", comment = "pressure loss", unit = "1mPa") +AddGlobal(name = "OutletFlux", comment = "pressure loss", unit = "1m2/s") +AddGlobal(name = "InletFlux", comment = "pressure loss", unit = "1m2/s") + +# Node types for boundaries +AddNodeType(name = "EPressure", group = "BOUNDARY") +AddNodeType(name = "WPressure", group = "BOUNDARY") + +AddNodeType(name = "NVelocity", group = "BOUNDARY") +AddNodeType(name = "SVelocity", group = "BOUNDARY") +AddNodeType(name = "WVelocity", group = "BOUNDARY") +AddNodeType(name = "EVelocity", group = "BOUNDARY") + +AddNodeType(name = "NSymmetry", group = "BOUNDARY") +AddNodeType(name = "SSymmetry", group = "BOUNDARY") + +AddNodeType(name = "WPressureWConcentration", group = "BOUNDARY") +AddNodeType(name = "EPressureEConcentration", group = "BOUNDARY") + +AddNodeType(name = "Inlet", group = "OBJECTIVE") +AddNodeType(name = "Outlet", group = "OBJECTIVE") +AddNodeType(name = "Solid", group = "BOUNDARY") +AddNodeType(name = "Wall", group = "BOUNDARY") +AddNodeType(name = "MRT", group = "COLLISION") +AddNodeType(name = "IncompressibleMRT", group = "COLLISION") + +# +AddNodeType(name = "Internal", group = "OBJECTIVE") \ No newline at end of file diff --git a/models/flow/d2q9_fracture/Dynamics.c.Rt b/models/flow/d2q9_fracture/Dynamics.c.Rt new file mode 100644 index 000000000..36c011ff7 --- /dev/null +++ b/models/flow/d2q9_fracture/Dynamics.c.Rt @@ -0,0 +1,525 @@ + + +// -- Macroscopic quantieties getters -- +CudaDeviceFunction real_t getRho(){ + return ; +} + +CudaDeviceFunction real_t getP(){ + return ( - 1.)/3.; +} + +CudaDeviceFunction vector_t getU(){ + real_t d = ; + vector_t u; + + u.x /= d; + u.y /= d; + + if (!IamBOUNDARY) { + u.x += BC[0]*0.5; + u.y += BC[1]*0.5; + } + + // half of the source added + u.x += GravitationX*0.5; + u.y += GravitationY*0.5; + + if ((NodeType & NODE_OBJECTIVE) == NODE_Internal) { + real_t nu = -(S2 + 1)/(6*S2 - 6); + real_t K = 12. * nu / h / h; + + u.x = u.x - K*u.x*0.5; + u.y = u.y - K*u.y*0.5; + } + + u.z = 0.0; + return u; +} + +CudaDeviceFunction real_t getH(){ + return h; +} + +CudaDeviceFunction real_t getC(){ + real_t C = ; + + // half of the source added + if ((NodeType & NODE_OBJECTIVE) == NODE_Internal) { + if(C < Saturation) { + // real_t ReactionConstant_eff = ReactionConstant * Diffusivity * Sherwood / (Diffusivity * Sherwood + 2*ReactionConstant*h); + real_t Q = ReactionConstant * (Saturation - C); + C += Q; + } + } + + return C; +} + +CudaDeviceFunction float2 Color() { + float2 ret; + vector_t u = getU(); + ret.x = sqrt(u.x*u.x + u.y*u.y); + if (NodeType == NODE_Solid){ + ret.y = 0; + } else { + ret.y = 1; + } + return ret; +} + +// -- Initialisation and running simulation -- +CudaDeviceFunction void SetEquilibrum(real_t rho, real_t Jx, real_t Jy, real_t C, real_t ux, real_t uy) { + switch (NodeType & NODE_COLLISION) { + case NODE_MRT: + + break; + case NODE_IncompressibleMRT: + + break; + } + + + + if ( IamBOUNDARY ) { + BC[0] = Jx / rho; + BC[1] = Jy / rho; + } else { + BC[0] = 0; + BC[1] = 0; + } + +} + +CudaDeviceFunction void Init() { + real_t rho, ux, uy, C; + rho = (1+Pressure*3); + ux = VelocityX; + uy = VelocityY; + C = Concentration; + if (IamWall){ + rho = 1; + ux = 0; + uy = 0; + } + SetEquilibrum( + rho, + ux*rho, + uy*rho, + C, + ux, + uy + ); + + h = Height; +} + +CudaDeviceFunction void InitFromFields() { + h = Height_parameter; +} + +CudaDeviceFunction void hEvolution() { + real_t C = ; + real_t Q = ReactionConstant * (Saturation - C); + h += hTimeStep*2*Q/SolidConcentration; +} + +CudaDeviceFunction void Run() { + doBC(); + doCollision(); +} + +CudaDeviceFunction void doBC() { + switch (NodeType & NODE_BOUNDARY) { + case NODE_Solid: + case NODE_Wall: + BounceBack(); + break; + case NODE_EVelocity: + EVelocity(); + break; + case NODE_WPressure: + WPressure(); + break; + case NODE_WVelocity: + WVelocity(); + break; + case NODE_EPressure: + EPressure(); + break; + case NODE_NVelocity: + NVelocity(); + break; + case NODE_SVelocity: + SVelocity(); + break; + case NODE_NSymmetry: + NSymmetry(); + break; + case NODE_SSymmetry: + SSymmetry(); + break; + case NODE_WPressureWConcentration: + WPressureWConcentration(); + break; + case NODE_EPressureEConcentration: + EPressureEConcentration(); + break; + } +} + +CudaDeviceFunction void doCollision() { + switch (NodeType & NODE_COLLISION) { + case NODE_MRT: + CollisionMRT(); + break; + case NODE_IncompressibleMRT: + CollisionIncompressibleMRT(); + break; + } + + ConcentrationCollisionSRT(); +} + +// -- BCs -- +CudaDeviceFunction void BounceBack() +{ + + if ( temp != 1+Pressure*3 ) { + + } else { + + } + +} + +CudaDeviceFunction void EVelocity() +{ + +} + +CudaDeviceFunction void WPressure() +{ + +} + +CudaDeviceFunction void WVelocity() +{ + +} + +CudaDeviceFunction void EPressure() +{ + +} + +CudaDeviceFunction void NVelocity() +{ + +} + +CudaDeviceFunction void SVelocity() +{ + +} + +CudaDeviceFunction void NSymmetry() +{ + +} + +CudaDeviceFunction void SSymmetry() +{ + +} + +CudaDeviceFunction void WPressureWConcentration() +{ + +//d2q5 +// c[1] = Concentration - c[0] - c[2] - c[3] - c[4]; +//d2q9 + real_t C = 6.*( Concentration - (c[0] + c[2] + c[4] + c[3] + c[7] + c[6]) ); + c[1] = (1./9.) * C; + c[5] = (1./36.) * C; + c[8] = (1./36.) * C; +} + +CudaDeviceFunction void EPressureEConcentration() +{ + +//d2q5 +// c[3] = Concentration - c[0] - c[1] - c[2] - c[4]; +//d2q9 + real_t C = 6.*( Concentration - ( c[0] + c[2] + c[4] + c[1] + c[5] + c[8] ) ); + c[3] = (1./ 9.) * C; + c[7] = (1./36.) * C; + c[6] = (1./36.) * C; +} + +// -- Colision and streaming -- +CudaDeviceFunction void CollisionMRT() +{ + 1 + R[!selR] = fMRT$Req[!selR] +?> + real_t ; + real_t Usq=0; + real_t omega_odd = 2.*(1.+S2)/((1.-S2)*(4.*MagicNumber-1.)+2.); + real_t S_odd = 1. - omega_odd; + + + switch (NodeType & NODE_OBJECTIVE) { + case NODE_Outlet: + + AddToOutletFlux(Jx/rho/rho); + AddToPressureLoss(-(Jx/rho)/rho*((rho-1.)/3. + Usq/rho/2.)); + break; + case NODE_Inlet: + + AddToInletFlux(Jx/rho/rho); + AddToPressureLoss((Jx/rho)/rho*((rho-1.)/3. + Usq/rho/2.)); + break; + } + + + + + if (!IamBOUNDARY) { + Jx = Jx + (GravitationX + BC[0])*rho ; + Jy = Jy + (GravitationY + BC[1])*rho ; + } + + Jx = Jx + GravitationX*rho; + Jy = Jy + GravitationY*rho; + + + if ((NodeType & NODE_OBJECTIVE) == NODE_Internal) { + real_t nu = -(S2 + 1)/(6*S2 - 6); + real_t K = 12. * nu / h / h; + Jx = Jx - K*Jx; + Jy = Jy - K*Jy; + } + +} + +CudaDeviceFunction void CollisionIncompressibleMRT() +{ + 1 + R[!selR] = fMRT_incomp$Req[!selR] +?> + real_t ; + real_t Usq=0; + + + switch (NodeType & NODE_OBJECTIVE) { + case NODE_Outlet: + + AddToOutletFlux(ux/rho); + AddToPressureLoss(-ux/rho*((rho-1.)/3. + Usq/rho/2.)); + break; + case NODE_Inlet: + + AddToInletFlux(ux/rho); + AddToPressureLoss(ux/rho*((rho-1.)/3. + Usq/rho/2.)); + break; + } + + + + + if (!IamBOUNDARY) { + ux = ux + GravitationX + BC[0]; + uy = uy + GravitationY + BC[1]; + } + + ux = ux + GravitationX*rho; + uy = uy + GravitationY*rho; + + + if ((NodeType & NODE_OBJECTIVE) == NODE_Internal) { + real_t nu = -(S2 + 1)/(6*S2 - 6); + real_t K = 12. * nu / h / h; + ux = ux - K*ux; + uy = uy - K*uy; + } + +} + +CudaDeviceFunction void ConcentrationCollisionMRT() { + 0 + RC[!selCR] = cMRT$Req[!selCR] +?> + real_t ; + real_t ux, uy; + ux = f[8] - f[7] - f[6] + f[5] - f[3] + f[1]; + uy = -f[8] - f[7] + f[6] + f[5] - f[4] + f[2]; + real_t omega_odd = 2.*(2.-ConcentrationRelaxationRate)/(ConcentrationRelaxationRate*(4.*MagicNumber_D-1.)+2.); + real_t SC_odd = 1. - omega_odd; + + if ((NodeType & NODE_OBJECTIVE) == NODE_Internal) { + if(C < Saturation) { + // real_t ReactionConstant_eff = ReactionConstant * Diffusivity * Sherwood / (Diffusivity * Sherwood + 2*ReactionConstant*h); + real_t Q = ReactionConstant * (Saturation - C); + C = C + 2*Q; + } else { + C = Saturation; + } + } + +} + +CudaDeviceFunction void ConcentrationCollisionSRT() { + real_t C = ; + real_t ux, uy; + ux = f[8]-f[7]-f[6]+f[5]-f[3]+f[1]; + uy = -f[8]-f[7]+f[6]+f[5]-f[4]+f[2]; + + real_t c_eq[9]; + + + + if ((NodeType & NODE_OBJECTIVE) == NODE_Internal) { + if(C < Saturation) { + // real_t ReactionConstant_eff = ReactionConstant * Diffusivity * Sherwood / (Diffusivity * Sherwood + 2*ReactionConstant*h); + real_t Q = ReactionConstant * (Saturation - C); + c[0] += 2*Q * 4./9.; + c[1] += 2*Q * 1./9; + c[2] += 2*Q * 1./9.; + c[3] += 2*Q * 1./9.; + c[4] += 2*Q * 1./9.; + c[5] += 2*Q * 1./36.; + c[6] += 2*Q * 1./36.; + c[7] += 2*Q * 1./36.; + c[8] += 2*Q * 1./36.; + // h += 2*Q/SolidConcentration; + } else { + c[0] = Saturation * 4./9.; + c[1] = Saturation * 1./9; + c[2] = Saturation * 1./9.; + c[3] = Saturation * 1./9.; + c[4] = Saturation * 1./9.; + c[5] = Saturation * 1./36.; + c[6] = Saturation * 1./36.; + c[7] = Saturation * 1./36.; + c[8] = Saturation * 1./36.; + } + } +} \ No newline at end of file diff --git a/models/flow/d2q9_fracture/conf.mk b/models/flow/d2q9_fracture/conf.mk new file mode 100644 index 000000000..3e710337e --- /dev/null +++ b/models/flow/d2q9_fracture/conf.mk @@ -0,0 +1,3 @@ +ADJOINT=0 +TEST=TRUE +OPT="bc*autosym" diff --git a/models/flow/d2q9_fracture_dissolution/Dynamics.R b/models/flow/d2q9_fracture_dissolution/Dynamics.R new file mode 100644 index 000000000..449b44dbd --- /dev/null +++ b/models/flow/d2q9_fracture_dissolution/Dynamics.R @@ -0,0 +1,104 @@ +# Densieties and fields +AddDensity( name="f[0]", dx= 0, dy= 0, group="f") +AddDensity( name="f[1]", dx= 1, dy= 0, group="f") +AddDensity( name="f[2]", dx= 0, dy= 1, group="f") +AddDensity( name="f[3]", dx=-1, dy= 0, group="f") +AddDensity( name="f[4]", dx= 0, dy=-1, group="f") +AddDensity( name="f[5]", dx= 1, dy= 1, group="f") +AddDensity( name="f[6]", dx=-1, dy= 1, group="f") +AddDensity( name="f[7]", dx=-1, dy=-1, group="f") +AddDensity( name="f[8]", dx= 1, dy=-1, group="f") + +AddDensity( name="c[0]", dx= 0, dy= 0, group="c") +AddDensity( name="c[1]", dx= 1, dy= 0, group="c") +AddDensity( name="c[2]", dx= 0, dy= 1, group="c") +AddDensity( name="c[3]", dx=-1, dy= 0, group="c") +AddDensity( name="c[4]", dx= 0, dy=-1, group="c") +AddDensity( name="c[5]", dx= 1, dy= 1, group="c") +AddDensity( name="c[6]", dx=-1, dy= 1, group="c") +AddDensity( name="c[7]", dx=-1, dy=-1, group="c") +AddDensity( name="c[8]", dx= 1, dy=-1, group="c") + +AddDensity( name="h", dx=0, dy=0, group="HZ") + +# Stages +AddDensity(name="Height_parameter", group="init", parameter=TRUE) +AddStage(name="BaseIteration", main="Run", load.densities=TRUE, save.fields=TRUE) +AddStage(name="BaseInitFromFields", main="InitFromFields", load.densities=TRUE, save.fields=TRUE) +AddAction(name="Iteration", "BaseIteration") +AddAction(name="InitFromFields", "BaseInitFromFields") + +# Output quantieties +AddQuantity(name="Rho", unit="kg/m3") +AddQuantity(name="P", unit="Pa") +AddQuantity(name="U", unit="m/s",vector=T) +AddQuantity(name="H", unit="m") +AddQuantity(name="C", unit="kg/m3") + +# Settings +AddSetting(name="RelaxationRate", S2='1-RelaxationRate', comment='one over relaxation time') +AddSetting(name="Viscosity", RelaxationRate='1.0/(3*Viscosity + 0.5)', default=0.166666667, comment='viscosity') +AddSetting(name="VelocityX", default=0, comment='inlet/outlet/init velocity', zonal=T) +AddSetting(name="VelocityY", default=0, comment='inlet/outlet/init velocity', zonal=T) +AddSetting(name="Pressure", default=0, comment='inlet/outlet/init density', zonal=T) +AddSetting(name="GravitationX",default=0, comment='body/external acceleration', zonal=TRUE) +AddSetting(name="GravitationY",default=0, comment='body/external acceleration', zonal=TRUE) + +AddSetting(name="S2", default="0", comment='MRT Sx') +AddSetting(name="S3", default="0", comment='MRT Sx') +AddSetting(name="S4", default="0", comment='MRT Sx') + +AddSetting(name="nubuffer", default=0.01, comment='Viscosity in the buffer layer (cumulant)') + +AddSetting(name="Height", default=1, zonal=TRUE) + +AddSetting(name="ConcentrationRelaxationRate", SC2="1-ConcentrationRelaxationRate", default=0, comment='one over concentration relaxation time') +AddSetting(name="Diffusivity", ConcentrationRelaxationRate="1.0/(3*Diffusivity+0.5)", default=0.166666667, comment='diffusivity') +AddSetting(name="Concentration", default=0, comment='inlet/outlet/init concentration', zonal=T) +AddSetting(name="Saturation", default=1) +AddSetting(name="k",default=1, comment='Reaction speed coefficient') +AddSetting(name="SolidConcentration", default=100, comment='Solid concentration') + +AddSetting(name="SC2", default="0", comment='MRT concentration Sx') +AddSetting(name="SC3", default="0", comment='MRT concentration Sx') +AddSetting(name="SC4", default="0", comment='MRT concentration Sx') +AddSetting(name="MagicNumber_D", default=0.25, comment='TRT concetration magic number') + + +# AddSetting(name="ConcentrationRelaxationRate_even", default=0, comment='TRT even concetration relaxation rate') +# AddSetting(name="ConcentrationRelaxationRate_odd", ConcentrationRelaxationRate_even="(1.0 - 0.5*ConcentrationRelaxationRate_odd)/((MagicNumber_D-0.25)*ConcentrationRelaxationRate_odd + 0.5)", default=0, comment='TRT odd concetration relaxation rate') + +# AddSetting(name="SC_odd", default="0", comment='MRT CSx') +# AddSetting(name="Diffusivity", ConcentrationRelaxationRate='1.0/(3*Diffusivity+0.5)', default=0.16666666, comment='Diffusivity') +# AddSetting(name="Diffusivity", ConcentrationRelaxationRate_odd="1.0/(3*Diffusivity+0.5)", default=0.16666666, comment='Diffusivity') + + +# Globals - table of global integrals that can be monitored and optimized +AddGlobal(name="PressureLoss", comment='pressure loss', unit="1mPa") +AddGlobal(name="OutletFlux", comment='pressure loss', unit="1m2/s") +AddGlobal(name="InletFlux", comment='pressure loss', unit="1m2/s") + +# Node types for boundaries +AddNodeType(name="EPressure", group="BOUNDARY") +AddNodeType(name="WPressure", group="BOUNDARY") + +AddNodeType(name="NVelocity", group="BOUNDARY") +AddNodeType(name="SVelocity", group="BOUNDARY") +AddNodeType(name="WVelocity", group="BOUNDARY") +AddNodeType(name="EVelocity", group="BOUNDARY") + +AddNodeType(name="WPressureWConcentration", group="BOUNDARY") +AddNodeType(name="EPressureEConcentration", group="BOUNDARY") +AddNodeType(name="WVelocityWConcentration", group="BOUNDARY") + +AddNodeType(name="NSymmetry", group="BOUNDARY") +AddNodeType(name="SSymmetry", group="BOUNDARY") + +AddNodeType(name="Inlet", group="OBJECTIVE") +AddNodeType(name="Outlet", group="OBJECTIVE") +AddNodeType(name="Solid", group="BOUNDARY") +AddNodeType(name="Wall", group="BOUNDARY") +AddNodeType(name="MRT", group="COLLISION") +AddNodeType(name="IncompressibleSRT", group="COLLISION") +AddNodeType(name="IncompressibleMRT", group="COLLISION") +AddNodeType(name="Reaction", group="SOURCE") diff --git a/models/flow/d2q9_fracture_dissolution/Dynamics.c.Rt b/models/flow/d2q9_fracture_dissolution/Dynamics.c.Rt new file mode 100644 index 000000000..e1749bbd7 --- /dev/null +++ b/models/flow/d2q9_fracture_dissolution/Dynamics.c.Rt @@ -0,0 +1,497 @@ + + +CudaDeviceFunction void Init() { + real_t rho = (1 + Pressure*3); + real_t C = Concentration; + vector_t u; + u.x = VelocityX; + u.y = VelocityY; + + SetEquilibrum_f(rho, rho*u.x, rho*u.y); + SetEquilibrum_c(C, C*u.x, C*u.y); + + h = Height; +} + +CudaDeviceFunction void InitFromFields() { + h = Height_parameter; +} + +CudaDeviceFunction void Run() { + switch (NodeType & NODE_BOUNDARY) { + case NODE_Solid: + case NODE_Wall: + BounceBack(); + break; + case NODE_EVelocity: + EVelocity(); + break; + case NODE_WPressure: + WPressure(); + break; + case NODE_WPressureWConcentration: + WPressureWConcentration(); + break; + case NODE_WVelocity: + WVelocity(); + break; + case NODE_WVelocityWConcentration: + WVelocityWConcentration(); + break; + case NODE_EPressure: + EPressure(); + break; + case NODE_EPressureEConcentration: + EPressureEConcentration(); + break; + case NODE_NVelocity: + NVelocity(); + break; + case NODE_SVelocity: + SVelocity(); + break; + case NODE_NSymmetry: + NSymmetry(); + break; + case NODE_SSymmetry: + SSymmetry(); + break; + } + switch (NodeType & NODE_COLLISION) + { + case NODE_MRT: + CollisionMRT(); + break; + case NODE_IncompressibleSRT: + CollisionIncompressibleSRT(); + break; + case NODE_IncompressibleMRT: + CollisionIncompressibleMRT(); + break; + } +} + +CudaDeviceFunction void SetEquilibrum_f(real_t rho, real_t Jx, real_t Jy) +{ + +} + +CudaDeviceFunction void SetEquilibrum_c(real_t C, real_t JCx, real_t JCy) +{ + +} + +CudaDeviceFunction void SetIncompressibleEquilibrium(real_t phi[9], real_t Phi, vector_t u){ + phi[0] = ( 2*Phi + ( -u.y*u.y - u.x*u.x )*3. )*2./9.; + phi[1] = ( 2*Phi + ( -u.y*u.y + ( 1 + u.x )*u.x*2. )*3. )/18.; + phi[2] = ( 2*Phi + ( -u.x*u.x + ( 1 + u.y )*u.y*2. )*3. )/18.; + phi[3] = ( 2*Phi + ( -u.y*u.y + ( -1 + u.x )*u.x*2. )*3. )/18.; + phi[4] = ( 2*Phi + ( -u.x*u.x + ( -1 + u.y )*u.y*2. )*3. )/18.; + phi[5] = ( Phi + ( ( 1 + u.y )*u.y + ( 1 + u.x + u.y*3. )*u.x )*3. )/36.; + phi[6] = ( Phi + ( ( 1 + u.y )*u.y + ( -1 + u.x - u.y*3. )*u.x )*3. )/36.; + phi[7] = ( Phi + ( ( -1 + u.y )*u.y + ( -1 + u.x + u.y*3. )*u.x )*3. )/36.; + phi[8] = ( Phi + ( ( -1 + u.y )*u.y + ( 1 + u.x - u.y*3. )*u.x )*3. )/36.; +} + +CudaDeviceFunction void SetEquilibrium(real_t phi[9], real_t Phi, vector_t u){ + phi[0] = ( 2. + ( -u.y*u.y - u.x*u.x )*3. )*Phi*2./9.; + phi[1] = ( 2. + ( -u.y*u.y + ( 1 + u.x )*u.x*2. )*3. )*Phi/18.; + phi[2] = ( 2. + ( -u.x*u.x + ( 1 + u.y )*u.y*2. )*3. )*Phi/18.; + phi[3] = ( 2. + ( -u.y*u.y + ( -1 + u.x )*u.x*2. )*3. )*Phi/18.; + phi[4] = ( 2. + ( -u.x*u.x + ( -1 + u.y )*u.y*2. )*3. )*Phi/18.; + phi[5] = ( 1. + ( ( 1 + u.y )*u.y + ( 1 + u.x + u.y*3. )*u.x )*3. )*Phi/36.; + phi[6] = ( 1. + ( ( 1 + u.y )*u.y + ( -1 + u.x - u.y*3. )*u.x )*3. )*Phi/36.; + phi[7] = ( 1. + ( ( -1 + u.y )*u.y + ( -1 + u.x + u.y*3. )*u.x )*3. )*Phi/36.; + phi[8] = ( 1. + ( ( -1 + u.y )*u.y + ( 1 + u.x - u.y*3. )*u.x )*3. )*Phi/36.; +} + +CudaDeviceFunction void BounceBack() +{ + +} + +CudaDeviceFunction void EVelocity() +{ + +} + +CudaDeviceFunction void WPressure() +{ + +} + +CudaDeviceFunction void WPressureWConcentration() +{ + +//d2q5 +//c[1] = Concentration - c[0] - c[2] - c[3] - c[4]; +//d2q9 + real_t C = 6.*( Concentration - (c[0] + c[2] + c[4] + c[3] + c[7] + c[6]) ); + c[1] = (1./9.) * C; + c[5] = (1./36.) * C; + c[8] = (1./36.) * C; +} + +CudaDeviceFunction void WVelocity() +{ + +} + +CudaDeviceFunction void WVelocityWConcentration() +{ + +//d2q5 +//c[1] = Concentration - c[0] - c[2] - c[3] - c[4]; +//d2q9 + real_t C = 6.*( Concentration - (c[0] + c[2] + c[4] + c[3] + c[7] + c[6]) ); + c[1] = (1./9.) * C; + c[5] = (1./36.) * C; + c[8] = (1./36.) * C; +} + +CudaDeviceFunction void EPressure() +{ + +} + +CudaDeviceFunction void EPressureEConcentration() +{ + +//d2q5 +//c[3] = Concentration - c[0] - c[1] - c[2] - c[4]; +//d2q9 + real_t C = 6.*( Concentration - ( c[0] + c[2] + c[4] + c[1] + c[5] + c[8] ) ); + c[3] = (1./ 9.) * C; + c[7] = (1./36.) * C; + c[6] = (1./36.) * C; +} + +CudaDeviceFunction void NVelocity() +{ + +} + +CudaDeviceFunction void SVelocity() +{ + +} + +CudaDeviceFunction void NSymmetry() +{ + +} + +CudaDeviceFunction void SSymmetry() +{ + +} + +// --- Collision operators --- + +CudaDeviceFunction void CollisionIncompressibleSRT() { + real_t d = getRho(); + vector_t u; + u.x = (( f[8]-f[7]-f[6]+f[5]-f[3]+f[1] )/d + GravitationX/ConcentrationRelaxationRate ); + u.y = ((-f[8]-f[7]+f[6]+f[5]-f[4]+f[2] )/d + GravitationY/ConcentrationRelaxationRate ); + real_t C = getC(); + + real_t phi_eq[9]; + SetIncompressibleEquilibrium(phi_eq, d, u); + for (int i=0; i< 9; i++) { + f[i] = f[i] + RelaxationRate*(phi_eq[i]-f[i]); + } + SetEquilibrium(phi_eq, C, u); + for (int i=0; i< 9; i++) { + c[i] = c[i] + ConcentrationRelaxationRate*(phi_eq[i]-c[i]); + } +} + +CudaDeviceFunction void CollisionMRT_f() { + 1 + R[!selR] = fMRT$Req[!selR] +?> + + real_t ; + real_t Usq=0; + + + real_t dJx = GravitationX*rho ; + real_t dJy = GravitationY*rho ; + + real_t nu = -(S2 + 1)/(6*S2 - 6) ; + real_t K = rho * 12. * nu / h / h; + + Jx = Jx + dJx - K*(Jx+dJx*0.5)/(1. + 0.5 * K) / rho; + Jy = Jy + dJy - K*(Jy+dJy*0.5)/(1. + 0.5 * K) / rho; + + +} + +CudaDeviceFunction void CollisionMRT() { + 1 + R[!selR] = fMRT$Req[!selR] +?> + + real_t ; + real_t Usq=0; + + + real_t dJx = GravitationX*rho ; + real_t dJy = GravitationY*rho ; + + real_t nu = -(S2 + 1)/(6*S2 - 6) ; + real_t K = rho * 12. * nu / h / h; + + Jx = Jx + dJx - K*(Jx+dJx*0.5)/(1. + 0.5 * K) / rho; + Jy = Jy + dJy - K*(Jy+dJy*0.5)/(1. + 0.5 * K) / rho; + + + + 1 + R[!selR] = cMRT$Req[!selR] +?> + real_t C, JCx, JCy; + real_t omega_even = 2.*(2.-ConcentrationRelaxationRate)/(ConcentrationRelaxationRate*(4.*MagicNumber_D-1.)+2.); + real_t SC_even = 1. - omega_even; + 0 + C(R[selR], (R - cMRT$Req)[selR]); + C(R[selR], (R * S)[selR]); +?> +// -- Reaction + if(C < Saturation) { + real_t Q = k * (Saturation - C); + real_t dC = 2*Q; + real_t dh = 2*Q/SolidConcentration; + C += dC; + h += dh; + } + + +} + +CudaDeviceFunction void CollisionIncompressibleMRT() { + 1 + R[!selR] = fMRT_incomp$Req[!selR] +?> + real_t ; + real_t Usq=0; + + real_t dux = GravitationX; + real_t duy = GravitationY; + + real_t nu = -(S2 + 1)/(6*S2 - 6) ; + real_t K = rho * 12. * nu / h / h; + + ux = ux + dux - K*(ux+dux*0.5)/(1. + 0.5 * K) / rho; + uy = uy + duy - K*(uy+duy*0.5)/(1. + 0.5 * K) / rho; + + + + 1 + R[!selR] = cMRT$Req[!selR] +?> + real_t C, JCx, JCy; + real_t omega_even = 2.*(2.-ConcentrationRelaxationRate)/(ConcentrationRelaxationRate*(4.*MagicNumber_D-1.)+2.); + real_t SC_even = 1. - omega_even; + 0 + C(R[selR], (R - cMRT$Req)[selR]); + C(R[selR], (R * S)[selR]); +?> +// -- Reaction + if(C < Saturation) { + real_t Q = k * (Saturation - C); + real_t dC = 2*Q; + real_t dh = 2*Q/SolidConcentration; + C += dC; + h += dh; + } + + +} + +// --- Macroscopic output quantieties --- + +CudaDeviceFunction real_t getRho(){ + return ; +} + +CudaDeviceFunction real_t getP(){ + return ( - 1.)/3.; +} + +CudaDeviceFunction real_t getC(){ + real_t C = ; + +// real_t Q = k * (Saturation - C); +// real_t dC = Q; +// C += dC; + + return C; +} + +CudaDeviceFunction real_t getH(){ + return h; +} + +CudaDeviceFunction vector_t getU(){ + real_t d = ; + vector_t u; + + + + u.x += d*GravitationX*0.5; + u.y += d*GravitationY*0.5; + + real_t nu = -(S2 + 1)/(6*S2 - 6); + + real_t K = d * 12. * nu / h / h; + + u.x = (u.x)/(1. + 0.5 * K) / d; + u.y = (u.y)/(1. + 0.5 * K) / d; + + u.z = 0.0; + return u; +} + +CudaDeviceFunction float2 Color() { + float2 ret; + vector_t u = getU(); + ret.x = sqrt(u.x*u.x + u.y*u.y); + if (NodeType == NODE_Solid){ + ret.y = 0; + } else { + ret.y = 1; + } + return ret; +} \ No newline at end of file diff --git a/models/flow/d2q9_fracture_dissolution/conf.mk b/models/flow/d2q9_fracture_dissolution/conf.mk new file mode 100644 index 000000000..3466dda56 --- /dev/null +++ b/models/flow/d2q9_fracture_dissolution/conf.mk @@ -0,0 +1 @@ +OPT="fields" \ No newline at end of file