Solution not updating in loop

General Discussion
1

So I want to have a simulation where I have a small circular source where the vector field oscillates up and down to see how the rest of the field responds over time.

There is a second order time derivative, I am using backwards finite difference formula for this.
In the code Eold is the prior iteration and Eoldold is two prior iterations.

The problem is my code is not updating with time, I have set the field to oscillate with time but it isn’t, any ideas why are much appreciated!

Edit: I think it might be something to do with the boundary conditions, if I change those I get some change but it is still completely broken.

Edit2: I think it might definitely be the boundary conditions, how to change them consistently between iterations is the problem. At the moment no vector field is generated from “Source” despite it being in the on() statements.

Edit3: definitely the boundary conditions! If I vary them with space I get a vector field, but varying them with t doesn’t work! any idea’s what I can do about this?

Edit 4: So the problem is with how to update Eold and Eoldold in the for loop, I can put the following lines in the for loop:

real[int] bcEold = varfEold(0,Xh);
real[int] bcEoldold = varfEoldold(0,Xh);

and it will update the values at the boundary but does not propagate the system, I.E. It’s like resetting the initial conditions and taking a single iteration each time, which is not what I am after.
So how do I propagate the whole system?

real T = 100000;
real dT = 0.06;
real t = 0;

real w = 2*pi*20;
// Set initial conditions
real ExInitialGradient = 1;
real EyInitialGradient = 1; 
real ExInitial = 0.5;
real EyInitial = 0;

// For mesh setup
int nn = 30;
int boxheight = 2;
int boxwidth = 4;
int source = 10;
int box = 1;
//constants
real mu0 = 1.25663706e-6;
real epsilon0 = 8.85e-12;



real a = 0.1; // radius of the circle;
// Set up mesh
border Source(t = 0,2*pi){x = a*cos(t)+0.3;y=a*sin(t)+0.3;label = source;}
border rightwall(t = 0,boxheight){x = boxwidth; y = t;label = box;}
border leftwall(t = 0,boxheight){x = 0;y = boxheight - t;label = box; }
border bottom(t = 0,boxwidth){x = t;y =0;label = box;}
border top(t = 0,boxwidth){x = boxwidth - t;y = boxheight;label = box;}
mesh Th = buildmesh(Source(10)+rightwall(10)+leftwall(10)+top(10)+bottom(10)); // Our mesh

fespace Vh(Th,P1); // For the scalar x and y components of the field
fespace Xh(Th,[P1,P1]); // For the x and y components, it is what we will actually solve for
// Declare everything
Vh Ex, Ey, Exold, Eyold, Exoldold, Eyoldold, fx, fy, intensity;

// Set initial conditions
Exoldold = ExInitial;
Eyoldold = EyInitial;
// interpolate for Eold;
Exold = Exoldold + dT*ExInitialGradient;
Eyold = Eyoldold + dT*EyInitialGradient;

func epsilonr = 10*exp(-10*y); // relative permitivity as a function of height;

//Next we set up are variable forms for A[Ex,Ey] = B[Exold,Eyold] + C[Exoldold,Eyoldold] + bcEold + bcEoldold;

// varf for A
varf aa([Ex,Ey],[fx,fy]) = int2d(Th)(
	(dy(fx)*dx(Ey)) + 
	(dx(fy)*dy(Ex)) +
	(dx(fx)*dx(Ex)) +
	(dy(fy)*dy(Ey))
) - int2d(Th)(
	mu0*epsilon0*epsilonr*(1/dT^2)*(
		Ex*fx + Ey*fy
	)
) + on(source,Ex = 0, Ey = sin(w*t));

// varf for B
varf bb([Exold,Eyold],[fx,fy]) = int2d(Th)(
	-2*(1/dT^2)*mu0*epsilon0*epsilonr*(
		Exold*fx + Eyold*fy
	)
);
//varf for C
varf cc([Exoldold,Eyoldold],[fx,fy]) = int2d(Th)(
	(1/dT^2)*mu0*epsilon0*epsilonr*(
		Exoldold*fx + Eyoldold*fy
	)
);

//Set up Bc's
varf varfEold([Exold,Eyold],[fx,fy]) = on(source, Exold = 0, Eyold = sin(w*(t-dT)));
varf varfEoldold([Exoldold,Eyoldold],[fx,fy]) = on(source, Exoldold = 0, Eyoldold = sin(w*(t-2*dT)));

//Set up matrices and vectors;
matrix A = aa(Xh,Xh,solver = UMFPACK);
matrix B = bb(Xh,Xh);
matrix C = cc(Xh,Xh);
//Set up vectors
real[int] sol(Xh.ndof), solOld(Xh.ndof),solOldOld(Xh.ndof); // This will be for E, Eold, and Eoldold respectively;
real[int] aux(Xh.ndof); // This will be our temporary store
//Boundary vectors
real[int] bcEold = varfEold(0,Xh);
real[int] bcEoldold = varfEoldold(0,Xh);


//Set up for the loop, we need some temporary variables
Xh [w1,w2] = [Exold,Eyold];
Xh [w3,w4] = [Exoldold,Eyoldold];
// Time loop!
solOld = w1[]; // For E(t-dt)
solOldOld = w3[]; // For E(t-2*dT)
int m, M = T/dT;
for (m = 0; m < M; m++){
	aux = B*solOld; aux += C*solOldOld; aux += bcEold; aux += bcEoldold;
	sol = A^-1 * aux; 
	w1[] = sol;
	Ex = w1; Ey = w2;
	intensity = Ex^2 + Ey^2;
	t += dT;
	Exold = Ex;
	Eyold = Ey;
	Exoldold = Exold;
	Eyoldold = Eyold;
	solOld = sol; solOldOld = solOld;
	cout << "\n\nNext loop, t =  " << t <<"\n";
	plot(intensity, [Ex,Ey],wait = true, fill = true);

	
}
2

the value Eyoldold is wrong !!!

1 Like