Skip to content

Commit

Permalink
Revised the A, B_row, and B_col matrices yet again...
Browse files Browse the repository at this point in the history
  • Loading branch information
@sap8b
sap8b committed Jun 20, 2021
1 parent 519f0e5 commit ac47b2c
Show file tree
Hide file tree
Showing 4 changed files with 114 additions and 96 deletions.
2 changes: 1 addition & 1 deletion Diffusion2D_Library.csproj
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@
<TargetFramework>net5.0</TargetFramework>
<PackageLicenseExpression>MIT OR Apache-2.0</PackageLicenseExpression>
<GeneratePackageOnBuild>true</GeneratePackageOnBuild>
<Version>0.1.8</Version>
<Version>0.1.9</Version>
<Description>Solves the parabolic partial differential equation for diffusion or heat ransfer in 1 or 2 dimensions.</Description>
<PackageIcon>2D_Region.png</PackageIcon>
<PackageIconUrl />
Expand Down
8 changes: 7 additions & 1 deletion Diffusion2D_Library.sln
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,9 @@ Microsoft Visual Studio Solution File, Format Version 12.00
# Visual Studio Version 16
VisualStudioVersion = 16.0.31402.337
MinimumVisualStudioVersion = 10.0.40219.1
Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "Diffusion2D_Library", "Diffusion2D_Library.csproj", "{0402B501-B293-47CA-AE2E-7939BD274CAE}"
Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "Diffusion2D_Library", "Diffusion2D_Library.csproj", "{0402B501-B293-47CA-AE2E-7939BD274CAE}"
EndProject
Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "TestProject_Diffusion", "..\TestProject_Diffusion\TestProject_Diffusion.csproj", "{11FBB911-0555-4649-85A8-A5F4CAA34DB5}"
EndProject
Global
GlobalSection(SolutionConfigurationPlatforms) = preSolution
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
{11FBB911-0555-4649-85A8-A5F4CAA34DB5}.Debug|Any CPU.ActiveCfg = Debug|Any CPU
{11FBB911-0555-4649-85A8-A5F4CAA34DB5}.Debug|Any CPU.Build.0 = Debug|Any CPU
{11FBB911-0555-4649-85A8-A5F4CAA34DB5}.Release|Any CPU.ActiveCfg = Release|Any CPU
{11FBB911-0555-4649-85A8-A5F4CAA34DB5}.Release|Any CPU.Build.0 = Release|Any CPU
EndGlobalSection
GlobalSection(SolutionProperties) = preSolution
HideSolutionNode = FALSE
Expand Down
192 changes: 98 additions & 94 deletions DiffusionSimulators_2D_MatrixD.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -185,9 +185,16 @@ public DiffusionSimulators_2D_MatrixD(RMatrix D, double dx, double dy, int nx, i
double nu0 = dt / (2 * Math.Pow(dx, 2));

RVector nu = new(nx);
RVector off_d_val_l = new(nx - 1);
RVector off_d_val_u = new(nx - 1);
RVector main = new(nx);
RVector off_d_val_l0 = new(nx - 1);
RVector off_d_val_u0 = new(nx - 1);
RVector main0 = new(nx);
RVector off_d_val_l1 = new(nx - 1);
RVector off_d_val_u1 = new(nx - 1);
RVector main1 = new(nx);
RVector off_d_val_l2 = new(nx - 1);
RVector off_d_val_u2 = new(nx - 1);
RVector main2 = new(nx);

double D_minus, D_plus;
for (int i = 0; i < nx; i++)
{
Expand All @@ -198,102 +205,97 @@ public DiffusionSimulators_2D_MatrixD(RMatrix D, double dx, double dy, int nx, i
//================================================
// A Matrix first
//================================================
D_plus = 0.5 * (D_row[1] - D_row[0]);
main[0] = main[0] = 1.0 - ((D_plus * nu0));
D_plus = 0.5 * (D_row[1] + D_row[0]);
main0[0] = 1.0 - ((D_plus * nu0));

D_minus = 0.5 * (D_row[nx - 2] - D_row[0]);
main[nx - 1] = 1.0 - ((D_minus * nu0));
D_minus = 0.5 * (D_row[nx - 2] + D_row[0]);
main0[nx - 1] = 1.0 - ((D_minus * nu0));

for (int j = 1; j < nx - 2; j++)
for (int j = 1; j < nx - 1; j++)
{
D_minus = 0.5 * (D_row[j - 1] - D_row[j]);
D_plus = 0.5 * (D_row[j + 1] - D_row[j]);
main[j] = 1.0 - ((D_minus * nu0) + (D_plus * nu0));
D_minus = 0.5 * (D_row[j - 1] + D_row[j]);
D_plus = 0.5 * (D_row[j + 1] + D_row[j]);
main0[j] = 1.0 - ((D_minus * nu0) + (D_plus * nu0));
}

//off_d_val_l[0] = 0.5 * (D_row[0] - D_row[1]) * nu0;
//off_d_val_l[nx - 2] = nu[nx - 1];
for (int j = 1; j < nx; j++)
{
D_minus = 0.5 * (D_row[j - 1] - D_row[j]);
off_d_val_l[j - 1] = nu0 * D_minus;
D_minus = 0.5 * (D_row[j] + D_row[j - 1]);
off_d_val_l0[j - 1] = nu0 * D_minus;
}

//off_d_val_u[0] = nu[0]; off_d_val_u[nx - 2] = nu[nx - 1];
for (int j = 0; j < nx - 1; j++)
{
D_plus = 0.5 * (D_row[j + 1] - D_row[j]);
off_d_val_u[j] = nu0 * D_plus;
D_plus = 0.5 * (D_row[j + 1] + D_row[j]);
off_d_val_u0[j] = nu0 * D_plus;
}

A[i] = new TridiagonalMatrix(main, off_d_val_l, off_d_val_u);
A[i] = new TridiagonalMatrix(main0, off_d_val_l0, off_d_val_u0);
////================================================
//// The B-row Matrix next
////================================================
//main = 1 + (2 * nu);

D_plus = 0.5 * (D_row[1] - D_row[0]);
main[0] = main[0] = 1.0 + ((D_plus * nu0));
D_plus = 0.5 * (D_row[1] + D_row[0]);
main1[0] = 1.0 + ((D_plus * nu0));

D_minus = 0.5 * (D_row[nx - 2] - D_row[0]);
main[nx - 1] = 1.0 + ((D_minus * nu0));
D_minus = 0.5 * (D_row[nx - 2] + D_row[0]);
main1[nx - 1] = 1.0 + ((D_minus * nu0));

for (int j = 1; j < nx - 2; j++)
for (int j = 1; j < nx - 1; j++)
{
D_minus = 0.5 * (D_row[j - 1] - D_row[j]);
D_plus = 0.5 * (D_row[j + 1] - D_row[j]);
main[j] = 1.0 + ((D_minus * nu0) + (D_plus * nu0));
D_minus = 0.5 * (D_row[j - 1] + D_row[j]);
D_plus = 0.5 * (D_row[j + 1] + D_row[j]);
main1[j] = 1.0 + ((D_minus * nu0) + (D_plus * nu0));
}

for (int j = 1; j < nx; j++)
{
D_minus = 0.5 * (D_row[j - 1] - D_row[j]);
off_d_val_l[j - 1] = -nu0 * D_minus;
D_minus = 0.5 * (D_row[j - 1] + D_row[j]);
off_d_val_l1[j - 1] = -nu0 * D_minus;
}

//off_d_val_u[0] = nu[0]; off_d_val_u[nx - 2] = nu[nx - 1];
for (int j = 0; j < nx - 1; j++)
{
D_plus = 0.5 * (D_row[j + 1] - D_row[j]);
off_d_val_u[j] = -nu0 * D_plus;
D_plus = 0.5 * (D_row[j + 1] + D_row[j]);
off_d_val_u1[j] = -nu0 * D_plus;
}

B_row[i] = new TridiagonalMatrix(main, off_d_val_l, off_d_val_u);
B_row[i] = new TridiagonalMatrix(main1, off_d_val_l1, off_d_val_u1);
////================================================
//// The B-column Matrix last
////================================================
RVector D_col = D.GetColVector(i);

D_plus = 0.5 * (D_col[1] - D_col[0]);
main[0] = main[0] = 1.0 + ((D_plus * nu0));
D_plus = 0.5 * (D_col[1] + D_col[0]);
main2[0] = 1.0 + ((D_plus * nu0));

D_minus = 0.5 * (D_col[nx - 2] - D_col[0]);
main[nx - 1] = 1.0 + ((D_minus * nu0));
D_minus = 0.5 * (D_col[nx - 2] + D_col[0]);
main2[nx - 1] = 1.0 + ((D_minus * nu0));

for (int j = 1; j < nx - 2; j++)
for (int j = 1; j < nx - 1; j++)
{
D_minus = 0.5 * (D_col[j - 1] - D_col[j]);
D_plus = 0.5 * (D_col[j + 1] - D_col[j]);
main[j] = 1.0 + ((D_minus * nu0) + (D_plus * nu0));
D_minus = 0.5 * (D_col[j - 1] + D_col[j]);
D_plus = 0.5 * (D_col[j + 1] + D_col[j]);
main2[j] = 1.0 + ((D_minus * nu0) + (D_plus * nu0));
}

for (int j = 1; j < nx; j++)
{
D_minus = 0.5 * (D_col[j - 1] - D_col[j]);
off_d_val_l[j - 1] = -nu0 * D_minus;
D_minus = 0.5 * (D_col[j - 1] + D_col[j]);
off_d_val_l2[j - 1] = -nu0 * D_minus;
}

//off_d_val_u[0] = nu[0]; off_d_val_u[nx - 2] = nu[nx - 1];
for (int j = 0; j < nx - 1; j++)
{
D_plus = 0.5 * (D_col[j + 1] - D_col[j]);
off_d_val_u[j] = -nu0 * D_plus;
D_plus = 0.5 * (D_col[j + 1] + D_col[j]);
off_d_val_u2[j] = -nu0 * D_plus;
}

////for (int j = 0; j < nx - 1; j++) { off_d_val_l[j] = nu[j]; }
////for (int j = 1; j < nx; j++) { off_d_val_u[j - 1] = nu[j]; }

B_col[i] = new TridiagonalMatrix(main, -off_d_val_l, -off_d_val_u);
B_col[i] = new TridiagonalMatrix(main2, off_d_val_l2, off_d_val_u2);
}

C_Initial = I0(CF_2D.xposition_matrix, CF_2D.yposition_matrix);
Expand Down Expand Up @@ -399,9 +401,16 @@ public DiffusionSimulators_2D_MatrixD(RMatrix D, double dx, double dy, int nx, i
double nu0 = dt / (2 * Math.Pow(dx, 2));

RVector nu = new(nx);
RVector off_d_val_l = new(nx - 1);
RVector off_d_val_u = new(nx - 1);
RVector main = new(nx);
RVector off_d_val_l0 = new(nx - 1);
RVector off_d_val_u0 = new(nx - 1);
RVector main0 = new(nx);
RVector off_d_val_l1 = new(nx - 1);
RVector off_d_val_u1 = new(nx - 1);
RVector main1 = new(nx);
RVector off_d_val_l2 = new(nx - 1);
RVector off_d_val_u2 = new(nx - 1);
RVector main2 = new(nx);

double D_minus, D_plus;
for (int i = 0; i < nx; i++)
{
Expand All @@ -412,102 +421,97 @@ public DiffusionSimulators_2D_MatrixD(RMatrix D, double dx, double dy, int nx, i
//================================================
// A Matrix first
//================================================
D_plus = 0.5 * (D_row[1] - D_row[0]);
main[0] = main[0] = 1.0 - ((D_plus * nu0));
D_plus = 0.5 * (D_row[1] + D_row[0]);
main0[0] = 1.0 - ((D_plus * nu0));

D_minus = 0.5 * (D_row[nx - 2] - D_row[0]);
main[nx - 1] = 1.0 - ((D_minus * nu0));
D_minus = 0.5 * (D_row[nx - 2] + D_row[0]);
main0[nx - 1] = 1.0 - ((D_minus * nu0));

for (int j = 1; j < nx - 2; j++)
for (int j = 1; j < nx - 1; j++)
{
D_minus = 0.5 * (D_row[j - 1] - D_row[j]);
D_plus = 0.5 * (D_row[j + 1] - D_row[j]);
main[j] = 1.0 - ((D_minus * nu0) + (D_plus * nu0));
D_minus = 0.5 * (D_row[j - 1] + D_row[j]);
D_plus = 0.5 * (D_row[j + 1] + D_row[j]);
main0[j] = 1.0 - ((D_minus * nu0) + (D_plus * nu0));
}

//off_d_val_l[0] = 0.5 * (D_row[0] - D_row[1]) * nu0;
//off_d_val_l[nx - 2] = nu[nx - 1];
for (int j = 1; j < nx; j++)
{
D_minus = 0.5 * (D_row[j - 1] - D_row[j]);
off_d_val_l[j - 1] = nu0 * D_minus;
D_minus = 0.5 * ( D_row[j] + D_row[j - 1]);
off_d_val_l0[j - 1] = nu0 * D_minus;
}

//off_d_val_u[0] = nu[0]; off_d_val_u[nx - 2] = nu[nx - 1];
for (int j = 0; j < nx - 1; j++)
{
D_plus = 0.5 * (D_row[j + 1] - D_row[j]);
off_d_val_u[j] = nu0 * D_plus;
D_plus = 0.5 * (D_row[j + 1] + D_row[j]);
off_d_val_u0[j] = nu0 * D_plus;
}

A[i] = new TridiagonalMatrix(main, off_d_val_l, off_d_val_u);
A[i] = new TridiagonalMatrix(main0, off_d_val_l0, off_d_val_u0);
////================================================
//// The B-row Matrix next
////================================================
//main = 1 + (2 * nu);

D_plus = 0.5 * (D_row[1] - D_row[0]);
main[0] = main[0] = 1.0 + ((D_plus * nu0));
D_plus = 0.5 * (D_row[1] + D_row[0]);
main1[0] = 1.0 + ((D_plus * nu0));

D_minus = 0.5 * (D_row[nx - 2] - D_row[0]);
main[nx - 1] = 1.0 + ((D_minus * nu0));
D_minus = 0.5 * (D_row[nx - 2] + D_row[0]);
main1[nx - 1] = 1.0 + ((D_minus * nu0));

for (int j = 1; j < nx - 2; j++)
for (int j = 1; j < nx - 1; j++)
{
D_minus = 0.5 * (D_row[j - 1] - D_row[j]);
D_plus = 0.5 * (D_row[j + 1] - D_row[j]);
main[j] = 1.0 + ((D_minus * nu0) + (D_plus * nu0));
D_minus = 0.5 * (D_row[j - 1] + D_row[j]);
D_plus = 0.5 * (D_row[j + 1] + D_row[j]);
main1[j] = 1.0 + ((D_minus * nu0) + (D_plus * nu0));
}

for (int j = 1; j < nx; j++)
{
D_minus = 0.5 * (D_row[j - 1] - D_row[j]);
off_d_val_l[j - 1] = -nu0 * D_minus;
D_minus = 0.5 * (D_row[j - 1] + D_row[j]);
off_d_val_l1[j - 1] = -nu0 * D_minus;
}

//off_d_val_u[0] = nu[0]; off_d_val_u[nx - 2] = nu[nx - 1];
for (int j = 0; j < nx - 1; j++)
{
D_plus = 0.5 * (D_row[j + 1] - D_row[j]);
off_d_val_u[j] = -nu0 * D_plus;
D_plus = 0.5 * (D_row[j + 1] + D_row[j]);
off_d_val_u1[j] = -nu0 * D_plus;
}

B_row[i] = new TridiagonalMatrix(main, off_d_val_l, off_d_val_u);
B_row[i] = new TridiagonalMatrix(main1, off_d_val_l1, off_d_val_u1);
////================================================
//// The B-column Matrix last
////================================================
RVector D_col = D.GetColVector(i);

D_plus = 0.5 * (D_col[1] - D_col[0]);
main[0] = main[0] = 1.0 + ((D_plus * nu0));
D_plus = 0.5 * (D_col[1] + D_col[0]);
main2[0] = 1.0 + ((D_plus * nu0));

D_minus = 0.5 * (D_col[nx - 2] - D_col[0]);
main[nx - 1] = 1.0 + ((D_minus * nu0));
D_minus = 0.5 * (D_col[nx - 2] + D_col[0]);
main2[nx - 1] = 1.0 + ((D_minus * nu0));

for (int j = 1; j < nx - 2; j++)
for (int j = 1; j < nx - 1; j++)
{
D_minus = 0.5 * (D_col[j - 1] - D_col[j]);
D_plus = 0.5 * (D_col[j + 1] - D_col[j]);
main[j] = 1.0 + ((D_minus * nu0) + (D_plus * nu0));
D_minus = 0.5 * (D_col[j - 1] + D_col[j]);
D_plus = 0.5 * (D_col[j + 1] + D_col[j]);
main2[j] = 1.0 + ((D_minus * nu0) + (D_plus * nu0));
}

for (int j = 1; j < nx; j++)
{
D_minus = 0.5 * (D_col[j - 1] - D_col[j]);
off_d_val_l[j - 1] = -nu0 * D_minus;
D_minus = 0.5 * (D_col[j - 1] + D_col[j]);
off_d_val_l2[j - 1] = -nu0 * D_minus;
}

//off_d_val_u[0] = nu[0]; off_d_val_u[nx - 2] = nu[nx - 1];
for (int j = 0; j < nx - 1; j++)
{
D_plus = 0.5 * (D_col[j + 1] - D_col[j]);
off_d_val_u[j] = -nu0 * D_plus;
D_plus = 0.5 * (D_col[j + 1] + D_col[j]);
off_d_val_u2[j] = -nu0 * D_plus;
}

////for (int j = 0; j < nx - 1; j++) { off_d_val_l[j] = nu[j]; }
////for (int j = 1; j < nx; j++) { off_d_val_u[j - 1] = nu[j]; }

B_col[i] = new TridiagonalMatrix(main, -off_d_val_l, -off_d_val_u);
B_col[i] = new TridiagonalMatrix(main2, off_d_val_l2, off_d_val_u2);
}

C_Initial = I0(CF_2D.xposition_matrix, CF_2D.yposition_matrix);
Expand Down
Loading

0 comments on commit ac47b2c

Please sign in to comment.