Your thoughts make sense and you're not missing an API. As I understand it, I think extracting the flux using the weights from the numerical scheme is the most accurate representation of the physical flux, as these account for the likely fact that the field does not vary linearly between cell centers. As you have found, we don't make it easy to get at this information. We'll have to think about a way to make it more accessible.
On Sep 28, 2015, at 8:59 AM, Urs Zimmermann <[email protected]> wrote: > Dear FiPy users and developers, > > first of all I would like to thank the developers for this great > software package. > > I am faced with a convection-diffusion equation > $$ > \frac{\partial \phi}{\partial t} = D \nabla^2 \phi + \nabla (\phi > \nabla V_{ext}) > $$ > with a strongly inhomogeneous concentration profile \phi in an external > potential V_ext. I would like to determine the flux. Mathematically, I > would have to calculate > $$ > j = -\nabla \phi - \phi \nabla V_{ext} > $$ > > However, a naive calculation of the flux in FiPy ("-phi.grad - V.grad * > phi" or "-phi.faceGrad - V.faceGrad * phi.faceValue") is strongly > deviating from what is expected, especially where the gradient in phi is > steep. See below for a minimal example. > > Considering the numerical finite volume scheme, I guess what I would > want to calculate is > $$ > (\vec{n} \cdot \vec{u})_f (\alpha_f \phi_P - (1-\alpha_f) \phi_A) + > \Gamma_f \frac{\phi_A - \phi_P}{d_{AP}}, > $$ > where \vec{u} is the gradient of V_{ext}. Summation of this field gives > the change of \phi in the cells. > Indeed, if I calculate the flux this way I get what I would expect. > However, to get the values for \alpha_f I extracted them from the > stencil in the ConvectionTerm. This is fine in the case of the simple > example given below but it becomes error-prone in a more complicated > context. > > So my questions are: > - Do these thoughts make sense? Is the physical flux well approximated > by the flux from the numerical scheme? > - Am I missing some API to get the flux directly? Is there a more > convenient way to get the flux as calculated from the numerical scheme? > > > Thanks in advance and best wishes, > > Urs > > > ##### Example ################### > # Equilibrium concentration profile > # in a harmonic potential > ############################## > > from fipy import * > > nx = 100 > dx = 1.0 > mesh = Grid1D(nx=nx, dx=dx) > x, = mesh.cellCenters() > D = 1.0 > L = nx*dx > > phi = CellVariable(name="solution variable", mesh=mesh, value=1.0, > hasOld=True) > V_ext = CellVariable(name="external potential", mesh=mesh, value=(x - > L/2.0)**2) > > eq = TransientTerm() == DiffusionTerm(coeff=D) + > ExponentialConvectionTerm(coeff=V_ext.faceGrad) > > # Solve for large timesteps to get to equilibrium > for step in range(10): > res = 10 > phi.updateOld() > while res > 1e-10: > res = eq.sweep(var=phi, dt=1000.0) > print(res) > print("next step") > > # Naive flux calculation > simple_flux1 = -phi.grad[0] - phi*(V_ext.grad[0]) > simple_flux2 = (-phi.faceGrad - phi.faceValue * V_ext.faceGrad)[0] > > # Calculate according to numeric CC-FVM scheme > # get stencil from ConvectionTerm to extract \alpha_f weight factor > stencil = eq.other.other.stencil['implicit'] > diffusive_flux = -phi.faceGrad[0] > conv_flux = V_ext.faceGrad[0][1:-1] * (stencil['cell 1 > diag'][1:-1]*phi.value[:-1] + stencil['cell 1 offdiag'][1:-1]*phi.value[1:]) > CCFVM_flux = diffusive_flux[1:-1] - conv_flux > > # Since this is the equilibrium configuration flux should be 0 > print("-------------------------------------------------") > print("Cell centered naive flux calculation: > max(|j|)={0}".format(abs(simple_flux1).max())) > print("Face centered naive flux calculation: > max(|j|)={0}".format(abs(simple_flux2).max())) > print("Face centered CC-FVM flux calculation: > max(|j|)={0}".format(abs(CCFVM_flux).max())) > > _______________________________________________ > fipy mailing list > [email protected] > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] _______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
