Trouble compiling weak form terms in 3d NSE simulation

General Discussion
1

Hello, I am trying to implement a numerical solver for the 3d Eithier-Steinman flow, and I am getting a curious error when I try to compile and run the code. The following is my .edp file

//Load necessary packages for our solver
load “UMFPACK64”
load “msh3”

//Define some globally used variables
real a = 1.25;
real d = 1.0;
real chi = 1;
real nu = 1.0e-4;
real delta = 0.1;
real t = 0.0;
real dt = 0.01;
int Gn = -1; //(-1 for NSE; 0,1, and 2 for Deconvolution order 0,1, and 2 respectively)
int nSteps = 100;
int NLerror;
int h = 10;

//Define initial condition / true solution functions and their derivatives
func real u1True(real t){
-a * (exp(a * x) * sin(a * y + d * z) + exp(a * z) * sin(a * x + d * y))*exp(-nu * d^2 * t);
}
func real u2True(real t){
-a * (exp(a * y) * sin(a * z + d * x) + exp(a * x) * sin(a * y + d * z))*exp(-nu * d^2 * t);
}
func real u3True(real t){
-a * (exp(a * z) * sin(a * x + d * y) + exp(a * y) * sin(a * z + d * x))*exp(-nu * d^2 * t);
}
func real pTrue(real t){
-a^2 / 2 * (exp(2 * a * x) + exp(2 * a * y) + exp(2 * a * z) + 2 * sin(a * x + d * y) * cos(a * z + d * x) * exp(a * (y + z)) + 2 * sin(a * y + d * z) * cos(a * x + d * y) * exp(a * (z + x))
+ 2 * sin(a * z + d * x) * cos(a * y + d * z) * exp(a * (x + y))) * exp(-nu * d^2 * t);
}
func real dxu1True(real t){
-a * (a * exp(a * x) * sin(a * y + d * z) - a * exp(a * z) * sin(a * x + d * y)) * exp(-nu * d^2 * t);
}
func real dyu1True(real t){
-a * (a * exp(a * x) * cos(a * y + d * z) - d * exp(a * z) * sin(a * x + d * y)) * exp(-nu * d^2 * t);
}
func real dzu1True(real t){
-a * (d * exp(a * x) * cos(a * y + d * z) + a * exp(a * z) * cos(a * x + d * y)) * exp(-nu * d^2 * t);
}
func real dxu2True(real t){
-a * (d * exp(a * y) * cos(a * z + d * x) + a * exp(a * x) * cos(a * y + d * z)) * exp(-nu * d^2 * t);
}
func real dyu2True(real t){
-a * (a * exp(a * y) * sin(a * z + d * x) - a * exp(a * x) * sin(a * y + d * z)) * exp(-nu * d^2 * t);
}
func real dzu2True(real t){
-a * (a * exp(a * y) * cos(a * z + d * x) - d * exp(a * x) * sin(a * x + d * y)) * exp(-nu * d^2 * t);
}
func real dxu3True(real t){
-a * (a * exp(a * z) * cos(a * x + d * y) - d * exp(a * y) * sin(a * z + d * z)) * exp(-nu * d^2 * t);
}
func real dyu3True(real t){
-a * (d * exp(a * z) * cos(a * x + d * y) + a * exp(a * y) * cos(a * z + d * z)) * exp(-nu * d^2 * t);
}
func real dzu3True(real t){
-a * (a * exp(a * z) * sin(a * x + d * y) - a * exp(a * y) * sin(a * z + d * y)) * exp(-nu * d^2 * t);
}

//Define the mesh and our finite element spaces
mesh3 Th = cube(h, h, h);
fespace XQh(Th, [P2, P2, P2, P1]);

//Define some convenient macros
macro Grad(u1, u2, u3) [dx(u1), dy(u1), dz(u1), dx(u2), dy(u2), dz(u2), dx(u3), dy(u3), dz(u3)] //

macro div(u1, u2, u3) (dx(u1) + dy(u2) + dz(u3)) //

macro b(u1, u2, u3, v1, v2, v3, w1, w2, w3) ((u1 * dx(v1) + u2 * dy(v1) + u3 * dz(v1)) * w1 + (u1 * dx(v2) + u2 * dy(v2) + u3 * dz(v2)) * w2 + (u1 * dx(v3) + u2 * dy(v3) + u3 * dz(v3)) * w3) //

macro bstar(u1, u2, u3, v1, v2, v3, w1, w2, w3) (0.5 * (b(u1, u2, u3, v1, v2, v3, w1, w2, w3) - b(u1, u2, u3, w1, w2, w3, v1, v2, v3)) //

macro c(u1, u2, u3, v1, v2, v3, w1, w2, w3) ((dx(u1) * v1 + 0.5 * (dx(u2) + dy(u1)) * v2 + 0.5 * (dx(u3) + dz(u1)) * v3) * w1 + (0.5 * (dx(u2) + dy(u1)) * v1 + dy(u2) * v2 + 0.5 * (dy(u3) + dz(u2)) * v3) * w2 + (0.5 * (dx(u3) + dz(u1)) * v1 + 0.5 * (dy(u3) + dz(u2)) * v2 + dz(u3) * v3) * w3 + div(u1,u2,u3) * (v1 * w1 + v2 * w2 + v3 * w3)) //

//Define some finite element space variables for our problems
XQh [u1, u2, u3, p],
[v1, v2, v3, q],
[u1NLold, u2NLold, u3NLold, pNLold],
[u1TLold, u2TLold, u3TLold, pTLold],
[e1, e2, e3, pe],
[ub1, ub2, ub3, pb],
[ubTest1, ubTest2, ubTest3, pbTest],
[ubb1, ubb2, ubb3, pbb],
[ubbTest1, ubbTest2, ubbTest3, pbbTest],
[ubbb1, ubbb2, ubbb3, pbbb],
[ubbbTest1, ubbbTest2, ubbbTest3, pbbbTest],
[uD1, uD2, uD3, pD],
[uDold1, uDold2, uDold3, pDold];

//Stokes Filter for u bar
//
problem StokesFP1([ub1, ub2, ub3, pb], [ubTest1, ubTest2, ubTest3, pbTest], solver = UMFPACK) =
int3d(Th)(
delta^2 * (Grad(ub1, ub2, ub3) '* Grad(ubTest1, ubTest2, ubTest3))
+ (ub1 * ubTest1 + ub2 * ubTest2 + ub3 * ubTest3)
- div(ubTest1, ubTest2, ubTest3) * pb
- div(ub1, ub2, ub3) * pbTest
+ pb * pbTest * (1e-10)
)
- int3d(Th)(
(u1 * ubTest1 + u2ubTest2 + u3ubTest3)
)
//Boundary Condition
+ on(1, 2, 3, 4, 5, 6, ub1 = u1True(t), ub2 = u2True(t), ub3 = u3True(t));
;

//Stokes Filter for u double bar
//
problem StokesFP2([ubb1, ubb2, ubb3, pbb], [ubbTest1, ubbTest2, ubbTest3, pbbTest], solver = UMFPACK) =
int3d(Th)(
delta^2 * (Grad(ubb1, ubb2, ubb3) '* Grad(ubbTest1, ubbTest2, ubbTest3))
+ (ubb1 * ubbTest1 + ubb2 * ubbTest2 + ubb3 * ubbTest3)
- div(ubbTest1, ubbTest2, ubbTest3) * pbb
- div(ubb1, ubb2, ubb3) * pbbTest
+ pbb * pbbTest * (1e-10)
)
- int3d(Th)(
(ub1 * ubbTest1 + ub2ubbTest2 + ub3ubbTest3)
)
//Boundary Condition
+ on(1, 2, 3, 4, 5, 6, ubb1 = u1True(t), ubb2 = u2True(t), ubb3 = u3True(t));
;

//Stokes Filter for u triple bar
//
problem StokesFP3([ubbb1, ubbb2, ubbb3, pbbb], [ubbbTest1, ubbbTest2, ubbbTest3, pbbbTest], solver = UMFPACK) =
int3d(Th)(
delta^2 * (Grad(ubbb1, ubbb2, ubbb3) '* Grad(ubbbTest1, ubbbTest2, ubbbTest3))
+ (ubbb1 * ubbbTest1 + ubbb2 * ubbbTest2 + ubbb3 * ubbbTest3)
- div(ubbbTest1, ubbbTest2, ubbbTest3) * pbbb
- div(ubbb1, ubbb2, ubbb3) * pbbbTest
+ pbbb * pbbbTest * (1e-10)
)
- int3d(Th)(
(ubb1 * ubbbTest1 + ubb2ubbbTest2 + ubb3ubbbTest3)
)
//Boundary Condition
+ on(1, 2, 3, 4, 5, 6, ubbb1 = u1True(t), ubbb2 = u2True(t), ubbb3 = u3True(t));
;

//Time Relaxation Model for SKEW
//
problem TRMSKEW([u1, u2, u3, p], [v1, v2, v3, q], solver = UMFPACK) =
int3d(Th)(
1 / dt * (u1 * v1 + u2 * v2 + u3 * v3)
+ 0.5 * bstar(u1NLold, u2NLold, u3NLold, u1, u2, u3, v1, v2, v3)
+ 0.5 * bstar(u1, u2, u3, u1TLold, u2TLold, u3TLold, v1, v2, v3)
+ 0.5 * bstar(u1TLold, u2TLold, u3TLold, u1, u2, u3, v1, v2, v3)
+ 0.5 * nu * (Grad(u1, u2, u3) '* Grad(v1, v2, v3))
- 0.5 * p * div(v1, v2, v3)
+ 0.5 * chi * (u1 * v1 + u2 * v2 + u3 * v3)
+ 0.5 * q * div(u1, u2, u3)
+ p * q * (1e-10)
)
- int3d(Th)(
1 / dt * (u1TLold * v1 + u2TLold * v2 + u3TLold * v3)
- 0.5 * bstar(u1TLold, u2TLold, u3TLold, u1, u2, u3, v1, v2, v3)
- 0.5 * nu * (Grad(u1TLold, u2TLold, u3TLold) '* Grad(v1, v2, v3))
- 0.5 * pTLold * div(v1, v2, v3)
- 0.5 * chi * (u1TLold * v1 + u2TLold * v2 + u3TLold * v3)
+ 0.5 * chi * ((uD1 + uDold1) * v1 + (uD2 + uDold2) * v2 + (uD3 + uDold3) * v3)
- 0.5 * q * div(u1TLold, u2TLold, u3TLold)
)
//Boundary Condition
+ on(1, 2, 3, 4, 5, 6, u1 = u1True(t), u2 = u2True(t), u3 = u3True(t));
;

//Initialize the time loop
[u1, u2, u3, p] = [u1True(0.0), u2True(0.0), u3True(0.0), pTrue(0.0)];

//TIME LOOP////////////////////////////////////////////////////////////////////////////////////

for(iter = 1; iter <= nSteps; iter++){

//Update the time, and set the previous time step variables to the old solution
t = t + dt;
[u1TLold, u2TLold, u3TLold, pTLold] = [u1, u2, u3, p];

//Count the time loop iterations
NLerror = 1;
cout << "Time loop iteration: " << iter << "------------------------" << endl;

//FIXED POINT ITERATION LOOP///////////////////////////////////////////////////////////////
for(itnl = 1; itnl <=MaxNlIt; itnl++){

	//Count the non-linear loop iterations, and set the previous non linear step variables to the old solution
	cout << "Non-Linear loop iteration: " << itnl << endl;
	[u1NLold, u2NLold, u3NLold, pNLold] = [u1, u2, u3, p];
	[uDold1, uDold2, uDold3, pDold] = [uD1, uD2, uD3, pD];
	
	//Solve the problem based on what order of deconvolution we want
	//NSE
	if(Gn == -1){
		chi =0.0;
		[uD1, uD2, uD3, pD] = [u1, u2, u3, p];
		TRMSKEW;
	}
	
	//TRM, N=0
	if(Gn == 0){
		StokesFP1;
		[uD1, uD2, uD3, pD] = [ub1, ub2, ub3, pb];
		TRMSKEW;
	}
	
	//TRM, N=1
	if(Gn == 1){
		StokesFP1;
		StokesFP2;
		[uD1, uD2, uD3, pD] = [2 * ub1 - ubb1, 2 * ub2 - ubb2, 2 * ub3 - ubb3, pb];
		TRMSKEW;
	}
	
	//TRM, N=2
	if(Gn == 2){
		StokesFP1;
		StokesFP2;
		StokesFP3;
		[uD1, uD2, uD3, pD] = [3 * ub1 - 3 * ubb1 + ubbb1, 3 * ub2 - 3 * ubb2 + ubbb2, 3 * ub3 - 3 * ubb3 + ubbb3, pb];
		TRMSKEW;
	}
	
	//Set the error values and output them to the user
	[e1, e2, e3, pe] = [abs(u1 - u1NLold), abs(u2 - u2NLold), abs(u3 - u3NLold), abs(p - pNLold)];
	cout << "e1[].max " << e1[].max << "     e2[].max " << e2[].max << endl;
	cout << "e1[].min " << e1[].min << "     e2[].min " << e2[].min << endl;
	cout << " " << endl;
	
	//Break the loop if the difference between consecutive iterations is small enough, otherwise reloop
	if ((e1[].max < 1.e-6) && (e2[].max < 1.e-6)) {
		break;
	}
	NLerror = NLerror + 1;
}

//Break out of the time loop if the non-linear iteration loop took too long to converge
cout << "Non-Linear loop iteration: " << NLerror << endl;
if(Nlerror >= MaxNlIt){
	cout << "Non-Linear iteration failed for time step iteration = " << iter << endl;
	break;
}

//Plot various iterations as the loop runs
if(iter == 20 || iter == 40 || iter == 60 || iter == 80 || iter == 100){
	plot(Th, [u1, u2, u3], ps = "ES_plot_" + iter + ".eps");
}

}

The error seems to be in the TRMSTEP problem, for some reason it will not get passed the gradient term in the second integral, which is weird to me since in the above integral the compiler takes in a nearly identical term with different inputs and compiles it no problem. Moreover, the StokesFP problems all have similar terms as well and they get compiled with no problems. Is there some sort of variable type mismatching going on here that I am not seeing? Any help provided is appreciated, thank you in advance.

2

Please surround the whole code between three single quotes ``` so that it’s easier for us to copy/paste.

Like this!

Thanks.

3

My apologies, here is the code again. This time it has been updated, I found a bug that solved my previous question, but now I am unsure why my compiler is giving me an error for the boundary condition on the TRMSTEP problem, especially considering the boundary condition for the other 3 problems is exactly the same and it compiles for those. Again any help is appreciated

//Load necessary packages for our solver
load "UMFPACK64"
load "msh3"

//Define some globally used variables
real a = 1.25;
real d = 1.0;
real chi = 1;
real nu = 1.0e-4;
real delta = 0.1;
real t = 0.0;
real dt = 0.01;
int Gn = -1;		//(-1 for NSE; 0,1, and 2 for Deconvolution order 0,1, and 2 respectively)
int nSteps = 100;
int NLerror;
int h = 10;

//Define initial condition / true solution functions and their derivatives
func real u1True(real t){
	-a * (exp(a * x) * sin(a * y + d * z) + exp(a * z) * sin(a * x + d * y))*exp(-nu * d^2 * t);
}
func real u2True(real t){
	-a * (exp(a * y) * sin(a * z + d * x) + exp(a * x) * sin(a * y + d * z))*exp(-nu * d^2 * t);
}
func real u3True(real t){
	-a * (exp(a * z) * sin(a * x + d * y) + exp(a * y) * sin(a * z + d * x))*exp(-nu * d^2 * t);
}
func real pTrue(real t){
	-a^2 / 2 * (exp(2 * a * x) + exp(2 * a * y) + exp(2 * a * z) + 2 * sin(a * x + d * y) * cos(a * z + d * x) * exp(a * (y + z)) + 2 * sin(a * y + d * z) * cos(a * x + d * y) * exp(a * (z + x)) 
	+ 2 * sin(a * z + d * x) * cos(a * y + d * z) * exp(a * (x + y))) * exp(-nu * d^2 * t);
}
func real dxu1True(real t){
	-a * (a * exp(a * x) * sin(a * y + d * z) - a * exp(a * z) * sin(a * x + d * y)) * exp(-nu * d^2 * t);
}
func real dyu1True(real t){
	-a * (a * exp(a * x) * cos(a * y + d * z) - d * exp(a * z) * sin(a * x + d * y)) * exp(-nu * d^2 * t);
}
func real dzu1True(real t){
	-a * (d * exp(a * x) * cos(a * y + d * z) + a * exp(a * z) * cos(a * x + d * y)) * exp(-nu * d^2 * t);
}
func real dxu2True(real t){
	-a * (d * exp(a * y) * cos(a * z + d * x) + a * exp(a * x) * cos(a * y + d * z)) * exp(-nu * d^2 * t);
}
func real dyu2True(real t){
	-a * (a * exp(a * y) * sin(a * z + d * x) - a * exp(a * x) * sin(a * y + d * z)) * exp(-nu * d^2 * t);
}
func real dzu2True(real t){
	-a * (a * exp(a * y) * cos(a * z + d * x) - d * exp(a * x) * sin(a * x + d * y)) * exp(-nu * d^2 * t);
}
func real dxu3True(real t){
	-a * (a * exp(a * z) * cos(a * x + d * y) - d * exp(a * y) * sin(a * z + d * z)) * exp(-nu * d^2 * t);
}
func real dyu3True(real t){
	-a * (d * exp(a * z) * cos(a * x + d * y) + a * exp(a * y) * cos(a * z + d * z)) * exp(-nu * d^2 * t);
}
func real dzu3True(real t){
	-a * (a * exp(a * z) * sin(a * x + d * y) - a * exp(a * y) * sin(a * z + d * y)) * exp(-nu * d^2 * t);
}

//Define the mesh and our finite element spaces
mesh3 Th = cube(h, h, h);
fespace XQh(Th, [P2, P2, P2, P1]);

//Define some convenient macros
macro Grad(u1, u2, u3) [dx(u1), dy(u1), dz(u1), dx(u2), dy(u2), dz(u2), dx(u3), dy(u3), dz(u3)] //

macro div(u1, u2, u3) (dx(u1) + dy(u2) + dz(u3)) //

macro b(u1, u2, u3, v1, v2, v3, w1, w2, w3) ((u1 * dx(v1) + u2 * dy(v1) + u3 * dz(v1)) * w1 + (u1 * dx(v2) + u2 * dy(v2) + u3 * dz(v2)) * w2 + (u1 * dx(v3) + u2 * dy(v3) + u3 * dz(v3)) * w3) //

macro bstar(u1, u2, u3, v1, v2, v3, w1, w2, w3) (0.5 * (b(u1, u2, u3, v1, v2, v3, w1, w2, w3) - b(u1, u2, u3, w1, w2, w3, v1, v2, v3)) //

macro c(u1, u2, u3, v1, v2, v3, w1, w2, w3) ((dx(u1) * v1 + 0.5 * (dx(u2) + dy(u1)) * v2 + 0.5 * (dx(u3) + dz(u1)) * v3) * w1 + (0.5 * (dx(u2) + dy(u1)) * v1 + dy(u2) * v2 + 0.5 * (dy(u3) + dz(u2)) * v3) * w2 + (0.5 * (dx(u3) + dz(u1)) * v1 + 0.5 * (dy(u3) + dz(u2)) * v2 + dz(u3) * v3) * w3 + div(u1,u2,u3) * (v1 * w1 + v2 * w2 + v3 * w3)) //

//Define some finite element space variables for our problems
XQh [u1, u2, u3, p],
	[v1, v2, v3, q],
	[u1NLold, u2NLold, u3NLold, pNLold],
	[u1TLold, u2TLold, u3TLold, pTLold],
	[e1, e2, e3, pe],
	[ub1, ub2, ub3, pb],
	[ubTest1, ubTest2, ubTest3, pbTest],
	[ubb1, ubb2, ubb3, pbb],
	[ubbTest1, ubbTest2, ubbTest3, pbbTest],
	[ubbb1, ubbb2, ubbb3, pbbb],
	[ubbbTest1, ubbbTest2, ubbbTest3, pbbbTest],
	[uD1, uD2, uD3, pD],
	[uDold1, uDold2, uDold3, pDold];
	
//Stokes Filter for u bar
//
problem StokesFP1([ub1, ub2, ub3, pb], [ubTest1, ubTest2, ubTest3, pbTest], solver = UMFPACK) = 
	int3d(Th)(
		delta^2 * (Grad(ub1, ub2, ub3) '* Grad(ubTest1, ubTest2, ubTest3))
		+ (ub1 * ubTest1 + ub2 * ubTest2 + ub3 * ubTest3)
		- div(ubTest1, ubTest2, ubTest3) * pb 
		- div(ub1, ub2, ub3) * pbTest
		+ pb * pbTest * (1.0e-10)
	)
	- int3d(Th)(
		(u1 * ubTest1 + u2*ubTest2 + u3*ubTest3)
	)
	//Boundary Condition
	+ on(1, 2, 3, 4, 5, 6, ub1 = u1True(t), ub2 = u2True(t), ub3 = u3True(t))
;

//Stokes Filter for u double bar
//
problem StokesFP2([ubb1, ubb2, ubb3, pbb], [ubbTest1, ubbTest2, ubbTest3, pbbTest], solver = UMFPACK) = 
	int3d(Th)(
		delta^2 * (Grad(ubb1, ubb2, ubb3) '* Grad(ubbTest1, ubbTest2, ubbTest3))
		+ (ubb1 * ubbTest1 + ubb2 * ubbTest2 + ubb3 * ubbTest3)
		- div(ubbTest1, ubbTest2, ubbTest3) * pbb 
		- div(ubb1, ubb2, ubb3) * pbbTest
		+ pbb * pbbTest * (1.0e-10)
	)
	- int3d(Th)(
		(ub1 * ubbTest1 + ub2*ubbTest2 + ub3*ubbTest3)
	)
	//Boundary Condition
	+ on(1, 2, 3, 4, 5, 6, ubb1 = u1True(t), ubb2 = u2True(t), ubb3 = u3True(t))
;

//Stokes Filter for u triple bar
//
problem StokesFP3([ubbb1, ubbb2, ubbb3, pbbb], [ubbbTest1, ubbbTest2, ubbbTest3, pbbbTest], solver = UMFPACK) = 
	int3d(Th)(
		delta^2 * (Grad(ubbb1, ubbb2, ubbb3) '* Grad(ubbbTest1, ubbbTest2, ubbbTest3))
		+ (ubbb1 * ubbbTest1 + ubbb2 * ubbbTest2 + ubbb3 * ubbbTest3)
		- div(ubbbTest1, ubbbTest2, ubbbTest3) * pbbb 
		- div(ubbb1, ubbb2, ubbb3) * pbbbTest
		+ pbbb * pbbbTest * (1.0e-10)
	)
	- int3d(Th)(
		(ubb1 * ubbbTest1 + ubb2*ubbbTest2 + ubb3*ubbbTest3)
	)
	//Boundary Condition
	+ on(1, 2, 3, 4, 5, 6, ubbb1 = u1True(t), ubbb2 = u2True(t), ubbb3 = u3True(t))
;

//Time Relaxation Model for SKEW
//
problem TRMSKEW([u1, u2, u3, p], [v1, v2, v3, q], solver = UMFPACK) =
	int3d(Th)(
		1 / dt * (u1 * v1 + u2 * v2 + u3 * v3)
		+ 0.5 * bstar(u1NLold, u2NLold, u3NLold, u1, u2, u3, v1, v2, v3)
		+ 0.5 * bstar(u1, u2, u3, u1TLold, u2TLold, u3TLold, v1, v2, v3)
		+ 0.5 * bstar(u1TLold, u2TLold, u3TLold, u1, u2, u3, v1, v2, v3)
		+ 0.5 * nu * (Grad(u1, u2, u3) '* Grad(v1, v2, v3))
		- 0.5 * p * div(v1, v2, v3)
		+ 0.5 * chi * (u1 * v1 + u2 * v2 + u3 * v3)
		+ 0.5 * q * div(u1, u2, u3)
		+ p * q * (1.0e-10)
	)
	- int3d(Th)(
		1 / dt * (u1TLold * v1 + u2TLold * v2 + u3TLold * v3)
		- 0.5 * bstar(u1TLold, u2TLold, u3TLold, u1TLold, u2TLold, u3TLold, v1, v2, v3)
		- 0.5 * nu * (Grad(u1TLold, u2TLold, u3TLold) '* Grad(v1, v2, v3))
		- 0.5 * pTLold * div(v1, v2, v3)
		- 0.5 * chi * (u1TLold * v1 + u2TLold * v2 + u3TLold * v3)
		+ 0.5 * chi * ((uD1 + uDold1) * v1 + (uD2 + uDold2) * v2 + (uD3 + uDold3) * v3)
		- 0.5 * q * div(u1TLold, u2TLold, u3TLold)
	)
	//Boundary Condition
	+ on(1, 2, 3, 4, 5, 6, u1 = u1True(t), u2 = u2True(t), u3 = u3True(t))
;

//Initialize the time loop
[u1, u2, u3, p] = [u1True(0.0), u2True(0.0), u3True(0.0), pTrue(0.0)];

//TIME LOOP////////////////////////////////////////////////////////////////////////////////////

for(iter = 1; iter <= nSteps; iter++){

	//Update the time, and set the previous time step variables to the old solution
	t = t + dt;
	[u1TLold, u2TLold, u3TLold, pTLold] = [u1, u2, u3, p];
	
	//Count the time loop iterations
	NLerror = 1;
	cout << "Time loop iteration: " << iter << "------------------------" << endl;
	
	//FIXED POINT ITERATION LOOP///////////////////////////////////////////////////////////////
	for(itnl = 1; itnl <=MaxNlIt; itnl++){
	
		//Count the non-linear loop iterations, and set the previous non linear step variables to the old solution
		cout << "Non-Linear loop iteration: " << itnl << endl;
		[u1NLold, u2NLold, u3NLold, pNLold] = [u1, u2, u3, p];
		[uDold1, uDold2, uDold3, pDold] = [uD1, uD2, uD3, pD];
		
		//Solve the problem based on what order of deconvolution we want
		//NSE
		if(Gn == -1){
			chi =0.0;
			[uD1, uD2, uD3, pD] = [u1, u2, u3, p];
			TRMSKEW;
		}
		
		//TRM, N=0
		if(Gn == 0){
			StokesFP1;
			[uD1, uD2, uD3, pD] = [ub1, ub2, ub3, pb];
			TRMSKEW;
		}
		
		//TRM, N=1
		if(Gn == 1){
			StokesFP1;
			StokesFP2;
			[uD1, uD2, uD3, pD] = [2 * ub1 - ubb1, 2 * ub2 - ubb2, 2 * ub3 - ubb3, pbb];
			TRMSKEW;
		}
		
		//TRM, N=2
		if(Gn == 2){
			StokesFP1;
			StokesFP2;
			StokesFP3;
			[uD1, uD2, uD3, pD] = [3 * ub1 - 3 * ubb1 + ubbb1, 3 * ub2 - 3 * ubb2 + ubbb2, 3 * ub3 - 3 * ubb3 + ubbb3, pbbb];
			TRMSKEW;
		}
		
		//Set the error values and output them to the user
		[e1, e2, e3, pe] = [abs(u1 - u1NLold), abs(u2 - u2NLold), abs(u3 - u3NLold), abs(p - pNLold)];
		cout << "e1[].max " << e1[].max << "     e2[].max " << e2[].max << endl;
		cout << "e1[].min " << e1[].min << "     e2[].min " << e2[].min << endl;
		cout << " " << endl;
		
		//Break the loop if the difference between consecutive iterations is small enough, otherwise reloop
		if ((e1[].max < 1.e-6) && (e2[].max < 1.e-6)) {
			break;
		}
		NLerror = NLerror + 1;
	}
	
	//Break out of the time loop if the non-linear iteration loop took too long to converge
	cout << "Non-Linear loop iteration: " << NLerror << endl;
	if(Nlerror >= MaxNlIt){
		cout << "Non-Linear iteration failed for time step iteration = " << iter << endl;
		break;
	}
	
	//Plot various iterations as the loop runs
	if(iter == 20 || iter == 40 || iter == 60 || iter == 80 || iter == 100){
		plot(Th, [u1, u2, u3], ps = "ES_plot_" + iter + ".eps");
	}
}
4

There is a problem with your macro bstar. If you comment all your calls to this macro, the code runs fine.

1 Like
5

Yep that was it. One missing parenthesis symbol, I found it and replaced it and now my code runs! Now to see if it gives me good results… Thank you for your help!

1 Like