Skip to content

Commit

Permalink
Working on the solver. Neumann BCs are still not working properly. An…
Browse files Browse the repository at this point in the history 
…d, I need to re-implement the Dirichlet BCs so that I'm accounting for the boundary nodes correctly...
  • Loading branch information
@sap8b
sap8b committed Jul 22, 2021
1 parent 8672a09 commit df1cbb2
Show file tree
Hide file tree
Showing 7 changed files with 114 additions and 77 deletions.
6 changes: 6 additions & 0 deletions Diffusion2D_Library.sln
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,8 @@ VisualStudioVersion = 16.0.31402.337
MinimumVisualStudioVersion = 10.0.40219.1
Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "Diffusion2D_Library", "Diffusion2D_Library.csproj", "{0402B501-B293-47CA-AE2E-7939BD274CAE}"
EndProject
Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "Test_Diffusion2D_Library", "..\Test_Diffusion2D_Library\Test_Diffusion2D_Library.csproj", "{43A63030-4802-4B22-AD0F-2A42A17C6C09}"
EndProject
Global
GlobalSection(SolutionConfigurationPlatforms) = preSolution
Debug|Any CPU = Debug|Any CPU
Expand All @@ -15,6 +17,10 @@ Global
{0402B501-B293-47CA-AE2E-7939BD274CAE}.Debug|Any CPU.Build.0 = Debug|Any CPU
{0402B501-B293-47CA-AE2E-7939BD274CAE}.Release|Any CPU.ActiveCfg = Release|Any CPU
{0402B501-B293-47CA-AE2E-7939BD274CAE}.Release|Any CPU.Build.0 = Release|Any CPU
{43A63030-4802-4B22-AD0F-2A42A17C6C09}.Debug|Any CPU.ActiveCfg = Debug|Any CPU
{43A63030-4802-4B22-AD0F-2A42A17C6C09}.Debug|Any CPU.Build.0 = Debug|Any CPU
{43A63030-4802-4B22-AD0F-2A42A17C6C09}.Release|Any CPU.ActiveCfg = Release|Any CPU
{43A63030-4802-4B22-AD0F-2A42A17C6C09}.Release|Any CPU.Build.0 = Release|Any CPU
EndGlobalSection
GlobalSection(SolutionProperties) = preSolution
HideSolutionNode = FALSE
Expand Down
10 changes: 5 additions & 5 deletions DiffusionSimulators_1D.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -111,7 +111,7 @@ public void OneD_FE(BoundaryCondition_Del Lbc, BoundaryCondition_Del Rbc, Initia
RVector b = new(n);

// Define the A matrix
double nu = (D * dt) / (Math.Pow(dx, 2));
double nu = D * dt / Math.Pow(dx, 2);
double off_d_val = nu;
double diag_val = 1 - (2 * nu);
TridiagonalMatrix A = new(n, diag_val, off_d_val, off_d_val);
Expand Down Expand Up @@ -158,8 +158,8 @@ public void OneD_FE(BoundaryCondition_Del Lbc, BoundaryCondition_Del Rbc, Initia
break;
case BoundaryConditions.Neumann:
xnew = A.Dot(xold) + g(position, t * dt, xold);
xnew[0] = (-2 * nu) * Lbc(t * dt);
xnew[n - 1] = (2 * nu) * Rbc(t * dt);
xnew[0] = -2 * nu * Lbc(t * dt);
xnew[n - 1] = 2 * nu * Rbc(t * dt);
break;
case BoundaryConditions.Mixed:
xnew = A.Dot(xold) + b;
Expand All @@ -184,7 +184,7 @@ public void OneD_BE(BoundaryCondition_Del Lbc, BoundaryCondition_Del Rbc, Initia
RVector xnew = new(n);

// Define the A matrix
double nu = (D * dt) / (Math.Pow(dx, 2));
double nu = D * dt / Math.Pow(dx, 2);
double off_d_val = -nu;
double diag_val = 1 + (2 * nu);
TridiagonalMatrix A = new(n, diag_val, off_d_val, off_d_val);
Expand Down Expand Up @@ -254,7 +254,7 @@ public void OneD_CN(BoundaryCondition_Del Lbc, BoundaryCondition_Del Rbc, Initia
RVector bj1 = new(n);

// Setup the A matrix
double nu = (D * dt) / (Math.Pow(dx, 2));
double nu = D * dt / Math.Pow(dx, 2);
double off_d_val, diag_val;
off_d_val = -nu / 2;
diag_val = 1.0 + nu;
Expand Down
28 changes: 14 additions & 14 deletions DiffusionSimulators_2D_ConstantD.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -255,7 +255,7 @@ public void Solve(int n_time_steps, int output_interval)
RVector CT = new(ncols);
RVector CB = new(ncols);

double nu = (d_coeff * dt) / (2 * Math.Pow(dx, 2));
double nu = d_coeff * dt / (2 * Math.Pow(dx, 2));

// Define the A matrix for the explicit steps
double off_d_val = nu;
Expand Down Expand Up @@ -329,9 +329,9 @@ public void Solve(int n_time_steps, int output_interval)
cUpper2 = C_Im2[nrows - 2, ncols - 2];
}
//C3 = BCs[1].BoundaryFunction(t * dt, ncols);
CR = ((-2 * nu) * C1) + ((2 * nu) * C2); // + C3
CR[0] = ((-2 * nu) * cLower1) + ((2 * nu) * cLower2);
CR[ncols - 1] = ((-2 * nu) * cUpper1) + ((2 * nu) * cUpper2);
CR = (-2 * nu * C1) + (2 * nu * C2); // + C3
CR[0] = (-2 * nu * cLower1) + (2 * nu * cLower2);
CR[ncols - 1] = (-2 * nu * cUpper1) + (2 * nu * cUpper2);
break;
default:
break;
Expand Down Expand Up @@ -365,9 +365,9 @@ public void Solve(int n_time_steps, int output_interval)
cUpper2 = C_Im2[nrows - 2, 1];
}
//C3 = BCs[1].BoundaryFunction(t * dt, ncols);
CL = ((-2 * nu) * C1) + ((2 * nu) * C2); // + C3
CL[0] = ((-2 * nu) * cLower1) + ((2 * nu) * cLower2);
CL[ncols - 1] = ((-2 * nu) * cUpper1) + ((2 * nu) * cUpper2);
CL = (-2 * nu * C1) + (2 * nu * C2); // + C3
CL[0] = (-2 * nu * cLower1) + (2 * nu * cLower2);
CL[ncols - 1] = (-2 * nu * cUpper1) + (2 * nu * cUpper2);
break;
default:
break;
Expand Down Expand Up @@ -404,15 +404,15 @@ public void Solve(int n_time_steps, int output_interval)
}

//f12 = (tau / 2.0) * (fn + f0);
f12 = (tau / 2.0) * fn;
f12 = tau / 2.0 * fn;

// BCs
switch (BCs_Functions[0].TypeBC)
{
case ABoundaryCondition.dirichlet:
RVector gn = BCs_Functions[0].BoundaryFunction((t + 1) * dt, BCs_Functions[0].PositionVaries, BCs_Functions[0].PositionFixed);
RVector g0 = BCs_Functions[0].BoundaryFunction(t * dt, BCs_Functions[0].PositionVaries, BCs_Functions[0].PositionFixed);
CT = ((0.5 * B.Dot(gn)) + (0.5 * A.Dot(g0))); ////((0.5 * Bm.Dot(gn)) + (0.5 * Bp.Dot(g0))); //nu *
CT = (0.5 * B.Dot(gn)) + (0.5 * A.Dot(g0)); ////((0.5 * Bm.Dot(gn)) + (0.5 * Bp.Dot(g0))); //nu *
break;
case ABoundaryCondition.neumann:
RVector C1, C2; //,C3;
Expand All @@ -427,7 +427,7 @@ public void Solve(int n_time_steps, int output_interval)
C2 = C_Im2.GetRowVector(nrows - 2);
}
//C3 = BCs[1].BoundaryFunction(t * dt, ncols);
CT = ((-2 * nu) * C1) + ((2 * nu) * C2); // + C3
CT = (-2 * nu * C1) + (2 * nu * C2); // + C3
break;
default:
break;
Expand All @@ -437,7 +437,7 @@ public void Solve(int n_time_steps, int output_interval)
case ABoundaryCondition.dirichlet:
RVector gn = BCs_Functions[3].BoundaryFunction((t + 1) * dt, BCs_Functions[3].PositionVaries, BCs_Functions[3].PositionFixed);
RVector g0 = BCs_Functions[3].BoundaryFunction(t * dt, BCs_Functions[3].PositionVaries, BCs_Functions[3].PositionFixed);
CB = ((0.5 * B.Dot(gn)) + (0.5 * A.Dot(g0))); //((0.5 * Bm.Dot(gn)) + (0.5 * Bp.Dot(g0))); //nu *
CB = (0.5 * B.Dot(gn)) + (0.5 * A.Dot(g0)); //((0.5 * Bm.Dot(gn)) + (0.5 * Bp.Dot(g0))); //nu *
break;
case ABoundaryCondition.neumann:
RVector C1, C2; //,C3;
Expand All @@ -452,7 +452,7 @@ public void Solve(int n_time_steps, int output_interval)
C2 = C_Im2.GetRowVector(1);
}
//C3 = BCs[1].BoundaryFunction(t * dt, ncols);
CB = ((-2 * nu) * C1) + ((2 * nu) * C2); // + C3
CB = (-2 * nu * C1) + (2 * nu * C2); // + C3
break;
default:
break;
Expand Down Expand Up @@ -496,7 +496,7 @@ public void Solve(int n_time_steps, int output_interval)
C2 = C_Im2.GetColVector(ncols - 2);
}
//C3 = BCs[1].BoundaryFunction(t * dt, ncols);
CR = ((-2 * nu) * C1) + ((2 * nu) * C2); // + C3
CR = (-2 * nu * C1) + (2 * nu * C2); // + C3
break;
default:
break;
Expand Down Expand Up @@ -524,7 +524,7 @@ public void Solve(int n_time_steps, int output_interval)
C2 = C_Im2.GetColVector(1);
}
//C3 = BCs[1].BoundaryFunction(t * dt, ncols);
CL = ((-2 * nu) * C1) + ((2 * nu) * C2); // + C3
CL = (-2 * nu * C1) + (2 * nu * C2); // + C3
break;
default:
break;
Expand Down
Loading

0 comments on commit df1cbb2

Please sign in to comment.