Can I used the CellVariable method to map some lists onto the grid and put it into the coefficients of the equation?

[LiuDongDong1989]

I used the CellVariable method to map the list onto the grid and put it into the coefficients of the equation. I am not sure if this approach is correct.

    def UpdateEquation(self, t_tmp):
        print("Update PDE Equations...")
        coeff1=CellVariable(mesh=self.mesh, value=self.h)
        coeff2=CellVariable(mesh=self.mesh, value=self.h * self.u)
        coeff3=CellVariable(mesh=self.mesh, value=self.h * self.DH)
        Tc=self.eval_Tc()
        Dc=self.eval_Dc(t_tmp)
        Sss_tmp=self.eval_Sss(Tc, Dc)
        Sss=CellVariable(mesh=self.mesh, value=Sss_tmp) 
        R = self.P * numerix.cos(self.slopeAngle) - self.i
        self.eqc = ( TransientTerm(var=self.c, coeff=coeff1) + 
                    ConvectionTerm(var=self.c, coeff=coeff2 * [[1]]) ) == (
                        DiffusionTerm (var=self.c, coeff=coeff3) + 
                        Sss -
                        self.c * R
                    )

In the above code, h, u, Sss_tmp are Python lists and have the same length as 1D Mesh. The code runs perfectly, but I don’t know if the equation is represented correctly, and are there any better methods to do this if I use 2D mesh?

[guyer]

Assuming that $h$, $u$, and maybe $D_H$, are field variables, you certainly should not write

        coeff1=CellVariable(mesh=self.mesh, value=self.h)
        coeff2=CellVariable(mesh=self.mesh, value=self.h * self.u)
        coeff3=CellVariable(mesh=self.mesh, value=self.h * self.DH)

This will create three static variables which do not change as self.h, self.u, and self.DH change.

Instead, write

        coeff1=self.h
        coeff2=self.h * self.u
        coeff3=self.h * self.DH

or just

        self.eqc = ( TransientTerm(var=self.c, coeff=self.h) + 
                    ConvectionTerm(var=self.c, coeff=self.h * self.u * [[1]]) ) == (
                        DiffusionTerm (var=self.c, coeff=self.h * self.DH) + 
                        Sss -
                        self.c * R
                    )

As to Sss, it all depends on what self.eval_Tc(), self.eval_Dc(), and self.eval_Sss() look like. There are probably better ways to write those, that involve FiPy CellVariable objects and not Python lists, but I have no way of knowing.

[LiuDongDong1989]

h, u, Sss_tmp are just Python lists, not the field variables , and their values were changed at each time step.
I used theself.UpdateEquation(self, t_tmp) to update the equation.

[guyer]

Okaaaaay, but if they have a value corresponding to every cell on mesh, then they are field variables. I'm suggesting that it might be possible to write self.eval_Tc(), self.eval_Dc(), and self.eval_Sss() such that that was all taken care of automatically. That would be how FiPy is intended to be used.

Creating a new equation at every time step is not a good design. If you absolutely can't just write variable expressions, then you would be better to declare, outside the scope of UpdateEquation():

        self.Tc = CellVariable(mesh=self.mesh)
        self.Dc = Variable(0.)
        self.Sss = CellVariable(mesh=self.mesh)
        R = self.P * numerix.cos(self.slopeAngle) - self.i
        self.eqc = ( TransientTerm(var=self.c, coeff=self.h) + 
                    ConvectionTerm(var=self.c, coeff=self.h * self.u * [[1]]) ) == (
                        DiffusionTerm (var=self.c, coeff=self.h * self.DH) + 
                        Sss -
                        self.c * R
                    )

and then

    def UpdateEquation(self, t_tmp):
        print("Update PDE coefficients...")
        self.Tc.value = self.eval_Tc()
        self.Dc.value = self.eval_Dc(t_tmp)
        self.Sss = self.eval_Sss(self.Tc, self.Dc)

The equation will update on its own.

There are probably still-better ways of doing this, but you'd need to tell me what self.eval_Tc(), self.eval_Dc(), and self.eval_Sss() actually do.

[LiuDongDong1989]

    
    def eval_Dc(self, t_tmp):
        C = self.eval_CumulativeRainfall(t_tmp)    
        if self.Kτ < self.Kτ_peak and not self.exceeded:
            self.Kτ = self.Kτ_initial * math.exp(self.λincrease * (C - self.Cp) * 100)
        else:
            if self.Cv is None:
                self.Cv = C
                self.exceeded = True
            self.Kτ = self.Kτ_steady + (self.Kτ_peak - self.Kτ_steady) * math.exp(-self.λdecrease * (C - self.Cv) * 100)
        
        Dc = self.Kτ * self.τ ** 1.5
        return Dc
    
    def eval_Tc(self): 
        Tc = [
            161 * numerix.tan(self.slopeAngle) * (q_i*60)**0.5 /60
            if (q_i * 60) <= 0.046 
            else 16320 * numerix.tan(self.slopeAngle) * (q_i*60)**2 /60
            for q_i in self.q
        ]
        # print('Tc=', Tc)
        return Tc

    def eval_Sss(self, Tc, Dc):
        Sss = [
            -self.U * (self.c.value[i] * self.q[i] - Tc_i) 
            if Tc_i < self.c.value[i] * self.q[i] 
            else Dc_i * (1 - self.c.value[i] * self.q[i] / Tc_i) 
            for i, (Tc_i, Dc_i) in enumerate(zip(Tc, Dc))
        ]
        # print('Sss=', Sss)
        return Sss

[guyer]

def eval_Kτ:
        C = self.eval_CumulativeRainfall(t_tmp)    
        if self.Kτ < self.Kτ_peak and not self.exceeded:
            Kτ = self.Kτ_initial * math.exp(self.λincrease * (C - self.Cp) * 100)
        else:
            if self.Cv is None:
                self.Cv = C
                self.exceeded = True
            Kτ = self.Kτ_steady + (self.Kτ_peak - self.Kτ_steady) * math.exp(-self.λdecrease * (C - self.Cv) * 100)
        return Kτ

def UpdateEquation(self, t_tmp):
    print("Update PDE coefficients...")
    self.Kτ.value = self.eval_Kτ(t_tmp)

self.Kτ = CellVariable(mesh=mesh)
self.Dc = self.Kτ * self.τ ** 1.5
self.Tc = ((161 * (self.q*60)**0.5) * ((self.q * 60) <= 0.046)
           + (16320  * (self.q*60)**2) * ((self.q * 60) > 0.046)) * numerix.tan(self.slopeAngle) / 60
self.Sss = ((-self.U * (self.c * self.q - Tc)) * (self.Tc < self.c * self.q)
            + (self.Dc * (1 - self.c * self.q / self.Tc) ) * (self.Tc >= self.c * self.q))

eval_Kτ() can probably be written as a variable expression, too, but I'd need to know more about the things in it.