diff --git a/models/flow/auto_porous_media/Dynamics.R b/models/flow/auto_porous_media/Dynamics.R new file mode 100644 index 000000000..b42e2646d --- /dev/null +++ b/models/flow/auto_porous_media/Dynamics.R @@ -0,0 +1,143 @@ + + + +IncludedADREModel = TRUE + +Qname = 'DissolutionReaction_ForImplicitSteadyState' +DREs <- ('PHI') +NumberOfODEs = 0 +Params <- c("LinearReactionRate") + +if (Options$d2q9) { + D3 = FALSE +} else { + D3 = TRUE + D3Q19 = FALSE #For advection solver +} +QIntegrator = 'Heun' +#QIntegrator = 'Trapezoid' +source("../models/reaction/d2q9_reaction_diffusion_system/Dynamics.R") + +if (Options$d2q9) { + x = c(0,1,-1); + P = expand.grid(x=0:2,y=0:2,z=0) + U = expand.grid(x,x,0) +} else { + x = c(0,1,-1); + P = expand.grid(x=0:2,y=0:2,z=0:2) + U = expand.grid(x,x,x) +} + + + +f_sel = rep(TRUE,nrow(U)) + +if (Options$d3q19) { + f_sel = rowSums(abs(U)) < 3 +} + +P=P[f_sel,] +U=U[f_sel,] +fname = paste("f",P$x,P$y,P$z,sep="") + +AddDensity( + name = fname, + dx = U[,1], + dy = U[,2], + dz = U[,3], + comment=paste("density",fname), + group="f" +) + + + +#AddDensity( name="fx", group="Force", parameter=FALSE) +#AddDensity( name="fy", group="Force", parameter=FALSE) +#AddDensity( name="fz", group="Force", parameter=FALSE) + + +AddDensity( name="InitialPorosity", group="Brinkman", parameter=TRUE) +#AddDensity( name="InitialPermability", group="Brinkman", parameter=TRUE) + +AddDensity( name="Porosity", group="Brinkman", parameter=TRUE) +AddDensity( name="Permability", group="Brinkman", parameter=TRUE) + +if (Options$GlobalTrapezoid) { + AddDensity( name="GlobalIntegration_1PreviousStep", group="Brinkman", parameter=TRUE) + AddDensity( name="GlobalIntegration_2PreviousStep", group="Brinkman", parameter=TRUE) + +} + + +AddQuantity(name="P",unit="Pa") +AddQuantity(name="U",unit="m/s",vector=T) + +AddQuantity(name="Porosity",unit="1") +AddQuantity(name="Permability",unit="1/m2") +AddQuantity(name="ReactiveFlux",unit="kg/s") +AddQuantity(name="BrinkmanForce",unit="N/m3",vector=TRUE) + + +AddStage(name="InitFromExternal", load.densities=TRUE, save.fields=TRUE) +AddAction(name="InitFromExternalAction", "InitFromExternal") + +AddStage(name="GlobasPorosityDissolutionTimeStep", load.densities=DensityAll$group == "Brinkman", save.fields=Fields$group == "Brinkman") +AddAction(name="GlobasPorosityDissolutionTimeStepAction", "GlobasPorosityDissolutionTimeStep") + + +AddSetting(name="Viscosity", default=0.16666666, comment='Viscosity') +AddSetting(name="Magic", default=3/16, comment='Magic parameter') +AddSetting(name="Velocity", default="0m/s", comment='Inlet velocity', zonal=TRUE) +AddSetting(name="Pressure", default="0Pa", comment='Inlet pressure', zonal=TRUE) + +AddSetting( name="SolidFluidReactionsRate", default="1", comment='Coefficent between fluid and solid mass flux, with dt ratio', zonal=TRUE) +AddSetting( name="ImpliciteReactionIntegration", default="0", comment='Fixed point implicit solution in terms of DPer/Dt.', zonal=FALSE) +AddSetting( name="KarmanKozenyCoefficient", default="1", comment='Permability = KKC * Porosity^3 / ((1-Porosity)^2+1E-8)', zonal=FALSE) + +AddSetting(name="GalileanCorrection",default=1.,comment='Galilean correction term') + +AddSetting(name="ForceX", default=0, comment='Force force X') +AddSetting(name="ForceY", default=0, comment='Force force Y') +AddSetting(name="ForceZ", default=0, comment='Force force Z') + +AddGlobal(name="Flux", comment='Volume flux') +AddGlobal(name="Concentration", comment='Volume flux') + +AddGlobal(name="PressureLoss", comment='pressure loss') +AddGlobal(name="OutletFlux", comment='pressure loss') +AddGlobal(name="InletFlux", comment='pressure loss') + +# AddNodeType(name="NVelocity", group="BOUNDARY") +# AddNodeType(name="SVelocity", group="BOUNDARY") +# AddNodeType(name="NPressure", group="BOUNDARY") +# AddNodeType(name="SPressure", group="BOUNDARY") +# AddNodeType(name="EPressure", group="BOUNDARY") +# AddNodeType(name="EVelocity", group="BOUNDARY") +AddNodeType(name="Wall", group="BOUNDARY") +if (Options$FlowInX || Options$d2q9) { + AddNodeType(name="WDirichlet", group="BOUNDARY") + AddNodeType(name="ENeuman", group="BOUNDARY") +} + +if (Options$FlowInZ) { + AddNodeType(name="IDirichlet", group="BOUNDARY") + AddNodeType(name="ONeuman", group="BOUNDARY") +} + +for (d in rows(DensityAll)) { + if (Options$FlowInX || Options$d2q9) { AddField( name=d$name, dx=-d$dx-1, dy=-d$dy, dz=-d$dz ) } + if (Options$FlowInZ) { AddField( name=d$name, dx=-d$dx, dy=-d$dy, dz=-d$dz-1 ) } +} + +# AddNodeType(name="WVelocity", group="BOUNDARY") + +# AddNodeType(name="W_PHI_DIRICHLET", group="BOUNDARY") +# AddNodeType(name="E_PHI_DIRICHLET", group="BOUNDARY") + +#for (f in fname) AddField(f,dx=0,dy=0,dz=0) # Make f accessible also in present node (not only streamed) + + +AddNodeType(name="Collision", group="COLLISION") + +AddNodeType(name="Inlet", group="OBJECTIVE") +AddNodeType(name="Outlet", group="OBJECTIVE") \ No newline at end of file diff --git a/models/flow/auto_porous_media/Dynamics.c.Rt b/models/flow/auto_porous_media/Dynamics.c.Rt new file mode 100644 index 000000000..cb5a48e0c --- /dev/null +++ b/models/flow/auto_porous_media/Dynamics.c.Rt @@ -0,0 +1,494 @@ + + + + +#define rho0 1 + +real_t AdvectX, AdvectY, AdvectZ; + +CudaDeviceFunction real_t getP() { + if (IamWall) { + return 1/0; + } else { + return (()-1.0)/3.0; + } +} + +CudaDeviceFunction vector_t getU() { + if(IamWall){ + vector_t u; + u.x = 1/0; + u.y = 1/0; + u.z = 1/0; + return u; + } else { + vector_t u; + + const real_t k = (rho0 + 0.5 * Viscosity / Permability); + u.x = (u.x + ForceX/2.)/k; + u.y = (u.y + ForceY/2.)/k; + u.z = (u.z + ForceZ/2.)/k; + return u; + } + +} + + + +CudaDeviceFunction real_t getPermability(){ + return Permability; +} + +CudaDeviceFunction real_t getReactiveFlux(){ + const real_t phi = getPHI(); + + const real_t EffectiveSurfaceArea = 4*Porosity*(1-Porosity); + real_t flux = EffectiveSurfaceArea*LinearReactionRate*(1 - phi); + + if (flux*SolidFluidReactionsRate > (1-InitialPorosity)) { + flux = - (1-InitialPorosity) / SolidFluidReactionsRate; + } + return flux; // negative signe when flux is limited +} + +CudaDeviceFunction real_t getPorosity(){ + if(IamWall){ + return 1/0; + } else { + return Porosity; + } +} + + + +CudaDeviceFunction real_t getInitialPorosity(){ + return InitialPorosity; +} + + +CudaDeviceFunction vector_t getBrinkmanForce(){ + + vector_t force = getU(); + + return force; +} + +CudaDeviceFunction float2 Color() { + float2 ret; + vector_t u = getU(); + ret.x = sqrt(u.x*u.x + u.y*u.y + u.z*u.z); + ret.y = 0; + return ret; +} + +CudaDeviceFunction void BounceBack() +{ + +} + + +CudaDeviceFunction void RunBoundaries() { + + + const real_t k = (rho0 + 0.5 * Viscosity / Permability); + real_t _phi_tmp; + switch (NodeType & NODE_BOUNDARY) { + case NODE_Wall: + BounceBack(); + DoADRE_BounceBack(); + break; + + case NODE_WDirichlet: + { + const real_t BCVelocity = Velocity*k - ForceX / 2; + + } + { + } + break; + + case NODE_ENeuman: + {} + + if (_phi_tmp > 1) { + + } + + break; + + + + case NODE_WDirichlet: + { + const real_t BCVelocity = Velocity*k - ForceX / 2; + + } + { + } + break; + + case NODE_ENeuman: + {} + + if (_phi_tmp > 1) { + + } + + break; + + + + + case NODE_IDirichlet: + { + const real_t BCVelocity = Velocity*k - ForceZ / 2; + + } + { + } + break; + + case NODE_ONeuman: + {} + + if (_phi_tmp > 1) { + + } + + + break; + + + } +} + +CudaDeviceFunction void Run() { + RunBoundaries(); + switch (NodeType & NODE_COLLISION) { + case NODE_COLLISION: + Collision_MTW_RT(); + DoADRE_Run(); + break; + } + +} + +CudaDeviceFunction void SetEquilibrum(real_t p, real_t Vx, real_t Vy, real_t Vz) +{ + + +} + +CudaDeviceFunction void Init() { + + Porosity = InitialPorosity; + Permability = InitialPorosity; + + SetEquilibrum(1.0,0,0,0); + DoADRE_Init(); + + + GlobalIntegration_1PreviousStep = 0; + +} + + +CudaDeviceFunction void InitFromExternal() { + SetEquilibrum(1.0,0,0,0); + DoADRE_InitFromExternal(); +} + + +CudaDeviceFunction void GlobasPorosityDissolutionTimeStep() { + + GlobalIntegration_2PreviousStep = GlobalIntegration_1PreviousStep; + + + InitialPorosity = Porosity; + Permability = KarmanKozenyCoefficient * Porosity*Porosity*Porosity / ((1-Porosity)*(1-Porosity)+1E-8); +} + + +CudaDeviceFunction void Collision_MTW_RT() +{ + + 2] = PV(1) + } + + #if (any(EQ$order < 2)) Omega[EQ$order < 2] = PV(1) + + + meq = EQ$Req + LambdaForce*(EQ_FORCE$feq %*% EQ$mat) + + + m = PV("m",1:ncol(EQ$mat)-1) + cat("real_t",paste(ToC(m),collapse=","),";\n") + C(m, f %*% EQ$mat) + #real_t omega = 1/(3*Viscosity+0.5); #symmetric relaxation + #real_t omega2 = (2.0 - omega)/(1 + 2*(Magic-0.25)*omega); #non-symmetric relaxation + + # Stokes‐Brinkman Flow Simulation Based on 3‐D μ‐CT Images of Porous Rock Using Grayscale Pore Voxel Permeability + # 10.1029/2018WR024179 + # after (Ginzburg et al., 2015) +?> + const real_t Lambda = Magic; + const real_t B = 1. / Permability; + const real_t LambdaPlus = 9. * (4. + B) / ( 4. *( 3.+2.*B*Lambda ) ) * Viscosity; + const real_t LambdaMinus = Lambda / 3. / Viscosity; + + // const real_t LambdaMinus = 3. * Viscosity; + //const real_t LambdaPlus = 3. * Viscosity; + + const real_t omega = 1./( LambdaPlus + 0.5 ); //symmetric relaxation + const real_t omega2 = 1./( LambdaMinus + 0.5 ); //non-symmetric relaxation + + + real_t p, Vx, Vy, Vz, dFX, dFY, dFZ; + const real_t k = 1./(rho0 + 0.5 * Viscosity / Permability); + + + +// Do Objectives/Logs + + switch (NodeType & NODE_OBJECTIVE) { + case NODE_Outlet: + AddToPressureLoss(-p); //only static pressure!!, simplified TODO + break; + case NODE_Inlet: + AddToPressureLoss(p); //only static pressure!!, simplified TODO + break; + } + + AddToFlux(); + +} + + + +// This should go at the end to avoid name conflicts in R + + + + + + CudaDeviceFunction void DoADRE_CalcPhi() + { +/** real_t _q, _phi, _phis; + _phis = ; + + DoADRE_CalcQ(&_phis, &_q); + _phi = _phis + 0.5 * _q; + + for (int i =0; i < 10; i++){ + DoADRE_CalcQ(&_phi, &_q); + _phi = _phis + 0.5 * _q; + //printf("%f %f %f\n", _phi, _phis, _q); + } + + phi[0] = _phi; +**/ + assert(1); + } + + CudaDeviceFunction void DoADRE_CalcQ(const real_t* _phi, real_t* _q) + { + + + if ((_phi[0] > 1) || (_phi[0] < 0)) { + Porosity = InitialPorosity ; + Permability = KarmanKozenyCoefficient * Porosity*Porosity*Porosity / ((1-Porosity)*(1-Porosity)+1E-8); + _q[0] = 0; + GlobalIntegration_1PreviousStep = 0; + return; + }// Limiter + + // Predictor + const real_t EffectiveSurfaceArea = 4*Porosity*(1-Porosity); + real_t flux = EffectiveSurfaceArea*LinearReactionRate*(1 - _phi[0]); + bool limited = 0; + + + if ( ((flux + GlobalIntegration_1PreviousStep) / 2.)*SolidFluidReactionsRate > (1-InitialPorosity)) { + + if ( flux*SolidFluidReactionsRate > (1-InitialPorosity)) { + + flux = (1.-InitialPorosity) / SolidFluidReactionsRate; + limited = 1; + } + + + _q[0] = flux; + GlobalIntegration_1PreviousStep = flux; + if (ImpliciteReactionIntegration == 1) { + + + + if ( limited ){ + Porosity = InitialPorosity + SolidFluidReactionsRate*flux; + } else { + Porosity = InitialPorosity + SolidFluidReactionsRate*(flux + GlobalIntegration_2PreviousStep) / 2.; + } + + + Porosity = InitialPorosity + SolidFluidReactionsRate*flux; + + + + Permability = KarmanKozenyCoefficient * Porosity*Porosity*Porosity / ((1-Porosity)*(1-Porosity)+1E-8); + } + + } + diff --git a/models/flow/auto_porous_media/conf.mk b/models/flow/auto_porous_media/conf.mk new file mode 100644 index 000000000..b8fe69ecb --- /dev/null +++ b/models/flow/auto_porous_media/conf.mk @@ -0,0 +1,3 @@ +OPT="(d3q19+d3q27+d3q15+d2q9)*(TRT)*(GinzburgEqOrd1)*(FlowInX+FlowInZ)*(GlobalTrapezoid)-1" +ADJOINT=0 +TEST=FALSE diff --git a/models/reaction/d2q9_reaction_diffusion_system/Dynamics.R b/models/reaction/d2q9_reaction_diffusion_system/Dynamics.R index 999f8c37e..49f5cd65b 100644 --- a/models/reaction/d2q9_reaction_diffusion_system/Dynamics.R +++ b/models/reaction/d2q9_reaction_diffusion_system/Dynamics.R @@ -1,39 +1,56 @@ -if (Options$AllenCahn) { - Qname = 'Allen-Cahn' - DREs <- ('PHI') - NumberOfODEs = 0 -# NumberOfAdditionalParams = 1 - Params <- c("Lambda") - -} else if (Options$SIR_SimpleLaplace) { - Qname = 'SIR_SimpleLaplace' - #NumberOfDREs = 3 - DREs <- c("S", "I", "R") - NumberOfODEs = 0 - Params <- c("Beta", "Gamma") - -} else if (Options$SIR_ModifiedPeng) { - Qname = 'SIR_ModifiedPeng' - # NumberOfDREs = 1 - # NumberOfODEs = 3+1 - # NumberOfAdditionalParams = 4 - - DREs <- c("W") - ODEs <- c("S", "I", "R", "N") - Params <- c("Beta", "Beta_w", "Gamma") - - -} else if (Options$SimpleDiffusion) { - Qname = 'SimpleDiffusion' - DREs <- ('PHI') - - NumberOfAdditionalParams = 0 - NumberOfODEs = 0 -} +if (!exists('IncludedADREModel')) { + IncludedADREModel = FALSE + D3 = FALSE + if (Options$AllenCahn) { + Qname = 'Allen-Cahn' + DREs <- ('PHI') + NumberOfODEs = 0 + # NumberOfAdditionalParams = 1 + Params <- c("Lambda") + + } else if (Options$SIR_SimpleLaplace) { + Qname = 'SIR_SimpleLaplace' + #NumberOfDREs = 3 + DREs <- c("S", "I", "R") + NumberOfODEs = 0 + Params <- c("Beta", "Gamma") + + } else if (Options$SIR_ModifiedPeng) { + Qname = 'SIR_ModifiedPeng' + # NumberOfDREs = 1 + # NumberOfODEs = 3+1 + # NumberOfAdditionalParams = 4 + + DREs <- c("W") + ODEs <- c("S", "I", "R", "N") + Params <- c("Beta", "Beta_w", "Gamma") + + Extra_Dynamics_C_Header = " + #define EIGEN_DEFAULT_DENSE_INDEX_TYPE int + #include + " + + } else if (Options$SimpleDiffusion) { + Qname = 'SimpleDiffusion' + DREs <- ('PHI') + + NumberOfAdditionalParams = 0 + NumberOfODEs = 0 + + } +} ##END MANUAL CONFIG - +if (!IncludedADREModel){ + if (Options$Heun) { + QIntegrator = 'Heun' + } else if (Options$Euler) { + QIntegrator = 'Euler' + } else { + QIntegrator = 'Trapezoid' + } + } if (exists('DREs')) { NumberOfDREs = length(DREs) @@ -81,20 +98,38 @@ if (NumberOfODEs > 0){ } # declaration of lattice (velocity) directions -x = c(0,1,-1); -P = expand.grid(x=0:2,y=0:2, z=0) -U = expand.grid(x,x,0) + +if (D3) { + x = c(0,1,-1); + P_ADRE = expand.grid(x=0:2,y=0:2, z=0:2) + U_ADRE = expand.grid(x,x,x) + + f_sel = rep(TRUE,nrow(U_ADRE)) + + if (D3Q19) { + f_sel = rowSums(abs(U_ADRE)) < 3 + } + + P_ADRE=P_ADRE[f_sel,] + U_ADRE=U_ADRE[f_sel,] + +} else { + x = c(0,1,-1); + P_ADRE = expand.grid(x=0:2,y=0:2, z=0) + U_ADRE = expand.grid(x,x,0) +} dre_loop(function(i){ # declaration of densities gname = paste("dre",i,sep="_") bname = paste(gname, 'f',sep="_") - fname = paste(bname,P$x,P$y,P$z,sep="") + fname = paste(bname,P_ADRE$x,P_ADRE$y,P_ADRE$z,sep="") AddDensity( name = fname, - dx = U[,1], - dy = U[,2], + dx = U_ADRE[,1], + dy = U_ADRE[,2], + dz = U_ADRE[,3], comment=paste("LB density fields ",gname), group=gname ) @@ -130,18 +165,20 @@ ode_loop(function(i) { AddSetting(name=bname, zonal=TRUE) }) -AddStage(name="InitFromExternal", load.densities=TRUE, save.fields=TRUE) -AddAction(name="InitFromExternalAction", "InitFromExternal") +if (!IncludedADREModel) { + AddStage(name="InitFromExternal", load.densities=TRUE, save.fields=TRUE) + AddAction(name="InitFromExternalAction", "InitFromExternal") -# Boundary things: -AddNodeType(name="Wall", group="BOUNDARY") + AddNodeType(name="Wall", group="BOUNDARY") +} AddNodeType(name="SRT_DF", group="COLLISION") AddNodeType(name="TRT_M", group="COLLISION") + # Inputs: Flow Properties -AddSetting(name="magic_parameter", default=1./6., comment='to control relaxation frequency of even moments in TRT collision kernel') +AddSetting(name="adre_magic_parameter", default=1./6., comment='to control relaxation frequency of even moments in TRT collision kernel') dre_loop(function(i) { bname = paste('Init', DREs[i], sep="_") @@ -160,10 +197,7 @@ if (NumberOfAdditionalParams > 0) { } } -Extra_Dynamics_C_Header = " -#define EIGEN_DEFAULT_DENSE_INDEX_TYPE int -#include -" + # see chapter 10.7.2, eq 10.48, p429 from 'The Lattice Boltzmann Method: Principles and Practice' # by T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, E.M. Viggen diff --git a/models/reaction/d2q9_reaction_diffusion_system/Dynamics.c.Rt b/models/reaction/d2q9_reaction_diffusion_system/Dynamics.c.Rt index ebaa3b55d..0603e159e 100644 --- a/models/reaction/d2q9_reaction_diffusion_system/Dynamics.c.Rt +++ b/models/reaction/d2q9_reaction_diffusion_system/Dynamics.c.Rt @@ -1,9 +1,11 @@ 0){ i = 1 @@ -33,18 +35,28 @@ } } - f_tmp = c(PV(paste('f_tmp[',seq(0,8),']', sep=''))) - qs = c(PV(paste('q[',seq(0,NumberOfODEs + NumberOfDREs),']', sep=''))) + U_ADE = as.matrix(Density[Density$group=='dre_1',c("dx","dy", "dz")]) + adre_f_tmp = c(PV(paste('adre_f_tmp[',seq(0,length(U_ADE[,1])-1),']', sep=''))) + adre_qs = c(PV(paste('q[',seq(0,NumberOfODEs + NumberOfDREs),']', sep=''))) ?> - //this is object-wide variable -real_t f_tmp[9]; +real_t adre_f_tmp[]; real_t q[]; real_t phi[]; +// --- Those functions should not exists if model is included (Actions, Color ETC.) +// --- For such cases Upper model developer is responsible for calling Run_Dispatch + + + +// turn off the advection of Reaction-Diffusion models +#define AdvectX 0 +#define AdvectY 0 +#define AdvectZ 0 + CudaDeviceFunction float2 Color() { float2 ret; ret.x = 0; @@ -52,6 +64,20 @@ CudaDeviceFunction float2 Color() { return ret; } + +CudaDeviceFunction void Init() { + DoADRE_Init(); +} + +CudaDeviceFunction void InitFromExternal() { + DoADRE_InitFromExternal(); +} + +CudaDeviceFunction void Run() { + DoADRE_Run(); +} + + // ------------------------ PARAVIEW OUTPUT BLOCK ------------------------ // Use this functions is only for vtk output. @@ -59,7 +85,7 @@ CudaDeviceFunction float2 Color() { dre_loop( function(i) { ?> CudaDeviceFunction real_t get() { - DispatchCalcPhi(); + DoADRE_DispatchCalcPhi(); return phi[]; } () { ode_loop( function(i) { ?> CudaDeviceFunction real_t get() { - DispatchCalcPhi(); + DoADRE_DispatchCalcPhi(); return phi[]; } () { // ------------------------ END OF PARAVIEW OUTPUT BLOCK ------------------------ -CudaDeviceFunction void Init() { +CudaDeviceFunction void DoADRE_Init() { - - DispatchCalcQ(); - - + DoADRE_DispatchCalcQ(); + - Init_eq(Init_ - 0.5*q[]); - DoADRE_Init_eq(Init_ - 0.5*q[]); DoADRE_Init_eq(Init_); DoADRE_Init_eq(Init_); - DispatchCalcPhi(); + DoADRE_DispatchCalcPhi(); } -CudaDeviceFunction void InitFromExternal() { +CudaDeviceFunction void DoADRE_InitFromExternal() { - - DispatchCalcQ(); - - + DoADRE_DispatchCalcQ(); + - Init_eq(Init__External(0,0) - 0.5*q[]); - DoADRE_Init_eq(Init__External(0,0) - 0.5*q[]); DoADRE_Init_eq(Init__External(0,0)); - DispatchCalcPhi(); + DoADRE_DispatchCalcPhi(); } -CudaDeviceFunction void Run() { - DispatchCalcPhi(); - DispatchCalcQ(); - +CudaDeviceFunction void DoADRE_Run() { + DoADRE_DispatchCalcPhi(); + DoADRE_DispatchCalcQ(); + switch (NodeType & NODE_BOUNDARY) { case NODE_Wall: - BounceBack(); + DoADRE_BounceBack(); return; break; } - + switch (NodeType & NODE_COLLISION) { - case NODE_TRT_M: + + case NODE_SRT_DF: - relax_and_collide_TRT_M(Diffusivity_, q[]); + DoADRE_relax_and_collide_SRT_DF(Diffusivity_, q[]); - break; - case NODE_SRT_DF: + break; + case NODE_TRT_M: default: + - relax_and_collide_SRT_DF(Diffusivity_, q[]); - , q[]); + + DoADRE_relax_and_collide_ADE_CM(Diffusivity_, q[]); + - break; + break; } @@ -205,7 +249,7 @@ CudaDeviceFunction void Run() { @@ -213,47 +257,299 @@ CudaDeviceFunction void Run() { -CudaDeviceFunction void BounceBack() +CudaDeviceFunction void DoADRE_BounceBack() { - + } - CudaDeviceFunction void Init_eq(real_t Init_Phi_Tilde) + CudaDeviceFunction void DoADRE_Init_eq(real_t Init_Phi_Tilde) { } -CudaDeviceFunction void relax_and_collide_SRT_DF(real_t diffusivity, const real_t q) +CudaDeviceFunction void DoADRE_relax_and_collide_TRT(real_t diffusivity, const real_t q) +{ + + + real_t omega_ade = 1.0/(3*diffusivity+0.5); + real_t omega_even = 2.*(2.-omega_ade)/(omega_ade*(4.*adre_magic_parameter-1.)+2.); + real_t tilde_phi = ; + +} + + +CudaDeviceFunction void DoADRE_relax_and_collide_ADE_CM(real_t diffusivity, const real_t q) +//(real_t rho, real_t omega_ade, vector_t u) +{ + const real_t omega_ade = 1.0/(3*diffusivity+0.5); + vector_t u; + + + const real_t rho = 1.; + + real_t Sigma2 = 1/3.; + + //=== THIS IS AUTOMATICALLY GENERATED CODE === + real_t uxuy = u.x*u.y; + real_t ux2 = u.x*u.x; + real_t uy2 = u.y*u.y; + real_t uxuz = u.x*u.z; + real_t uyuz = u.y*u.z; + real_t uz2 = u.z*u.z; + real_t H = h000 + h001 + h002 + h010 + h011 + h012 + h020 + h021 + h022 + h100 + h101 + h102 + h110 + h111 + h112 + h120 + h121 + h122 + h200 + h201 + h202 + h210 + h211 + h212 + h220 + h221 + h222; + real_t temp000 = h000; + real_t temp100 = h100; + real_t temp010 = h010; + real_t temp001 = h001; + real_t temp110 = h110; + real_t temp101 = h101; + real_t temp011 = h011; + real_t temp200 = h200; + real_t temp020 = h020; + real_t temp002 = h002; + real_t temp120 = h120; + real_t temp102 = h102; + real_t temp210 = h210; + real_t temp201 = h201; + real_t temp012 = h012; + real_t temp021 = h021; + real_t temp111 = h111; + real_t temp220 = h220; + real_t temp202 = h202; + real_t temp022 = h022; + real_t temp211 = h211; + real_t temp121 = h121; + real_t temp112 = h112; + real_t temp122 = h122; + real_t temp212 = h212; + real_t temp221 = h221; + real_t temp222 = h222; + //raw moments from density-probability functions + h000 = temp000 + temp001 + temp002 + temp010 + temp011 + temp012 + temp020 + temp021 + temp022 + temp100 + temp101 + temp102 + temp110 + temp111 + temp112 + temp120 + temp121 + temp122 + temp200 + temp201 + temp202 + temp210 + temp211 + temp212 + temp220 + temp221 + temp222; + h100 = temp100 + temp101 + temp102 + temp110 + temp111 + temp112 + temp120 + temp121 + temp122 - temp200 - temp201 - temp202 - temp210 - temp211 - temp212 - temp220 - temp221 - temp222; + h010 = temp010 + temp011 + temp012 - temp020 - temp021 - temp022 + temp110 + temp111 + temp112 - temp120 - temp121 - temp122 + temp210 + temp211 + temp212 - temp220 - temp221 - temp222; + h001 = temp001 - temp002 + temp011 - temp012 + temp021 - temp022 + temp101 - temp102 + temp111 - temp112 + temp121 - temp122 + temp201 - temp202 + temp211 - temp212 + temp221 - temp222; + h110 = temp110 + temp111 + temp112 - temp120 - temp121 - temp122 - temp210 - temp211 - temp212 + temp220 + temp221 + temp222; + h101 = temp101 - temp102 + temp111 - temp112 + temp121 - temp122 - temp201 + temp202 - temp211 + temp212 - temp221 + temp222; + h011 = temp011 - temp012 - temp021 + temp022 + temp111 - temp112 - temp121 + temp122 + temp211 - temp212 - temp221 + temp222; + h200 = temp100 + temp101 + temp102 + temp110 + temp111 + temp112 + temp120 + temp121 + temp122 + temp200 + temp201 + temp202 + temp210 + temp211 + temp212 + temp220 + temp221 + temp222; + h020 = temp010 + temp011 + temp012 + temp020 + temp021 + temp022 + temp110 + temp111 + temp112 + temp120 + temp121 + temp122 + temp210 + temp211 + temp212 + temp220 + temp221 + temp222; + h002 = temp001 + temp002 + temp011 + temp012 + temp021 + temp022 + temp101 + temp102 + temp111 + temp112 + temp121 + temp122 + temp201 + temp202 + temp211 + temp212 + temp221 + temp222; + h120 = temp110 + temp111 + temp112 + temp120 + temp121 + temp122 - temp210 - temp211 - temp212 - temp220 - temp221 - temp222; + h102 = temp101 + temp102 + temp111 + temp112 + temp121 + temp122 - temp201 - temp202 - temp211 - temp212 - temp221 - temp222; + h210 = temp110 + temp111 + temp112 - temp120 - temp121 - temp122 + temp210 + temp211 + temp212 - temp220 - temp221 - temp222; + h201 = temp101 - temp102 + temp111 - temp112 + temp121 - temp122 + temp201 - temp202 + temp211 - temp212 + temp221 - temp222; + h012 = temp011 + temp012 - temp021 - temp022 + temp111 + temp112 - temp121 - temp122 + temp211 + temp212 - temp221 - temp222; + h021 = temp011 - temp012 + temp021 - temp022 + temp111 - temp112 + temp121 - temp122 + temp211 - temp212 + temp221 - temp222; + h111 = temp111 - temp112 - temp121 + temp122 - temp211 + temp212 + temp221 - temp222; + h220 = temp110 + temp111 + temp112 + temp120 + temp121 + temp122 + temp210 + temp211 + temp212 + temp220 + temp221 + temp222; + h202 = temp101 + temp102 + temp111 + temp112 + temp121 + temp122 + temp201 + temp202 + temp211 + temp212 + temp221 + temp222; + h022 = temp011 + temp012 + temp021 + temp022 + temp111 + temp112 + temp121 + temp122 + temp211 + temp212 + temp221 + temp222; + h211 = temp111 - temp112 - temp121 + temp122 + temp211 - temp212 - temp221 + temp222; + h121 = temp111 - temp112 + temp121 - temp122 - temp211 + temp212 - temp221 + temp222; + h112 = temp111 + temp112 - temp121 - temp122 - temp211 - temp212 + temp221 + temp222; + h122 = temp111 + temp112 + temp121 + temp122 - temp211 - temp212 - temp221 - temp222; + h212 = temp111 + temp112 - temp121 - temp122 + temp211 + temp212 - temp221 - temp222; + h221 = temp111 - temp112 + temp121 - temp122 + temp211 - temp212 + temp221 - temp222; + h222 = temp111 + temp112 + temp121 + temp122 + temp211 + temp212 + temp221 + temp222; + //central moments from raw moments + temp000 = h000; + temp100 = -h000*u.x + h100; + temp010 = -h000*u.y + h010; + temp001 = -h000*u.z + h001; + temp110 = h000*uxuy - h010*u.x - h100*u.y + h110; + temp101 = h000*uxuz - h001*u.x - h100*u.z + h101; + temp011 = h000*uyuz - h001*u.y - h010*u.z + h011; + temp200 = h000*ux2 - 2.*h100*u.x + h200; + temp020 = h000*uy2 - 2.*h010*u.y + h020; + temp002 = h000*uz2 - 2.*h001*u.z + h002; + temp120 = -h000*u.x*uy2 + 2.*h010*uxuy - h020*u.x + h100*uy2 - 2.*h110*u.y + h120; + temp102 = -h000*u.x*uz2 + 2.*h001*uxuz - h002*u.x + h100*uz2 - 2.*h101*u.z + h102; + temp210 = -h000*ux2*u.y + h010*ux2 + 2.*h100*uxuy - 2.*h110*u.x - h200*u.y + h210; + temp201 = -h000*ux2*u.z + h001*ux2 + 2.*h100*uxuz - 2.*h101*u.x - h200*u.z + h201; + temp012 = -h000*u.y*uz2 + 2.*h001*uyuz - h002*u.y + h010*uz2 - 2.*h011*u.z + h012; + temp021 = -h000*uy2*u.z + h001*uy2 + 2.*h010*uyuz - 2.*h011*u.y - h020*u.z + h021; + temp111 = -h000*uxuy*u.z + h001*uxuy + h010*uxuz - h011*u.x + h100*uyuz - h101*u.y - h110*u.z + h111; + temp220 = h000*ux2*uy2 - 2.*h010*ux2*u.y + h020*ux2 - 2.*h100*u.x*uy2 + 4.*h110*uxuy - 2.*h120*u.x + h200*uy2 - 2.*h210*u.y + h220; + temp202 = h000*ux2*uz2 - 2.*h001*ux2*u.z + h002*ux2 - 2.*h100*u.x*uz2 + 4.*h101*uxuz - 2.*h102*u.x + h200*uz2 - 2.*h201*u.z + h202; + temp022 = h000*uy2*uz2 - 2.*h001*uy2*u.z + h002*uy2 - 2.*h010*u.y*uz2 + 4.*h011*uyuz - 2.*h012*u.y + h020*uz2 - 2.*h021*u.z + h022; + temp211 = h000*ux2*uyuz - h001*ux2*u.y - h010*ux2*u.z + h011*ux2 - 2.*h100*uxuy*u.z + 2.*h101*uxuy + 2.*h110*uxuz - 2.*h111*u.x + h200*uyuz - h201*u.y - h210*u.z + h211; + temp121 = h000*u.x*uy2*u.z - h001*u.x*uy2 - 2.*h010*uxuy*u.z + 2.*h011*uxuy + h020*uxuz - h021*u.x - h100*uy2*u.z + h101*uy2 + 2.*h110*uyuz - 2.*h111*u.y - h120*u.z + h121; + temp112 = h000*uxuy*uz2 - 2.*h001*uxuy*u.z + h002*uxuy - h010*u.x*uz2 + 2.*h011*uxuz - h012*u.x - h100*u.y*uz2 + 2.*h101*uyuz - h102*u.y + h110*uz2 - 2.*h111*u.z + h112; + temp122 = -h000*u.x*uy2*uz2 + 2.*h001*u.x*uy2*u.z - h002*u.x*uy2 + 2.*h010*uxuy*uz2 - 4.*h011*uxuy*u.z + 2.*h012*uxuy - h020*u.x*uz2 + 2.*h021*uxuz - h022*u.x + h100*uy2*uz2 - 2.*h101*uy2*u.z + h102*uy2 - 2.*h110*u.y*uz2 + 4.*h111*uyuz - 2.*h112*u.y + h120*uz2 - 2.*h121*u.z + h122; + temp212 = -h000*ux2*u.y*uz2 + 2.*h001*ux2*uyuz - h002*ux2*u.y + h010*ux2*uz2 - 2.*h011*ux2*u.z + h012*ux2 + 2.*h100*uxuy*uz2 - 4.*h101*uxuy*u.z + 2.*h102*uxuy - 2.*h110*u.x*uz2 + 4.*h111*uxuz - 2.*h112*u.x - h200*u.y*uz2 + 2.*h201*uyuz - h202*u.y + h210*uz2 - 2.*h211*u.z + h212; + temp221 = -h000*ux2*uy2*u.z + h001*ux2*uy2 + 2.*h010*ux2*uyuz - 2.*h011*ux2*u.y - h020*ux2*u.z + h021*ux2 + 2.*h100*u.x*uy2*u.z - 2.*h101*u.x*uy2 - 4.*h110*uxuy*u.z + 4.*h111*uxuy + 2.*h120*uxuz - 2.*h121*u.x - h200*uy2*u.z + h201*uy2 + 2.*h210*uyuz - 2.*h211*u.y - h220*u.z + h221; + temp222 = h000*ux2*uy2*uz2 - 2.*h001*ux2*uy2*u.z + h002*ux2*uy2 - 2.*h010*ux2*u.y*uz2 + 4.*h011*ux2*uyuz - 2.*h012*ux2*u.y + h020*ux2*uz2 - 2.*h021*ux2*u.z + h022*ux2 - 2.*h100*u.x*uy2*uz2 + 4.*h101*u.x*uy2*u.z - 2.*h102*u.x*uy2 + 4.*h110*uxuy*uz2 - 8.*h111*uxuy*u.z + 4.*h112*uxuy - 2.*h120*u.x*uz2 + 4.*h121*uxuz - 2.*h122*u.x + h200*uy2*uz2 - 2.*h201*uy2*u.z + h202*uy2 - 2.*h210*u.y*uz2 + 4.*h211*uyuz - 2.*h212*u.y + h220*uz2 - 2.*h221*u.z + h222; + //collision in central moments space + //collide + h000 = H; + h100 = -temp100*(omega_ade - 1.); + h010 = -temp010*(omega_ade - 1.); + h001 = -temp001*(omega_ade - 1.); + h110 = 0; + h101 = 0; + h011 = 0; + h200 = H*Sigma2; + h020 = H*Sigma2; + h002 = H*Sigma2; + h120 = 0; + h102 = 0; + h210 = 0; + h201 = 0; + h012 = 0; + h021 = 0; + h111 = 0; + h220 = H*Sigma2*Sigma2; + h202 = H*Sigma2*Sigma2; + h022 = H*Sigma2*Sigma2; + h211 = 0; + h121 = 0; + h112 = 0; + h122 = 0; + h212 = 0; + h221 = 0; + h222 = H*Sigma2*Sigma2*Sigma2; + //back to raw moments + temp000 = h000; + temp100 = h000*u.x + h100; + temp010 = h000*u.y + h010; + temp001 = h000*u.z + h001; + temp110 = h000*uxuy + h010*u.x + h100*u.y + h110; + temp101 = h000*uxuz + h001*u.x + h100*u.z + h101; + temp011 = h000*uyuz + h001*u.y + h010*u.z + h011; + temp200 = h000*ux2 + 2.*h100*u.x + h200; + temp020 = h000*uy2 + 2.*h010*u.y + h020; + temp002 = h000*uz2 + 2.*h001*u.z + h002; + temp120 = h000*u.x*uy2 + 2.*h010*uxuy + h020*u.x + h100*uy2 + 2.*h110*u.y + h120; + temp102 = h000*u.x*uz2 + 2.*h001*uxuz + h002*u.x + h100*uz2 + 2.*h101*u.z + h102; + temp210 = h000*ux2*u.y + h010*ux2 + 2.*h100*uxuy + 2.*h110*u.x + h200*u.y + h210; + temp201 = h000*ux2*u.z + h001*ux2 + 2.*h100*uxuz + 2.*h101*u.x + h200*u.z + h201; + temp012 = h000*u.y*uz2 + 2.*h001*uyuz + h002*u.y + h010*uz2 + 2.*h011*u.z + h012; + temp021 = h000*uy2*u.z + h001*uy2 + 2.*h010*uyuz + 2.*h011*u.y + h020*u.z + h021; + temp111 = h000*uxuy*u.z + h001*uxuy + h010*uxuz + h011*u.x + h100*uyuz + h101*u.y + h110*u.z + h111; + temp220 = h000*ux2*uy2 + 2.*h010*ux2*u.y + h020*ux2 + 2.*h100*u.x*uy2 + 4.*h110*uxuy + 2.*h120*u.x + h200*uy2 + 2.*h210*u.y + h220; + temp202 = h000*ux2*uz2 + 2.*h001*ux2*u.z + h002*ux2 + 2.*h100*u.x*uz2 + 4.*h101*uxuz + 2.*h102*u.x + h200*uz2 + 2.*h201*u.z + h202; + temp022 = h000*uy2*uz2 + 2.*h001*uy2*u.z + h002*uy2 + 2.*h010*u.y*uz2 + 4.*h011*uyuz + 2.*h012*u.y + h020*uz2 + 2.*h021*u.z + h022; + temp211 = h000*ux2*uyuz + h001*ux2*u.y + h010*ux2*u.z + h011*ux2 + 2.*h100*uxuy*u.z + 2.*h101*uxuy + 2.*h110*uxuz + 2.*h111*u.x + h200*uyuz + h201*u.y + h210*u.z + h211; + temp121 = h000*u.x*uy2*u.z + h001*u.x*uy2 + 2.*h010*uxuy*u.z + 2.*h011*uxuy + h020*uxuz + h021*u.x + h100*uy2*u.z + h101*uy2 + 2.*h110*uyuz + 2.*h111*u.y + h120*u.z + h121; + temp112 = h000*uxuy*uz2 + 2.*h001*uxuy*u.z + h002*uxuy + h010*u.x*uz2 + 2.*h011*uxuz + h012*u.x + h100*u.y*uz2 + 2.*h101*uyuz + h102*u.y + h110*uz2 + 2.*h111*u.z + h112; + temp122 = h000*u.x*uy2*uz2 + 2.*h001*u.x*uy2*u.z + h002*u.x*uy2 + 2.*h010*uxuy*uz2 + 4.*h011*uxuy*u.z + 2.*h012*uxuy + h020*u.x*uz2 + 2.*h021*uxuz + h022*u.x + h100*uy2*uz2 + 2.*h101*uy2*u.z + h102*uy2 + 2.*h110*u.y*uz2 + 4.*h111*uyuz + 2.*h112*u.y + h120*uz2 + 2.*h121*u.z + h122; + temp212 = h000*ux2*u.y*uz2 + 2.*h001*ux2*uyuz + h002*ux2*u.y + h010*ux2*uz2 + 2.*h011*ux2*u.z + h012*ux2 + 2.*h100*uxuy*uz2 + 4.*h101*uxuy*u.z + 2.*h102*uxuy + 2.*h110*u.x*uz2 + 4.*h111*uxuz + 2.*h112*u.x + h200*u.y*uz2 + 2.*h201*uyuz + h202*u.y + h210*uz2 + 2.*h211*u.z + h212; + temp221 = h000*ux2*uy2*u.z + h001*ux2*uy2 + 2.*h010*ux2*uyuz + 2.*h011*ux2*u.y + h020*ux2*u.z + h021*ux2 + 2.*h100*u.x*uy2*u.z + 2.*h101*u.x*uy2 + 4.*h110*uxuy*u.z + 4.*h111*uxuy + 2.*h120*uxuz + 2.*h121*u.x + h200*uy2*u.z + h201*uy2 + 2.*h210*uyuz + 2.*h211*u.y + h220*u.z + h221; + temp222 = h000*ux2*uy2*uz2 + 2.*h001*ux2*uy2*u.z + h002*ux2*uy2 + 2.*h010*ux2*u.y*uz2 + 4.*h011*ux2*uyuz + 2.*h012*ux2*u.y + h020*ux2*uz2 + 2.*h021*ux2*u.z + h022*ux2 + 2.*h100*u.x*uy2*uz2 + 4.*h101*u.x*uy2*u.z + 2.*h102*u.x*uy2 + 4.*h110*uxuy*uz2 + 8.*h111*uxuy*u.z + 4.*h112*uxuy + 2.*h120*u.x*uz2 + 4.*h121*uxuz + 2.*h122*u.x + h200*uy2*uz2 + 2.*h201*uy2*u.z + h202*uy2 + 2.*h210*u.y*uz2 + 4.*h211*uyuz + 2.*h212*u.y + h220*uz2 + 2.*h221*u.z + h222; + //back to density-probability functions + h000 = temp000 - temp002 - temp020 + temp022 - temp200 + temp202 + temp220 - temp222; + h100 = 1/2.*temp100 - 1/2.*temp102 - 1/2.*temp120 + 1/2.*temp122 + 1/2.*temp200 - 1/2.*temp202 - 1/2.*temp220 + 1/2.*temp222; + h200 = -1/2.*temp100 + 1/2.*temp102 + 1/2.*temp120 - 1/2.*temp122 + 1/2.*temp200 - 1/2.*temp202 - 1/2.*temp220 + 1/2.*temp222; + h010 = 1/2.*temp010 - 1/2.*temp012 + 1/2.*temp020 - 1/2.*temp022 - 1/2.*temp210 + 1/2.*temp212 - 1/2.*temp220 + 1/2.*temp222; + h110 = 1/4.*temp110 - 1/4.*temp112 + 1/4.*temp120 - 1/4.*temp122 + 1/4.*temp210 - 1/4.*temp212 + 1/4.*temp220 - 1/4.*temp222; + h210 = -1/4.*temp110 + 1/4.*temp112 - 1/4.*temp120 + 1/4.*temp122 + 1/4.*temp210 - 1/4.*temp212 + 1/4.*temp220 - 1/4.*temp222; + h020 = -1/2.*temp010 + 1/2.*temp012 + 1/2.*temp020 - 1/2.*temp022 + 1/2.*temp210 - 1/2.*temp212 - 1/2.*temp220 + 1/2.*temp222; + h120 = -1/4.*temp110 + 1/4.*temp112 + 1/4.*temp120 - 1/4.*temp122 - 1/4.*temp210 + 1/4.*temp212 + 1/4.*temp220 - 1/4.*temp222; + h220 = 1/4.*temp110 - 1/4.*temp112 - 1/4.*temp120 + 1/4.*temp122 - 1/4.*temp210 + 1/4.*temp212 + 1/4.*temp220 - 1/4.*temp222; + h001 = 1/2.*temp001 + 1/2.*temp002 - 1/2.*temp021 - 1/2.*temp022 - 1/2.*temp201 - 1/2.*temp202 + 1/2.*temp221 + 1/2.*temp222; + h101 = 1/4.*temp101 + 1/4.*temp102 - 1/4.*temp121 - 1/4.*temp122 + 1/4.*temp201 + 1/4.*temp202 - 1/4.*temp221 - 1/4.*temp222; + h201 = -1/4.*temp101 - 1/4.*temp102 + 1/4.*temp121 + 1/4.*temp122 + 1/4.*temp201 + 1/4.*temp202 - 1/4.*temp221 - 1/4.*temp222; + h011 = 1/4.*temp011 + 1/4.*temp012 + 1/4.*temp021 + 1/4.*temp022 - 1/4.*temp211 - 1/4.*temp212 - 1/4.*temp221 - 1/4.*temp222; + h111 = 1/8.*temp111 + 1/8.*temp112 + 1/8.*temp121 + 1/8.*temp122 + 1/8.*temp211 + 1/8.*temp212 + 1/8.*temp221 + 1/8.*temp222; + h211 = -1/8.*temp111 - 1/8.*temp112 - 1/8.*temp121 - 1/8.*temp122 + 1/8.*temp211 + 1/8.*temp212 + 1/8.*temp221 + 1/8.*temp222; + h021 = -1/4.*temp011 - 1/4.*temp012 + 1/4.*temp021 + 1/4.*temp022 + 1/4.*temp211 + 1/4.*temp212 - 1/4.*temp221 - 1/4.*temp222; + h121 = -1/8.*temp111 - 1/8.*temp112 + 1/8.*temp121 + 1/8.*temp122 - 1/8.*temp211 - 1/8.*temp212 + 1/8.*temp221 + 1/8.*temp222; + h221 = 1/8.*temp111 + 1/8.*temp112 - 1/8.*temp121 - 1/8.*temp122 - 1/8.*temp211 - 1/8.*temp212 + 1/8.*temp221 + 1/8.*temp222; + h002 = -1/2.*temp001 + 1/2.*temp002 + 1/2.*temp021 - 1/2.*temp022 + 1/2.*temp201 - 1/2.*temp202 - 1/2.*temp221 + 1/2.*temp222; + h102 = -1/4.*temp101 + 1/4.*temp102 + 1/4.*temp121 - 1/4.*temp122 - 1/4.*temp201 + 1/4.*temp202 + 1/4.*temp221 - 1/4.*temp222; + h202 = 1/4.*temp101 - 1/4.*temp102 - 1/4.*temp121 + 1/4.*temp122 - 1/4.*temp201 + 1/4.*temp202 + 1/4.*temp221 - 1/4.*temp222; + h012 = -1/4.*temp011 + 1/4.*temp012 - 1/4.*temp021 + 1/4.*temp022 + 1/4.*temp211 - 1/4.*temp212 + 1/4.*temp221 - 1/4.*temp222; + h112 = -1/8.*temp111 + 1/8.*temp112 - 1/8.*temp121 + 1/8.*temp122 - 1/8.*temp211 + 1/8.*temp212 - 1/8.*temp221 + 1/8.*temp222; + h212 = 1/8.*temp111 - 1/8.*temp112 + 1/8.*temp121 - 1/8.*temp122 - 1/8.*temp211 + 1/8.*temp212 - 1/8.*temp221 + 1/8.*temp222; + h022 = 1/4.*temp011 - 1/4.*temp012 - 1/4.*temp021 + 1/4.*temp022 - 1/4.*temp211 + 1/4.*temp212 + 1/4.*temp221 - 1/4.*temp222; + h122 = 1/8.*temp111 - 1/8.*temp112 - 1/8.*temp121 + 1/8.*temp122 + 1/8.*temp211 - 1/8.*temp212 - 1/8.*temp221 + 1/8.*temp222; + h222 = -1/8.*temp111 + 1/8.*temp112 + 1/8.*temp121 - 1/8.*temp122 + 1/8.*temp211 - 1/8.*temp212 - 1/8.*temp221 + 1/8.*temp222; + + +} + +CudaDeviceFunction void DoADRE_relax_and_collide_SRT_DF(real_t diffusivity, const real_t q) { real_t omega_ade = 1.0/(3*diffusivity+0.5); - real_t tilde_phi = ; + real_t tilde_phi = ; } -CudaDeviceFunction void relax_and_collide_TRT_M(real_t diffusivity, real_t Q) +CudaDeviceFunction void DoADRE_relax_and_collide_TRT_M(real_t diffusivity, real_t Q) { // see eq. 7 from @@ -262,23 +558,23 @@ CudaDeviceFunction void relax_and_collide_TRT_M(real_t diffusivity, real_t Q) // by I. Ginzburg, D. d’Humières, A. Kuzmin, 2010 real_t omega_ade = 1.0/(3*diffusivity+0.5); - real_t omega_even = 2.*(2.-omega_ade)/(omega_ade*(4.*magic_parameter-1.)+2.); + real_t omega_even = 2.*(2.-omega_ade)/(omega_ade*(4.*adre_magic_parameter-1.)+2.); // real_t omega_even = 1.; // **to match notations ** vector_t u; // disable advection - u.x = 0; - u.y = 0; - u.z = 0; - real_t f000 = f_tmp[0]; - real_t f100 = f_tmp[1]; - real_t f200 = f_tmp[2]; - real_t f010 = f_tmp[3]; - real_t f110 = f_tmp[4]; - real_t f210 = f_tmp[5]; - real_t f020 = f_tmp[6]; - real_t f120 = f_tmp[7]; - real_t f220 = f_tmp[8]; + u.x = AdvectX; + u.y = AdvectY; + u.z = AdvectZ; + real_t f000 = adre_f_tmp[0]; + real_t f100 = adre_f_tmp[1]; + real_t f200 = adre_f_tmp[2]; + real_t f010 = adre_f_tmp[3]; + real_t f110 = adre_f_tmp[4]; + real_t f210 = adre_f_tmp[5]; + real_t f020 = adre_f_tmp[6]; + real_t f120 = adre_f_tmp[7]; + real_t f220 = adre_f_tmp[8]; // end of **to match notations ** //=== THIS IS AUTOMATICALLY GENERATED CODE === @@ -343,19 +639,19 @@ CudaDeviceFunction void relax_and_collide_TRT_M(real_t diffusivity, real_t Q) f220 = 1/4.*m_star_110 - 1/4.*m_star_120 - 1/4.*m_star_210 + 1/4.*m_star_220; // **to match notations ** - f_tmp[0] = f000; - f_tmp[1] = f100; - f_tmp[2] = f200; - f_tmp[3] = f010; - f_tmp[4] = f110; - f_tmp[5] = f210; - f_tmp[6] = f020; - f_tmp[7] = f120; - f_tmp[8] = f220; + adre_f_tmp[0] = f000; + adre_f_tmp[1] = f100; + adre_f_tmp[2] = f200; + adre_f_tmp[3] = f010; + adre_f_tmp[4] = f110; + adre_f_tmp[5] = f210; + adre_f_tmp[6] = f020; + adre_f_tmp[7] = f120; + adre_f_tmp[8] = f220; // end of **to match notations ** } -CudaDeviceFunction void DispatchCalcPhi() +CudaDeviceFunction void DoADRE_DispatchCalcPhi() { switch (NodeType & NODE_BOUNDARY) { @@ -382,9 +678,24 @@ CudaDeviceFunction void DispatchCalcPhi() break; } - CalcPhi(); + + DoADRE_CalcPhi(); + + + } -CudaDeviceFunction void DispatchCalcQ() +CudaDeviceFunction void DoADRE_DispatchCalcQ() { switch (NodeType & NODE_BOUNDARY) { case NODE_Wall: @@ -408,232 +719,55 @@ CudaDeviceFunction void DispatchCalcQ() break; } - CalcQ(); -} - - - - - - - - CudaDeviceFunction void CalcPhi() - { - real_t tilde_phi = ; - - const real_t dt = 1.; - const real_t x0 = cbrt(3.) ; - const real_t x1 = 1/(dt*Lambda) ; - const real_t x2 = dt*Lambda - 2 ; - const real_t x3 = cbrt(-9*tilde_phi*x1 + sqrt(3.)*sqrt((27*pow(tilde_phi, 2) - x1*pow(x2, 3))/(pow(dt, 2)*pow(Lambda, 2)))) ; - phi[0] = -1.0/3.0*x0*(x0*x1*x2 + pow(x3, 2))/x3 ; - } - CudaDeviceFunction void CalcQ() - { - q[0] = Lambda*phi[0]*(1 - pow(phi[0], 2)) ; - } - - - - - - - - CudaDeviceFunction void CalcPhi() - { - const real_t x0 = ; - const real_t x1 = ; - const real_t x2 = ; - - // x0 = S^\star - // x1 = I^\star - // x2 = R^\star - const real_t x3 = Beta; - const real_t x4 = Gamma; - const real_t x5 = 1; - // Opers0 = 410 - const real_t x6 = x4*x5 ; // 1 - const real_t x7 = x3*x5 ; // 1 - const real_t x8 = x0*x7 ; // 1 - const real_t x9 = x1*x7 ; // 1 - const real_t x10 = pow(x4, 2) ; // 1 - const real_t x11 = pow(x5, 2) ; // 1 - const real_t x12 = 2*x4 ; // 1 - const real_t x13 = x11*x12*x3 ; // 2 - const real_t x14 = x11*pow(x3, 2) ; // 2 - const real_t x15 = sqrt(pow(x0, 2)*x14 + 2*x0*x1*x14 - x0*x13 + pow(x1, 2)*x14 + x1*x13 + x10*x11 + 4*x6 - 4*x8 + 4*x9 + 4) ; // 24 - const real_t x16 = x6 + 2 ; // 1 - const real_t x17 = x15 + x16 ; // 1 - const real_t x18 = x8 + x9 ; // 1 - const real_t x19 = 1.0/x3 ; // 1 - const real_t x20 = x19/x5 ; // 1 - const real_t x21 = (1.0/2.0)*x20 ; // 1 - const real_t x22 = 1.0/x16 ; // 1 - const real_t x23 = x20*x22 ; // 1 - const real_t x24 = x15*x4 ; // 1 - const real_t x25 = x2*x3 ; // 1 - const real_t x26 = x3*x6 ; // 1 - const real_t x27 = x0*x26 + x1*x26 - x10*x5 - x12 + 2*x25*x6 + 4*x25 ; // 11 - const real_t x28 = (1.0/2.0)*x19*x22 ; // 2 - const real_t Sp = x21*(x17 + x18) ; // 2 - const real_t Ip = x23*(-x17 + x8 + x9) ; // 3 - const real_t Rp = x28*(-x24 + x27) ; // 2 - const real_t Sm = x21*(-x15 + x16 + x18) ; // 3 - const real_t Im = x23*(x15 + x18 - x6 - 2) ; // 4 - const real_t Rm = x28*(x24 + x27) ; // 2 - // Opers = 104 - - - const real_t eps = -1E-12; - if ( Sp >=eps && Ip >=eps && Rp >=eps ) { - phi[0] = Sp < 0 ? 0 : Sp; - phi[1] = Ip < 0 ? 0 : Ip; - phi[2] = Rp < 0 ? 0 : Rp; - } else if ( Sm >=eps&& Im >=eps && Rm >=eps ) { - phi[0] = Sm < 0 ? 0 : Sm; - phi[1] = Im < 0 ? 0 : Im; - phi[2] = Rm < 0 ? 0 : Rm; - } else { - printf("x3 x4 %f %f \n", x3, x4); - printf("s %e %e %e \n", x0, x1, x2); - printf("p %e %e %e \n", Sp, Ip, Rp); - printf("m %e %e %e \n", Sm, Im, Rm); - assert(0); - } - - - } - CudaDeviceFunction void CalcQ() - { - const real_t s = phi[0] ; - const real_t i = phi[1] ; - const real_t r = phi[2] ; - q[0] = -Beta*s*i ; - q[1] = Beta*s*i - Gamma*i; - q[2] = Gamma*i ; - - // printf("Phi %f %f %f \n", phi[0], phi[1], phi[2]); - // printf("Q %f %f %f \n", q[0], q[1], q[2]); - + + DoADRE_CalcQ(phi,q); + - - - - CudaDeviceFunction void CalcPhi() - { - - const real_t x0 = ; - const real_t x1 = (0,0); - const real_t x2 = (0,0); - const real_t x3 = (0,0); - const real_t N = (0,0); - - phi[4] = N; - - - Eigen::Matrix3d Jacobian; - - Eigen::Vector3d Xn, dX, X0, F; - Xn(0) = x0/N; - Xn(1) = x1/N; - Xn(2) = x2/N; - - X0 = Xn; - - Jacobian(0,0) = -1.0/2.0*Beta_w - 1 ; - Jacobian(1,0) = 0 ; - Jacobian(2,0) = (1.0/2.0)*Beta_w ; - Jacobian(2,1) = 0 ; - Jacobian(2,2) = -1.0/2.0*Gamma - 1 ; - - for (int i=0; i < 20; i++) { - - F(0) = (1.0/2.0)*Beta_w*(-Xn(0) + Xn(2)) - Xn(0) + X0(0) ; - F(1) = -1.0/2.0*Beta*Xn(0)*Xn(1) - Xn(1) + X0(1) ; - F(2) = (1.0/2.0)*Beta*Xn(0)*Xn(1) - 1.0/2.0*Xn(2)*Gamma - Xn(2) + X0(2) ; - if (F.norm() < 1E-5) { - break; - } - - Jacobian(0,1) = -1.0/2.0*Beta*Xn(1) ; - Jacobian(0,2) = (1.0/2.0)*Beta*Xn(1) ; - Jacobian(1,1) = -1.0/2.0*Beta*Xn(0) - 1 ; - Jacobian(1,2) = (1.0/2.0)*Beta*Xn(0) ; + if (QIntegrator == 'Heun') { + ?> + real_t q_tmp[]; + real_t phi_euler[]; + DoADRE_CalcQ(phi,q_tmp); + + phi_euler[] = phi[] + q_tmp[]; + + phi_euler[] = phi[] + q_tmp[]; + + DoADRE_CalcQ(phi_euler,q); + + q[] = ( q[] + q_tmp[] ) / 2; + + q[] = ( q[] + q_tmp[] ) / 2; + +} - CudaDeviceFunction void CalcQ() - { - - - const real_t W = phi[0]; - const real_t S = phi[1]; - const real_t I = phi[2]; - const real_t R = phi[3]; - - const real_t N = phi[4]; - - - q[0] = Beta_w*(I-W); - q[1] = -Beta*S*W / N; - q[2] = (Beta*S*W/N - Gamma*I); - q[3] = Gamma*I; - q[4] = 0; //N - //printf("CalcQ Phi %lf %lf %lf %lf\n", phi[0], phi[1], phi[2], phi[3]); - //printf("CalcQ Q %lf %lf %lf %lf\n", q[0], q[1], q[2], q[3]); - } - - - - - CudaDeviceFunction void CalcPhi() - { - phi[0] = ; - } - CudaDeviceFunction void CalcQ() - { - q[0] = 0; - } - + diff --git a/models/reaction/d2q9_reaction_diffusion_system/ReactionTerms.c.Rt b/models/reaction/d2q9_reaction_diffusion_system/ReactionTerms.c.Rt new file mode 100644 index 000000000..5963ba67c --- /dev/null +++ b/models/reaction/d2q9_reaction_diffusion_system/ReactionTerms.c.Rt @@ -0,0 +1,241 @@ + + + + CudaDeviceFunction void DoADRE_CalcPhi() + { + real_t tilde_phi = ; + + const real_t dt = 1.; + const real_t x0 = cbrt(3.) ; + const real_t x1 = 1/(dt*Lambda) ; + const real_t x2 = dt*Lambda - 2 ; + const real_t x3 = cbrt(-9*tilde_phi*x1 + sqrt(3.)*sqrt((27*pow(tilde_phi, 2) - x1*pow(x2, 3))/(pow(dt, 2)*pow(Lambda, 2)))) ; + phi[0] = -1.0/3.0*x0*(x0*x1*x2 + pow(x3, 2))/x3 ; + } + CudaDeviceFunction void DoADRE_CalcQ(const real_t* _phi, real_t* _q) + { + _q[0] = Lambda*_phi[0]*(1 - pow(_phi[0], 2)) ; + } + + + + + + + + CudaDeviceFunction void DoADRE_CalcPhi() + { + const real_t x0 = ; + const real_t x1 = ; + const real_t x2 = ; + + // x0 = S^\star + // x1 = I^\star + // x2 = R^\star + const real_t x3 = Beta; + const real_t x4 = Gamma; + const real_t x5 = 1; + // Opers0 = 410 + const real_t x6 = x4*x5 ; // 1 + const real_t x7 = x3*x5 ; // 1 + const real_t x8 = x0*x7 ; // 1 + const real_t x9 = x1*x7 ; // 1 + const real_t x10 = pow(x4, 2) ; // 1 + const real_t x11 = pow(x5, 2) ; // 1 + const real_t x12 = 2*x4 ; // 1 + const real_t x13 = x11*x12*x3 ; // 2 + const real_t x14 = x11*pow(x3, 2) ; // 2 + const real_t x15 = sqrt(pow(x0, 2)*x14 + 2*x0*x1*x14 - x0*x13 + pow(x1, 2)*x14 + x1*x13 + x10*x11 + 4*x6 - 4*x8 + 4*x9 + 4) ; // 24 + const real_t x16 = x6 + 2 ; // 1 + const real_t x17 = x15 + x16 ; // 1 + const real_t x18 = x8 + x9 ; // 1 + const real_t x19 = 1.0/x3 ; // 1 + const real_t x20 = x19/x5 ; // 1 + const real_t x21 = (1.0/2.0)*x20 ; // 1 + const real_t x22 = 1.0/x16 ; // 1 + const real_t x23 = x20*x22 ; // 1 + const real_t x24 = x15*x4 ; // 1 + const real_t x25 = x2*x3 ; // 1 + const real_t x26 = x3*x6 ; // 1 + const real_t x27 = x0*x26 + x1*x26 - x10*x5 - x12 + 2*x25*x6 + 4*x25 ; // 11 + const real_t x28 = (1.0/2.0)*x19*x22 ; // 2 + const real_t Sp = x21*(x17 + x18) ; // 2 + const real_t Ip = x23*(-x17 + x8 + x9) ; // 3 + const real_t Rp = x28*(-x24 + x27) ; // 2 + const real_t Sm = x21*(-x15 + x16 + x18) ; // 3 + const real_t Im = x23*(x15 + x18 - x6 - 2) ; // 4 + const real_t Rm = x28*(x24 + x27) ; // 2 + // Opers = 104 + + + const real_t eps = -1E-12; + if ( Sp >=eps && Ip >=eps && Rp >=eps ) { + phi[0] = Sp < 0 ? 0 : Sp; + phi[1] = Ip < 0 ? 0 : Ip; + phi[2] = Rp < 0 ? 0 : Rp; + } else if ( Sm >=eps&& Im >=eps && Rm >=eps ) { + phi[0] = Sm < 0 ? 0 : Sm; + phi[1] = Im < 0 ? 0 : Im; + phi[2] = Rm < 0 ? 0 : Rm; + } else { + printf("x3 x4 %f %f \n", x3, x4); + printf("s %e %e %e \n", x0, x1, x2); + printf("p %e %e %e \n", Sp, Ip, Rp); + printf("m %e %e %e \n", Sm, Im, Rm); + assert(0); + } + + + } + CudaDeviceFunction void DoADRE_CalcQ(const real_t* _phi, real_t* _q) + { + const real_t s = _phi[0] ; + const real_t i = _phi[1] ; + const real_t r = _phi[2] ; + _q[0] = -Beta*s*i ; + _q[1] = Beta*s*i - Gamma*i; + _q[2] = Gamma*i ; + + // printf("Phi %f %f %f \n", phi[0], phi[1], phi[2]); + // printf("Q %f %f %f \n", q[0], q[1], q[2]); + + + } + + + + + + + CudaDeviceFunction void DoADRE_CalcPhi() + { + + const real_t x0 = ; + const real_t x1 = (0,0); + const real_t x2 = (0,0); + const real_t x3 = (0,0); + const real_t N = (0,0); + + phi[4] = N; + + + Eigen::Matrix3d Jacobian; + + Eigen::Vector3d Xn, dX, X0, F; + Xn(0) = x0/N; + Xn(1) = x1/N; + Xn(2) = x2/N; + + X0 = Xn; + + Jacobian(0,0) = -1.0/2.0*Beta_w - 1 ; + Jacobian(1,0) = 0 ; + Jacobian(2,0) = (1.0/2.0)*Beta_w ; + Jacobian(2,1) = 0 ; + Jacobian(2,2) = -1.0/2.0*Gamma - 1 ; + + for (int i=0; i < 20; i++) { + + F(0) = (1.0/2.0)*Beta_w*(-Xn(0) + Xn(2)) - Xn(0) + X0(0) ; + F(1) = -1.0/2.0*Beta*Xn(0)*Xn(1) - Xn(1) + X0(1) ; + F(2) = (1.0/2.0)*Beta*Xn(0)*Xn(1) - 1.0/2.0*Xn(2)*Gamma - Xn(2) + X0(2) ; + if (F.norm() < 1E-5) { + break; + } + + Jacobian(0,1) = -1.0/2.0*Beta*Xn(1) ; + Jacobian(0,2) = (1.0/2.0)*Beta*Xn(1) ; + Jacobian(1,1) = -1.0/2.0*Beta*Xn(0) - 1 ; + Jacobian(1,2) = (1.0/2.0)*Beta*Xn(0) ; + + dX = Jacobian.ldlt().solve(-F); + + Xn = Xn + dX; + + } + Xn = N*Xn; + for(int k=0; k < 3; k++){ + Xn(k) = fabs(Xn(k)); + } + phi[0] = Xn(0); + phi[1] = Xn(1); + phi[2] = Xn(2); + phi[3] = Xn(2)*Gamma/2. + x3; + //printf("CalcPhi PhiTilde %lf %lf %lf %lf\n", x0, x1, x2, x3); + //printf("CalcPhi Phi %lf %lf %lf %lf\n", phi[0], phi[1], phi[2], phi[3]); + + return; + + } + + + + + + CudaDeviceFunction void DoADRE_CalcQ(const real_t* _phi, real_t* _q) + { + _q[0] + + + const real_t W = _phi[0]; + const real_t S = _phi[1]; + const real_t I = _phi[2]; + const real_t R = _phi[3]; + + const real_t N = _phi[4]; + + + _q[0] = Beta_w*(I-W); + _q[1] = -Beta*S*W / N; + _q[2] = (Beta*S*W/N - Gamma*I); + _q[3] = Gamma*I; + _q[4] = 0; //N + //printf("CalcQ Phi %lf %lf %lf %lf\n", phi[0], phi[1], phi[2], phi[3]); + //printf("CalcQ Q %lf %lf %lf %lf\n", q[0], q[1], q[2], q[3]); + } + + + + + CudaDeviceFunction void DoADRE_CalcPhi() + { + phi[0] = ; + } + CudaDeviceFunction void DoADRE_CalcQ(const real_t* _phi, real_t* _q) + { + _q[0] = 0; + } + + + + + + CudaDeviceFunction void DoADRE_CalcPhi() + { + assert(0); // not implemented + } + CudaDeviceFunction void DoADRE_CalcQ(const real_t* _phi, real_t* _q) + { + _q[0] = LinearReactionRate*_phi[0]; + } + diff --git a/models/reaction/d2q9_reaction_diffusion_system/TestWSIRReversingStep.cpp b/models/reaction/d2q9_reaction_diffusion_system/TestWSIRReversingStep.cpp new file mode 100644 index 000000000..b19bfecb2 --- /dev/null +++ b/models/reaction/d2q9_reaction_diffusion_system/TestWSIRReversingStep.cpp @@ -0,0 +1,634 @@ +#define CudaDeviceFunction +#define real_t double +#include +#include +#include +//x4 x5 x6 x7 +// 4.269327243786159798e+00 -2.705207483582461315e-13 4.313636855125160707e+00 1.000000000000000000e+01 8 +// p -2.047304e+05 4.115177e+07 -4.115166e+07 -1.028791e+02 +// m 4.269548e+00 -1.338658e-09 4.313626e+00 4.313648e-05 + +// const real_t w0 = 0.1;//4.269327243786159798e+00; +// const real_t s0 = 10;// -2.705207483582461315e-13; +// const real_t i0 = 0.1;//4.313636855125160707e+00; +// const real_t r0 = 1e-15;///3.235241e-05; +// const real_t n0 = s0+r0+i0;//4.313669207536073635e+00; + +// const real_t C_1 = 1;//4.213999999999999972e-05; +// const real_t C_2 = 1;//1.000000000000000021e-02; +// const real_t C_3 = 1;//5.000000000000000409e-06; + + +const real_t w0 = 4.269327243786159798e+00; +const real_t s0 = -2.705207483582461315e-13; +const real_t i0 = 4.313636855125160707e+00; +const real_t r0 = 3.235241e-05; +const real_t n0 = 4.313669207536073635e+00; + +const real_t C_1 = 4.213999999999999972e-05; +const real_t C_2 = 1.000000000000000021e-02; +const real_t C_3 = 5.000000000000000409e-06; + +CudaDeviceFunction void localCalcQ(real_t* localPhi, real_t* localQ, real_t N) +{ + + const real_t Beta = C_1; + const real_t Beta_w = C_2; + const real_t Gamma = C_3; + + + + const real_t W = localPhi[0]; + const real_t S = localPhi[1]; + const real_t I = localPhi[2]; + //const real_t R = localPhi[3]; + + //const real_t N = localPhi[4]; + + localQ[0] = Beta_w*(I-W); + localQ[1] = -Beta*S*W / N; + localQ[2] = (Beta*S*W/N - Gamma*I); + //localQ[3] = Gamma*I; + //localQ[4] = 0; //N + //printf("Phi %f %f %f %f\n", localPhi[0], localPhi[1], localPhi[2], localPhi[3]); + //printf("Q %f %f %f %f\n", q[0], q[1], q[2], q[3]); +} + +#include +#include + +using namespace std; +using namespace Eigen; + +CudaDeviceFunction int main() +{ + + const real_t x0 = w0; + const real_t x1 = s0; + const real_t x2 = i0; + real_t X1[3], X2[3], f1[3], f2[3]; + + // X1[0] = w0; + // X1[1] = s0; + // X1[2] = i0; + + // localCalcQ(X1, f1, n0); + // X1[0] = X1[0] - f1[0]/2.; + // X1[1] = X1[1] - f1[1]/2.; + // X1[2] = X1[2] - f1[2]/2.; + + // const real_t x0 = X1[0]; + // const real_t x1 = X1[1]; + // const real_t x2 = X1[2]; + + + const real_t Beta = C_1; + const real_t Beta_w = C_2; + const real_t gamma = C_3; + const real_t N = n0; + + Matrix3d Jacobian; + + Vector3d X, dX, X0, F; + X(0) = x0/N; + X(1) = x1/N; + X(2) = x2/N; + + X0 = X; + cout << "Newton Init " << X.transpose() << endl; + + for (int i=0; i < 100; i++) { + Jacobian(0,0) = -1.0/2.0*Beta_w - 1 ; + Jacobian(0,1) = -1.0/2.0*Beta*X(1) ; + Jacobian(0,2) = (1.0/2.0)*Beta*X(1) ; + // Jacobian(0,3) = 0 ; + Jacobian(1,0) = 0 ; + Jacobian(1,1) = -1.0/2.0*Beta*X(0) - 1 ; + Jacobian(1,2) = (1.0/2.0)*Beta*X(0) ; + // Jacobian(1,3) = 0 ; + Jacobian(2,0) = (1.0/2.0)*Beta_w ; + Jacobian(2,1) = 0 ; + Jacobian(2,2) = -1.0/2.0*gamma - 1 ; + // Jacobian(2,3) = (1.0/2.0)*gamma ; + // Jacobian(3,0) = 0 ; + // Jacobian(3,1) = 0 ; + // Jacobian(3,2) = 0 ; + // Jacobian(3,3) = -1 ; + + F(0) = (1.0/2.0)*Beta_w*(-X(0) + X(2)) - X(0) + X0(0) ; + F(1) = -1.0/2.0*Beta*X(0)*X(1) - X(1) + X0(1) ; + F(2) = (1.0/2.0)*Beta*X(0)*X(1) - 1.0/2.0*X(2)*gamma - X(2) + X0(2) ; + // F(3) = (1.0/2.0)*X(2)*gamma - X(3) + X0(3) ; + + dX = Jacobian.householderQr().solve(-F); + // cout << "Newton Error " << F.norm() << endl; + X = X + dX; + if (F.norm() < 1E-16) { + cout << "Converged in "<< i << endl; + break; + } + cout << "Newton Step " << X.transpose() << endl; + + // cout << "Newton dX " << dX.transpose() << endl; + + // cout << "Newton Step " << X.transpose() << endl; + } + cout << "Newton Error " << F.norm() << endl; + cout << "Newton Result " << N*X.transpose() << endl<< endl<< endl; + // x4 x5 x6 x7 4.313669e+00 4.214000e-05 1.000000e-02 5.000000e-06 + // s 4.269327e+00 -2.705207e-13 4.313637e+00 3.235241e-05 + + + // Try Secant + + + + + + + + // END FAKE INTI + + X1[0] = fabs(x0); + X1[1] = fabs(x1); + X1[2] = fabs(x2); + + localCalcQ(X1, f1, n0); + + X2[0] = fabs(X1[0] + f1[0]/2.); + X2[1] = fabs(X1[1] + f1[1]/2.); + X2[2] = fabs(X1[2] + f1[2]/2.); + + printf("Secant 0: X1 %e %e %e \n", X1[0], X1[1], X1[2]); + printf("Secant 0: X2 %e %e %e \n", X2[0], X2[1], X2[2]); + + for(int i =0; i < 20; i++) { + localCalcQ(X1, f1, n0); + // printf("Newton Q1 %e %e %e %e \n", f1[0], f1[1], f1[2]); + + f1[0] = x0 + f1[0]/2. - X1[0]; + f1[1] = x1 + f1[1]/2. - X1[1]; + f1[2] = x2 + f1[2]/2. - X1[2]; + + + real_t epsilon_step = 0; + for(int k=0; k < 3; k++){ + if ( fabs(f1[k]) > epsilon_step ) { + epsilon_step = fabs(f1[k]); + } + } + if (epsilon_step < 1E-12){ + + for(int k=0; k < 3; k++){ + X1[k] = fabs(X1[k]); + } + printf("Secant result: %e %e %e EPS: %e\n", X1[0], X1[1], X1[2], epsilon_step); + printf("CONVERGED\n\n"); + break; + } + + + // printf("Newton F1 %e %e %e \n", f1[0], f1[1], f1[2]); + + localCalcQ(X2, f2, n0); + // printf("Newton Q2 %e %e %e \n", f2[0], f2[1], f2[2]); + f2[0] = x0 + f2[0]/2. - X2[0]; + f2[1] = x1 + f2[1]/2. - X2[1]; + f2[2] = x2 + f2[2]/2. - X2[2]; + + // printf("Newton F2 %e %e %e \n", f2[0], f2[1], f2[2]); + + //printf("Iter %e %e %e %e \n", f1[0], f1[1], f1[2], f1[3]); + // printf("Newton %e %e %e %e \n", X2[0]-x0, X2[1]-x1, X2[2]-x2, X2[3]-x3); + + for(int k=0; k < 3; k++){ + const real_t tt = X1[k]; + X1[k] = X1[k] - f1[k] * (X1[k] - X2[k]) / (f1[k] - f2[k]); + X2[k] = tt; + } + + printf("Secant X1 %e %e %e EPS: %e\n", X1[0], X1[1], X1[2], epsilon_step); + + printf("STEP\n\n"); + + + } + + + //printf("x4 x5 x6 x7 %e %e %e %e\n", x4, x5, x6, x7); + //printf("m %e %e %e %e \n", X2[0]-x0, X2[1]-x1, X2[2]-x2, X2[3]-x3); + + + + return 0; + +} + + +// CudaDeviceFunction void CalcQ() +// { + +// const real_t Beta = C_1; +// const real_t Beta_w = C_2; +// const real_t Gamma = C_3; + + + +// const real_t W = phi[0]; +// const real_t S = phi[1]; +// const real_t I = phi[2]; +// const real_t R = phi[3]; + +// const real_t N = phi[4]; + + +// q[0] = Beta_w*(I-W); +// q[1] = -Beta*S*W / N; +// q[2] = (Beta*S*W/N - Gamma*I); +// q[3] = Gamma*I; +// q[4] = 0; //N +// //printf("Phi %f %f %f %f\n", phi[0], phi[1], phi[2], phi[3]); +// //printf("Q %f %f %f %f\n", q[0], q[1], q[2], q[3]); +// } + + +// x4 x5 x6 x7 4.313669e+00 4.214000e-05 1.000000e-02 5.000000e-06 +// s 4.269327e+00 -2.705207e-13 4.313637e+00 1.000000e+01 +// p -2.047304e+05 4.115177e+07 -4.115166e+07 -1.028791e+02 +// m 4.269548e+00 -1.338658e-09 4.313626e+00 4.313648e-05 + + +// x4 x5 x6 x7 4.313669e+00 4.214000e-05 1.000000e-02 5.000000e-06 +// s 4.269327e+00 -2.705207e-13 4.313637e+00 1.000000e+01 +// p -2.047304e+05 4.115177e+07 -4.115166e+07 -1.028791e+02 +// m 4.269547e+00 2.459523e-09 4.313626e+00 4.313648e-05 + + +// // phi[4] = x4; // N +// // x0 = W^\star +// // x1 = S^\star +// // x2 = I^\star +// // x3 = R^\star + +// // x4 = N +// //const real_t x4 = x1 + x2 + x3; + +// // x5 = Beta +// const real_t x5 = C_1; + +// // x6 = Beta_w +// const real_t x6 = C_2; + +// // x7 = gamma +// const real_t x7 = C_3; + + +// const real_t x8 = x0*x5 ; // 1 +// const real_t x9 = 2*x8 ; // 1 +// const real_t x10 = x7*x8 ; // 1 +// const real_t x11 = pow(x4, 2) ; // 1 +// const real_t x12 = 16*x11 ; // 1 +// const real_t x13 = x4*x8 ; // 1 +// const real_t x14 = x5*x6 ; // 1 +// const real_t x15 = x1*x14 ; // 1 +// const real_t x16 = 8*x4 ; // 1 +// const real_t x17 = x14*x2 ; // 1 +// const real_t x18 = x6*x7 ; // 1 +// const real_t x19 = x18*x4 ; // 1 +// const real_t x20 = 4*x4 ; // 1 +// const real_t x21 = x18*x5 ; // 1 +// const real_t x22 = x1*x21 ; // 1 +// const real_t x23 = x2*x21 ; // 1 +// const real_t x24 = pow(x5, 2) ; // 1 +// const real_t x25 = pow(x0, 2)*x24 ; // 2 +// const real_t x26 = 4*x25 ; // 1 +// const real_t x27 = pow(x6, 2) ; // 1 +// const real_t x28 = 4*x11 ; // 1 +// const real_t x29 = x27*x28 ; // 1 +// const real_t x30 = pow(x7, 2) ; // 1 +// const real_t x31 = x28*x30 ; // 1 +// const real_t x32 = x0*x24 ; // 1 +// const real_t x33 = x1*x32 ; // 1 +// const real_t x34 = 4*x6 ; // 1 +// const real_t x35 = x2*x32 ; // 1 +// const real_t x36 = x30*x8 ; // 1 +// const real_t x37 = x27*x5 ; // 1 +// const real_t x38 = x20*x37 ; // 1 +// const real_t x39 = 2*x6 ; // 1 +// const real_t x40 = x39*x7 ; // 1 +// const real_t x41 = x39*x4 ; // 1 +// const real_t x42 = 2*x7 ; // 1 +// const real_t x43 = x4*x42 ; // 1 +// const real_t x44 = x37*x43 ; // 1 +// const real_t x45 = x24*x27 ; // 1 +// const real_t x46 = sqrt(pow(x1, 2)*x45 + 2*x1*x2*x45 - x1*x38 - x1*x44 + 16*x10*x4 + x11*x27*x30 + x12*x18 + x12*x6 + x12*x7 + x12 + 8*x13*x6 + 16*x13 - x15*x16 + x16*x17 + 8*x19*x8 + pow(x2, 2)*x45 + x2*x38 + x2*x44 - x20*x22 + x20*x23 + x20*x36 + x25*x30 + x26*x7 + x26 + x29*x7 + x29 + x31*x6 + x31 + x33*x34 + x33*x40 + x34*x35 + x35*x40 + x36*x41) ; // 71 +// const real_t x47 = x19 + x20 + x41 + x43 ; // 3 +// const real_t x48 = x46 + x47 ; // 1 +// const real_t x49 = -x15 - x17 ; // 2 +// const real_t x50 = x48 + x49 ; // 1 +// const real_t x51 = 1.0/x5 ; // 1 +// const real_t x52 = x51/(x18 + x39 + x42 + 4) ; // 4 +// const real_t x53 = x10 + x9 ; // 1 +// const real_t x54 = x15 + x17 + x53 ; // 2 +// const real_t x55 = x51/x6 ; // 1 +// const real_t x56 = (1.0/2.0)*x55 ; // 1 +// const real_t x57 = 1.0/(x7 + 2) ; // 2 +// const real_t x58 = x55*x57 ; // 1 +// const real_t x59 = x20*x7 ; // 1 +// const real_t x60 = x42*x8 ; // 1 +// const real_t x61 = 4*x14*x3 ; // 2 +// const real_t x62 = x41*x7 ; // 1 +// const real_t x63 = x30*x4 ; // 1 +// const real_t x64 = 2*x63 ; // 1 +// const real_t x65 = x6*x63 ; // 1 +// const real_t x66 = x3*x40*x5 ; // 2 +// const real_t x67 = x46*x7 ; // 1 +// const real_t x68 = x56*x57 ; // 1 +// const real_t x69 = -x46 + x47 ; // 1 +// const real_t Wp = x52*(x10 - x50 + x9) ; // 3 +// const real_t Sp = x56*(x48 + x54) ; // 2 +// const real_t Ip = -x58*(x50 + x53) ; // 3 +// //const real_t Rp = x68*(x22 + x23 - x36 - x59 - x60 + x61 - x62 - x64 - x65 + x66 - x67) ; // 11 +// const real_t Wm = x52*(-x19 - x20 - x41 - x43 + x46 + x54) ; // 6 +// const real_t Sm = x56*(x54 + x69) ; // 2 +// const real_t Im = -x58*(x49 + x53 + x69) ; // 4 +// //const real_t Rm = x68*(x22 + x23 - x36 - x59 - x60 + x61 - x62 - x64 - x65 + x66 + x67) ; // 11 +// // Opers = 255 + +// const real_t Rs = r0; +// const real_t Rp = Ip * x7 / 2. + Rs; +// const real_t Rm = Im * x7 / 2. + Rs; + +// const real_t eps = -1E-10; + +// printf("x4 x5 x6 x7 %e %e %e %e\n", x4, x5, x6, x7); +// printf("s %e %e %e %e \n", x0, x1, x2, x3); +// printf("p %e %e %e %e \n", Wp, Sp, Ip, Rp); +// printf("m %e %e %e %e \n", Wm, Sm, Im, Rm); + + + + + + +/** + CudaDeviceFunction void CalcPhi() + { + const real_t x0 = ; + const real_t x1 = (0,0); + const real_t x2 = (0,0); + const real_t x3 = 10.; + const real_t x4 = (0,0); + + phi[4] = x4; // N + // x0 = W^\star + // x1 = S^\star + // x2 = I^\star + // x3 = R^\star + + // x4 = N + //const real_t x4 = x1 + x2 + x3; + + // x5 = Beta + const real_t x5 = C_1; + + // x6 = Beta_w + const real_t x6 = C_2; + + // x7 = gamma + const real_t x7 = C_3; + + + const real_t x8 = x0*x5 ; // 1 + const real_t x9 = 2*x8 ; // 1 + const real_t x10 = x7*x8 ; // 1 + const real_t x11 = pow(x4, 2) ; // 1 + const real_t x12 = 16*x11 ; // 1 + const real_t x13 = x4*x8 ; // 1 + const real_t x14 = x5*x6 ; // 1 + const real_t x15 = x1*x14 ; // 1 + const real_t x16 = 8*x4 ; // 1 + const real_t x17 = x14*x2 ; // 1 + const real_t x18 = x6*x7 ; // 1 + const real_t x19 = x18*x4 ; // 1 + const real_t x20 = 4*x4 ; // 1 + const real_t x21 = x18*x5 ; // 1 + const real_t x22 = x1*x21 ; // 1 + const real_t x23 = x2*x21 ; // 1 + const real_t x24 = pow(x5, 2) ; // 1 + const real_t x25 = pow(x0, 2)*x24 ; // 2 + const real_t x26 = 4*x25 ; // 1 + const real_t x27 = pow(x6, 2) ; // 1 + const real_t x28 = 4*x11 ; // 1 + const real_t x29 = x27*x28 ; // 1 + const real_t x30 = pow(x7, 2) ; // 1 + const real_t x31 = x28*x30 ; // 1 + const real_t x32 = x0*x24 ; // 1 + const real_t x33 = x1*x32 ; // 1 + const real_t x34 = 4*x6 ; // 1 + const real_t x35 = x2*x32 ; // 1 + const real_t x36 = x30*x8 ; // 1 + const real_t x37 = x27*x5 ; // 1 + const real_t x38 = x20*x37 ; // 1 + const real_t x39 = 2*x6 ; // 1 + const real_t x40 = x39*x7 ; // 1 + const real_t x41 = x39*x4 ; // 1 + const real_t x42 = 2*x7 ; // 1 + const real_t x43 = x4*x42 ; // 1 + const real_t x44 = x37*x43 ; // 1 + const real_t x45 = x24*x27 ; // 1 + const real_t x46 = sqrt(pow(x1, 2)*x45 + 2*x1*x2*x45 - x1*x38 - x1*x44 + 16*x10*x4 + x11*x27*x30 + x12*x18 + x12*x6 + x12*x7 + x12 + 8*x13*x6 + 16*x13 - x15*x16 + x16*x17 + 8*x19*x8 + pow(x2, 2)*x45 + x2*x38 + x2*x44 - x20*x22 + x20*x23 + x20*x36 + x25*x30 + x26*x7 + x26 + x29*x7 + x29 + x31*x6 + x31 + x33*x34 + x33*x40 + x34*x35 + x35*x40 + x36*x41) ; // 71 + const real_t x47 = x19 + x20 + x41 + x43 ; // 3 + const real_t x48 = x46 + x47 ; // 1 + const real_t x49 = -x15 - x17 ; // 2 + const real_t x50 = x48 + x49 ; // 1 + const real_t x51 = 1.0/x5 ; // 1 + const real_t x52 = x51/(x18 + x39 + x42 + 4) ; // 4 + const real_t x53 = x10 + x9 ; // 1 + const real_t x54 = x15 + x17 + x53 ; // 2 + const real_t x55 = x51/x6 ; // 1 + const real_t x56 = (1.0/2.0)*x55 ; // 1 + const real_t x57 = 1.0/(x7 + 2) ; // 2 + const real_t x58 = x55*x57 ; // 1 + const real_t x59 = x20*x7 ; // 1 + const real_t x60 = x42*x8 ; // 1 + const real_t x61 = 4*x14*x3 ; // 2 + const real_t x62 = x41*x7 ; // 1 + const real_t x63 = x30*x4 ; // 1 + const real_t x64 = 2*x63 ; // 1 + const real_t x65 = x6*x63 ; // 1 + const real_t x66 = x3*x40*x5 ; // 2 + const real_t x67 = x46*x7 ; // 1 + const real_t x68 = x56*x57 ; // 1 + const real_t x69 = -x46 + x47 ; // 1 + const real_t Wp = x52*(x10 - x50 + x9) ; // 3 + const real_t Sp = x56*(x48 + x54) ; // 2 + const real_t Ip = -x58*(x50 + x53) ; // 3 + //const real_t Rp = x68*(x22 + x23 - x36 - x59 - x60 + x61 - x62 - x64 - x65 + x66 - x67) ; // 11 + const real_t Wm = x52*(-x19 - x20 - x41 - x43 + x46 + x54) ; // 6 + const real_t Sm = x56*(x54 + x69) ; // 2 + const real_t Im = -x58*(x49 + x53 + x69) ; // 4 + //const real_t Rm = x68*(x22 + x23 - x36 - x59 - x60 + x61 - x62 - x64 - x65 + x66 + x67) ; // 11 + // Opers = 255 + + const real_t Rs = (0,0); + const real_t Rp = Ip * x7 / 2. + Rs; + const real_t Rm = Im * x7 / 2. + Rs; + + const real_t eps = -1E-10; + if ( Sp >=eps && Ip >=eps && Rp >=eps && Wp >=eps) { + phi[1] = Sp < 0 ? 0 : Sp; + phi[2] = Ip < 0 ? 0 : Ip; + phi[3] = Rp < 0 ? 0 : Rp; + phi[0] = Wp < 0 ? 0 : Wp; + + } else if ( Sm >=eps&& Im >=eps && Rm >=eps && Wm >=eps ) { + phi[1] = Sm < 0 ? 0 : Sm; + phi[2] = Im < 0 ? 0 : Im; + phi[3] = Rm < 0 ? 0 : Rm; + phi[0] = Wm < 0 ? 0 : Wm; + + } else { + printf("x4 x5 x6 x7 %.18e %.18e %.18e %.18e\n", x4, x5, x6, x7); + printf("s %.18e %.18e %.18e %.18e %ld\n", x0, x1, x2, x3,sizeof(real_t)); + printf("p %e %e %e %e \n", Wp, Sp, Ip, Rp); + printf("m %e %e %e %e \n", Wm, Sm, Im, Rm); + + // Try Secant + + real_t X1[5], X2[5], f1[5], f2[5]; + + X1[0] = x0; + X1[1] = x1; + X1[2] = x2; + X1[3] = 1; + X1[4] = phi[4]; + + X2[0] = x0+0.001; + X2[1] = x1; + X2[2] = x2; + X2[3] = 1; + X2[4] = phi[4]; + + printf("Newton 0 %e %e %e %e \n", X2[0], X2[1], X2[2], X2[3]); + + for(int i =0; i < 50; i++) { + localCalcQ(X1, f1); + localCalcQ(X2, f2); + f1[0] = x0 - f1[0]/2. - X1[0]; + f1[1] = x1 - f1[1]/2. - X1[1]; + f1[2] = x2 - f1[2]/2. - X1[2]; + + f2[0] = x0 - f2[0]/2. - X2[0]; + f2[1] = x1 - f2[1]/2. - X2[1]; + f2[2] = x2 - f2[2]/2. - X2[2]; + + printf("Iter %e %e %e %e \n", f1[0], f1[1], f1[2], f1[3]); + + for(int k=0; k < 3; k++){ + const real_t tt = X1[k]; + X1[k] = X2[k] - f1[k] * (X2[k] - X1[k]) / (f2[k] - f1[k]); + X2[k] = tt; + } + + + + + printf("m %e %e %e %e \n", X2[0], X2[1], X2[2], X2[3]); + + } + + + printf("x4 x5 x6 x7 %e %e %e %e\n", x4, x5, x6, x7); + printf("m %e %e %e %e \n", X2[0], X2[1], X2[2], X2[3]); + + assert(0); + } + + + } + + CudaDeviceFunction void localCalcQ(real_t* localPhi, real_t* localQ) + { + + const real_t Beta = C_1; + const real_t Beta_w = C_2; + const real_t Gamma = C_3; + + + + const real_t W = localPhi[0]; + const real_t S = localPhi[1]; + const real_t I = localPhi[2]; + const real_t R = localPhi[3]; + + const real_t N = localPhi[4]; + + localQ[0] = Beta_w*(I-W); + localQ[1] = -Beta*S*W / N; + localQ[2] = (Beta*S*W/N - Gamma*I); + localQ[3] = Gamma*I; + localQ[4] = 0; //N + //printf("Phi %f %f %f %f\n", localPhi[0], localPhi[1], localPhi[2], localPhi[3]); + //printf("Q %f %f %f %f\n", q[0], q[1], q[2], q[3]); + } +**/ + + + +/** + // poor man`s secant method + real_t X1[3], X2[3], f1[3], f2[3]; + + + X1[0] = fabs(x0); + X1[1] = fabs(x1); + X1[2] = fabs(x2); + + localCalcQ(X1, f1, x4); + + X2[0] = fabs(X1[0] + f1[0]/2.); + X2[1] = fabs(X1[1] + f1[1]/2.); + X2[2] = fabs(X1[2] + f1[2]/2.); + + for(int i =0; i < 50; i++) { + localCalcQ(X1, f1, x4); + + f1[0] = x0 + f1[0]/2. - X1[0]; + f1[1] = x1 + f1[1]/2. - X1[1]; + f1[2] = x2 + f1[2]/2. - X1[2]; + + real_t epsilon_step = 0; + for(int k=0; k < 3; k++){ + if ( fabs(f1[k]) > epsilon_step ) { + epsilon_step = fabs(f1[k]); + } + } + if (epsilon_step < 1E-8){ + + for(int k=0; k < 3; k++){ + X1[k] = fabs(X1[k]); + } + phi[0] = X1[0]; + phi[1] = X1[1]; + phi[2] = X1[2]; + phi[3] = X1[2]*C_3/2. + x3; + + return; + } + + localCalcQ(X2, f2, x4); + + f2[0] = x0 + f2[0]/2. - X2[0]; + f2[1] = x1 + f2[1]/2. - X2[1]; + f2[2] = x2 + f2[2]/2. - X2[2]; + + for(int k=0; k < 3; k++){ + const real_t tt = X1[k]; + X1[k] = X1[k] - f1[k] * (X1[k] - X2[k]) / (f1[k] - f2[k]); + X2[k] = tt; + } + + } +*/ + \ No newline at end of file diff --git a/models/reaction/d2q9_reaction_diffusion_system/a.out b/models/reaction/d2q9_reaction_diffusion_system/a.out new file mode 100755 index 000000000..50ff6c304 Binary files /dev/null and b/models/reaction/d2q9_reaction_diffusion_system/a.out differ diff --git a/models/reaction/d2q9_reaction_diffusion_system/conf.mk b/models/reaction/d2q9_reaction_diffusion_system/conf.mk index 547da5d5d..e78e86d8e 100644 --- a/models/reaction/d2q9_reaction_diffusion_system/conf.mk +++ b/models/reaction/d2q9_reaction_diffusion_system/conf.mk @@ -1,3 +1,3 @@ ADJOINT=0 TEST=FALSE -OPT="(AllenCahn+SIR_ModifiedPeng+SIR_SimpleLaplace+SimpleDiffusion)-1" +OPT="(AllenCahn+SIR_ModifiedPeng+SIR_SimpleLaplace+SimpleDiffusion)*(Heun+Euler)-1" diff --git a/src/lib/boundary.R b/src/lib/boundary.R index f9a607f43..404d36066 100644 --- a/src/lib/boundary.R +++ b/src/lib/boundary.R @@ -140,7 +140,7 @@ ZouHeNew = function(EQ, f, direction, sign, order, group=f, known="rho",mom) { C(f[sel],fs[sel]) } -ZouHeRewrite = function(EQ, f, n, type=c("velocity","pressure","do nothing"), rhs) { +ZouHeRewrite = function(EQ, f, n, type=c("velocity","pressure","do nothing","first_moment"), rhs) { # --- Prepare arguments type=match.arg(type) @@ -183,6 +183,9 @@ ZouHeRewrite = function(EQ, f, n, type=c("velocity","pressure","do nothing"), rh # --- Set all velocity components eqn = V( fs %*% EQ2$U %*% diag(d)) rhs = rhs * sum(fs) + } else if (type == "first_moment") { + # --- Set all velocity components + eqn = V( fs %*% EQ2$U %*% diag(d)) } else stop("Unknown type in ZouHe") cat("/********* ", type, "-type Zou He boundary condition ****************/\n",sep=""); eqn = eqn - rhs; diff --git a/src/lib/feq.R b/src/lib/feq.R index 5f83f45f9..94a5317a5 100644 --- a/src/lib/feq.R +++ b/src/lib/feq.R @@ -108,3 +108,17 @@ WMRT_mat = function(U) { M = M %*% solve(R) M } + + +get.M.matrix = function(a) { + an = lapply(a@vec,names) + an = unique(do.call(c,an)) + an = setdiff(an,".M") + ret = lapply(a@vec,function(p) { nan = setdiff(an,names(p)); p[,nan] = 0; p}) + n = length(ret) + for (i in seq_len(n)) names(ret[[i]])[names(ret[[i]]) == ".M"] = paste0(".M",i) + m=ret[[1]] + for (i in seq_len(n-1)+1) m = merge(x=m,y=ret[[i]],all=TRUE,suffixes = c("a","b")) + m[is.na(m)] = 0 + as.matrix(m[,paste0(".M",1:n)]) +}