[Python] Robin boundary conditions for PDE - fipy

I'm trying to use fipy to solve a PDE system. In this case, I have a coupled of PDEs and I need to solve for two components of a vector field (m1 and m2). However, I dont need to solve these simultaneously. In the procedure, I solve first for m1 while using the previously computed values of m2. Also, these systems have Robin type boundary conditions. I was wondering if anybody had some experience with setting these type of boundary conditions (Robin) with fipy? I was taking a look to the documentation but the type of Robin BC shown there is different from mine (I dont have a convection term in the BC)

The expression for the PDE associated to m1 is given as PDE system.

The associated boundary conditions to this system are of Robin type and the expression is here boundary conditions

As mentioned, in this case I dont have to solve for m2 in this single system (that is done in the next step). That means, I use some pre-computed values for m2 while solving this specific system.

This is the way I have set my system using fipy:

n_mesh = self.local_mesh.faceNormals
betta = (1 - self.alpha) / (self.alpha)
alpha_f = FaceVariable(name="alpha_f", mesh=self.local_mesh, rank=1)
alpha_f.setValue(betta * n_mesh, where=self.bf_mask == True)
c11 = CellVariable(name="c11",mesh=self.local_mesh)
c11_f = FaceVariable(name="c11_f",mesh=self.local_mesh)

c1_d_c2 = CellVariable(name="c1_d_c2",mesh=self.local_mesh)
c1_d_c2_f = FaceVariable(name="c1_d_c2_f",mesh=self.local_mesh)

m2_previous = CellVariable(name="m2_previous",mesh=self.local_mesh)
m1 = CellVariable(name="m1",mesh=self.local_mesh)
C_T_P_1 = CellVariable(name="c_t_p_1",mesh=self.local_mesh,rank=1)
C_T_P_1_f = FaceVariable(name="c_t_p_1",mesh=self.local_mesh,rank=1)
m_1_b = ImplicitSourceTerm(coeff=alpha_f.divergence)
b1_f = FaceVariable(name="b1_f", mesh=self.local_mesh, rank=1)

c11.setValue(c11_val[0:self.N_cells])
c11_f.setValue(c11.arithmeticFaceValue)
c11_f.setValue(0.0,where=self.bf_mask)

c1_d_c2.setValue(c1_d_c2_val[0:self.N_cells])
c1_d_c2_f.setValue(c1_d_c2.arithmeticFaceValue)
c1_d_c2_f.setValue(0.0,where=self.bf_mask)

m2_previous.setValue(self.m2_iter[0:self.N_cells])
C_T_P_1.setValue(np.stack((r11[0:self.N_cells],r12[0:self.N_cells]),axis=0))
C_T_P_1_f.setValue(C_T_P_1.arithmeticFaceValue)
C_T_P_1_f.setValue(0.0,where=self.bf_mask)

b1_f_np = np.zeros(self.bf_mask.shape)
b1_f_np[self.bf_mask == True] = b_array[self.N_cells:, 0, 0]

r1_vec = np.zeros((2,self.bf_mask.shape[0]))
r1_vec[0,self.bf_mask==True] = r11[self.N_cells:]
r1_vec[1,self.bf_mask==True] = r12[self.N_cells:]
b1_f.setValue(betta * b1_f_np * n_mesh , where=self.bf_mask)

##We set the equations -m1
m1_eq = -DiffusionTerm(coeff=c11_f) == -(c1_d_c2_f * m2_previous.faceGrad).divergence + (C_T_P_1_f).divergence - b1_f.divergence + m_1_b


The first term DiffusionTerm is just the diffusion term associated to the LHS in the system I am trying to solve. The -(c1_d_c2_f * m2_previous.faceGrad).divergence corresponds to the same operation but with m2 (based on the documentation from fipy, this is how I should compute this soure term (https://www.ctcms.nist.gov/fipy/documentation/FAQ.html)

Then (C_T_P_1_f).divergence corresponds to the RHS in the original PDE expression.

Finally, b1_f.divergence and m_1_b corresponds to the corrections that I add to the boundary faces to account for the Robin boundary conditions.

The setup for m2 is exactly the same but of course by changing the role between m1 and m2.

Now the thing is that in principle I can solve the system via m1_eq.solve(var=m1). However, first at all the solutions that I obtain are not exactly what I would be expecting.

Also, this is part of an iterative algorithm and right after the first iteration, the solutions I get are completely off (really large values over the solution domain for both m1 and m2) which then tells me that there is something wrong with my system.

I have the feeling that the main problem is related to how I am setting the boundary contidions. For this, my idea was basically to take the original boundary condition, express the gradient of m1 as a function of the rest of terms in the boundary condition expression and then replace this as the flux associated to the boundary faces in the flux computation for the boundary cells.

This is the reason why once I assign the values via setValue I ue the mask self.bf_mask so that all terms associated to the standard expression are set to 0 for the boundary faces and all terms associated to the boundary conditions accordingly (in this case, just b1_f and m_1_b have non-zero entries corresponding to the boundary faces. Also, the expression that I am using for the flux at the boundaries from the boundary conditions is this: flux at boundaries

Also, in the examples given in fipy, the diffusion term is always given on the RHS, however in my case, I dont have any other non-source term. I was wondering if internally it flips the sign to the elements of the generated matrix? (I noticed that if I dont include the minus sign in front of the DifussionTerm then the solution that I obtain is flipped along x and y with respect to what I would be expecting)

Maybe anybody has some idea about what could be going on or what could I check for?

In advance, thanks for the comments and help.

Continue reading...