Question regarding power level for milonga

36 views
Skip to first unread message

Ferraro, Diego Ernesto (INR)

unread,
Sep 3, 2018, 9:11:06 AM9/3/18
to was...@seamplex.com

Hi guys,

 I’ve a 2 short questions related with units in milonga that I suppose Is quite straightforward to you.

 

1 - I’ve a 3-D problem and I would like to set up the total power of the system to scale neutron fluxes. I understand that I should set the power and eFiss XS as:

MATERIAL xxx {

eSigmaF_1=yyy eSigmaF_2=zzz }

and

power = ppp

 

The problem is that I cannot figure out if the if eSigmaF is really the key I should use and which units I should use

à I suppose that the units are something like (per group)  Flux [n/cm2s] = Power[W/cm3] / SigmaF [1/cm] / Efiss [J/fiss] , thus Power should be in Watts and Efiss should be in W per fission (per group). Thus I’m using

eSigmaF_X = SigmaF_X [1/cm] * efiss [J/fiss]  (with power in Watts).

Anyway, is there any input card reference where I could get the available options for MATERIAL, power, etc (I cannot find it in the repository)?

 

2-I have some doubts with the scattering matrix. I’m using the P0 condensed data for the scattering matrix, but I suppose that the self-scattering (i.e. diagonal) should be transport corrected. So I just changed to P1 for the diagonal.  Nevertheless, when I change the S1.1 and S2.2 values I see no effect in keff. I’m missing something?

 

Thanks,

Diego

 

 

César Pablo Camusso

unread,
Sep 3, 2018, 9:43:42 AM9/3/18
to was...@seamplex.com
Diego:
1)Yes eSigmaF is the macroscopic fission cross section multiplied by the energy per fission.

power is the power in wats.
2)SigmaS_1.1 and SigmaS_2.2 are the diagonal elements.
Do you use SigmaA_g or SigmaT_g for the absortion cross sections?
It can change how the cross section are read.
from diffusion_elements.c:

        // absorcion
        xi = gsl_matrix_get(A, g, g);
        if (material_xs->SigmaT[g] != NULL && material_xs->SigmaT[g]->n_tokens != 0) {
          xi += wasora_evaluate_expression(material_xs->SigmaT[g]);
        } else {
          xi += wasora_evaluate_expression(material_xs->SigmaA[g]);
          for (g_prime = 0; g_prime < milonga.groups; g_prime++) {
            xi += wasora_evaluate_expression(material_xs->SigmaS0[g][g_prime]);
          }
        }


--
You received this message because you are subscribed to the Google Groups "wasora" group.
To unsubscribe from this group and stop receiving emails from it, send an email to wasora+un...@seamplex.com.
To post to this group, send email to was...@seamplex.com.
Visit this group at https://groups.google.com/a/seamplex.com/group/wasora/.
To view this discussion on the web visit https://groups.google.com/a/seamplex.com/d/msgid/wasora/70ed3ac1a6d84c36a51c3241c18ea76b%40kit-msx-25.kit.edu.
For more options, visit https://groups.google.com/a/seamplex.com/d/optout.

Ferraro, Diego Ernesto (INR)

unread,
Sep 3, 2018, 9:59:25 AM9/3/18
to was...@seamplex.com

I was using SigmaA, I suppose I should use Sigma T then?

Thanks,

Diego

BTW, there is a Manual  “MILONGA version 0.1 , a free nuclear reactor core analysis code”, that I can only find in scribd. Is there any other link available?

César Pablo Camusso

unread,
Sep 3, 2018, 10:06:19 AM9/3/18
to was...@seamplex.com
No, I think that SigmaA and SigmaS are enough.
There are more references in:
https://www.seamplex.com/milonga/publications.html
In this, chapter 5 there are examples and explanations:
https://www.seamplex.com/docs/imef-2013-12-15.pdf

jeremy theler

unread,
Sep 3, 2018, 10:57:51 AM9/3/18
to was...@seamplex.com
On Mon, Sep 3, 2018 at 10:11 AM Ferraro, Diego Ernesto (INR) <diego....@kit.edu> wrote:

Hi guys,

 I’ve a 2 short questions related with units in milonga that I suppose Is quite straightforward to you.

 

1 - I’ve a 3-D problem and I would like to set up the total power of the system to scale neutron fluxes. I understand that I should set the power and eFiss XS as:

MATERIAL xxx {

eSigmaF_1=yyy eSigmaF_2=zzz }

and

power = ppp


yes

 

The problem is that I cannot figure out if the if eSigmaF is really the key I should use and which units I should use


yes
 

à I suppose that the units are something like (per group)  Flux [n/cm2s] = Power[W/cm3] / SigmaF [1/cm] / Efiss [J/fiss] , thus Power should be in Watts and Efiss should be in W per fission (per group). Thus I’m using

eSigmaF_X = SigmaF_X [1/cm] * efiss [J/fiss]  (with power in 

Watts).


correct

what defines the units is the mesh, in the following way:
(see table 3.1 in page 46 in the old milonga v0.1 doc):

All XSs are in 1/length
D is in length
eSigmaF is in energy/length

power depends on the dimension of the problem:
in 3d it is energy/time
in 2d it is energy/(time length)
in 1d it is energy/(time length^2)

and fluxes are in neutrons/(time length^2)

choose your units as you wish, my recommendation is
length = cm
time = s
energy = joules


Anyway, is there any input card reference where I could get the available options for MATERIAL, power, etc (I cannot find it in the repository)?


no, recall that I do not get paid for this and I used to have a regular family with kids and stuff
now I am not even doing any paid work, so this will not be done in any nearly future

I used to think that maybe the community would take over (nerdy IB students with no kids and family yet?)
But I see milonga being highly under-used because I cannot transmit the spirit and the design basis, like when I showed how to analyze the transition between a cube and a sphere... I think there is no other piece of neutronic software that allows users to do such analysis...
 

 

2-I have some doubts with the scattering matrix. I’m using the P0 condensed data for the scattering matrix, but I suppose that the self-scattering (i.e. diagonal) should be transport corrected. So I just changed to P1 for the diagonal.  Nevertheless, when I change the S1.1 and S2.2 values I see no effect in keff. I’m missing something?


I did not follow you this time.
The self-scattering usually appears in both sides of the equation and thus can be discarded, I do not remember the actual equations but all the benchmarks I tested were fine.

Do you have MWE?
 

 

--
You received this message because you are subscribed to the Google Groups "wasora" group.
To unsubscribe from this group and stop receiving emails from it, send an email to wasora+un...@seamplex.com.
To post to this group, send email to was...@seamplex.com.
Visit this group at https://groups.google.com/a/seamplex.com/group/wasora/.
To view this discussion on the web visit https://groups.google.com/a/seamplex.com/d/msgid/wasora/70ed3ac1a6d84c36a51c3241c18ea76b%40kit-msx-25.kit.edu.
For more options, visit https://groups.google.com/a/seamplex.com/d/optout.
--
--
jeremy theler
www.seamplex.com

jeremy theler

unread,
Sep 3, 2018, 11:01:07 AM9/3/18
to was...@seamplex.com
You can use either SigmaA and SigmaS or SigmaT
That document is out of date, but probably some information is still valid.
I wish I had time and skills to update it (my original plan was to have a working startup which allowed me to make a living as well as being able do that stuff for free, but all my plans changed some weeks ago)

As I told you, read the part around table 3.1 in page 46 of that document

--
--
jeremy theler
www.seamplex.com

jeremy theler

unread,
Sep 3, 2018, 11:30:09 AM9/3/18
to was...@seamplex.com
I added that old 2011 doc in markdown (converted from the old latex source) in the doc/ dir of the repo:



Would you write
 1. a list of topics you might want to have addressed in milonga's documentation?
 2. the documentation that addresses the issues in 1?

--
--
jeremy theler
www.seamplex.com

jeremy theler

unread,
Sep 3, 2018, 12:44:27 PM9/3/18
to was...@seamplex.com
one more thing to keep in mind, if you give esigmaf as the prompt energy, power ought to be the prompt power

if esigmaf is the total energy, power ough to be (yeah you guessed) the total power

this is for the steady-state case, I do not know what pablo took into account for the transient case

César Pablo Camusso

unread,
Sep 3, 2018, 1:46:32 PM9/3/18
to was...@seamplex.com
For the transient I used the prompt power.
I had to add a new variable: current_power which is the power as function of time.

jeremy theler

unread,
Sep 3, 2018, 5:19:03 PM9/3/18
to was...@seamplex.com
I would leave it as "power" and use the initial value "power_0" for the initial power.
For each variable defined in wasora (say a), it defines two more:

a_initial is the value of a when t=0 (in_static) and step_static = 1 (in_static_first)
a_0 is the value of a when t=0 (in_static) and step_static = static_steps (in_static_last)

in steady-state, you can set either power or power_0 or power_initial

In any case, a way of computing the delay power should be added. Probably directly as a sum of convolution integrals between the prompt power and N exponentials of time each one with a particular decay constant (not to confuse with the delayed neutron precursor concentration, but the equations are the same).

Different situations arise if eSigmaF is just the prompt or the total energy in the transient case. I have still not thought about how to handle them.

Any thoughts?




To view this discussion on the web visit https://groups.google.com/a/seamplex.com/d/msgid/wasora/1320720072.402736.1535996784749%40mail.yahoo.com.

For more options, visit https://groups.google.com/a/seamplex.com/d/optout.

Ferraro, Diego Ernesto (INR)

unread,
Sep 4, 2018, 5:08:42 AM9/4/18
to was...@seamplex.com

Thanks guys for all the answers!

Just a couple of comments:

1-      I’ve started with a “critical” calculation for my 3x3 core. This configuration should be critical (I’m using 2 group XS obtained from the same 3D model Serpent, where I include Sa, nuSf, Diff and Scattering Matrix with P1 in the diagonal and vacuum condition in the outside). When I analyze the results I get keff 1.00917951 (+900 pcm) which is quite good, but the flux shape is quite strange (some negative values in the top/bottom, some strange results in the central zone which I cannot figure out where came from). I attach 3 plots (the geometry, the boundaries and the flux g1). Any typical hint of what I’m doing wrong?

2-      Regarding the transient, my idea is just to adjust the absortion XS in the central FA as a function of time. Nevertheless I will have to depart from a critical initial state. How should I deal with that?

3-      Regarding the delayed power, I cannot figure out how to deal with that. Usually total power is considered and everything is just scaled to this initial value. But that’s probable not a very good approximation.

Thanks in advance

Diego

geometry.png
outside.png
flux.png

jeremy theler

unread,
Sep 4, 2018, 6:01:02 AM9/4/18
to was...@seamplex.com
try first-order elements first

To view this discussion on the web visit https://groups.google.com/a/seamplex.com/d/msgid/wasora/c4f4e2ba78544186a14e4d27a71b729f%40kit-msx-25.kit.edu.

For more options, visit https://groups.google.com/a/seamplex.com/d/optout.

César Pablo Camusso

unread,
Sep 4, 2018, 8:00:10 AM9/4/18
to was...@seamplex.com
Dear German,
I have added an issue with your advice (power_0 and power_initial).
I used only one power which is calculated as the benchmarks ask.

Dear Diego,
1) When I have those kind of problems, usually the problem is in the mesh or geometry. You can check it by putting the same material in all the physical volumes, 2 groups and null or void in all the faces, so you should get a cosine in the 3 directions. You can see it with plot over line in paraview.
You also can try finite volumes because milonga checks the neighbors and finds mistakes in the mesh.

2) You have to calculate keff, then divide the fission cross section into keff. At this point it is critical, so you have to change the properties of the material of the control rod as function of time.

I did it so:

static_steps = 2 #two static steps.
.
.
.
.
mik_init = 1.0  # initial keff
IF in_static
mik = keff
ENDIF
.
.
.
#The nuSigmaF are divided into mik.
MATERIAL centre {
 D_1   1.0   SigmaA_1   0.02-0.01     nuSigmaF_1  0.005/mik SigmaS_1.2   0.01     eSigmaF_1  0.005  vel_1 1e7
 D_2   0.5   SigmaA_2   0.08          nuSigmaF_2  0.099/mik                       eSigmaF_2  0.099  vel_2 3e5
}

#The change of properties as function of time:
MATERIAL zone1 {
 D_1   1.5   SigmaA_1   0.026-0.015   nuSigmaF_1   0.01/mik  SigmaS_1.2   0.015   eSigmaF_1  0.01   vel_1 1e7
 D_2   0.5   SigmaA_2   SigA(t)       nuSigmaF_2   0.20/mik                       eSigmaF_2  0.20   vel_2 3e5
}

#SigA(t) was defined before:
FUNCTION SigA(t) INTERPOLATION linear DATA {
0     0.18*1
1     0.18*1.03
2     0.18*1.03
}

3) I think you should calculate the power in the same way in both codes to compare the same things.


Ferraro, Diego Ernesto (INR)

unread,
Sep 4, 2018, 10:25:02 AM9/4/18
to was...@seamplex.com

Thanks, I’ll check that.

Actually I changed to volumes and then I got some mismatch in faces, which could explain the issues. But no worries, I’ll fix that.

Diego

Hector Lestani

unread,
Sep 4, 2018, 3:39:45 PM9/4/18
to was...@seamplex.com
With respect to your point number 3....
I think your approach (direct scaling of total power linear to neutron flux) is perfect at any power level BUT some time AFTER achieving the new steady state. During the transient, however, the prompt power follows the neutron flux and the delayed power is related to the power history, so you should model the concentration of the isotopes whose radioactive decay produces power.
As German said, that (delayed power) leads to a sum of exponentials (each one representing the lifetime of an isotope or a group of similar isotopes).
I hope this means something useful for you.

I don't understand your 2nd question, perhaps Pablo does...

Regards.

jeremy theler

unread,
Sep 4, 2018, 4:46:14 PM9/4/18
to was...@seamplex.com
On Tue, 2018-09-04 at 16:37 -0300, Hector Lestani wrote:
> With respect to your point number 3....
> I think your approach (direct scaling of total power linear to
> neutron flux) is perfect at any power level BUT some time AFTER
> achieving the new steady state. During the transient, however, the
> prompt power follows the neutron flux and the delayed power is
> related to the power history, so you should model the concentration
> of the isotopes whose radioactive decay produces power.
> As German said, that (delayed power) leads to a sum of exponentials
> (each one representing the lifetime of an isotope or a group of
> similar isotopes).
> I hope this means something useful for you.

well, no need to do actually that...
the way to handle decay power I always have seen in different codes is
more a fit than a physical model. One gives a set of N (most of the
times N=6 but this has nothing to do with neutron precursor groups)
pairs (amplitude, decay constant) and then convolute the prompt power
history with exponential kernels (a differential way of coding this is
more efficient). It is just a matter of having the right set of
amplitudes and constants. Just to make you a picture, the sum of all
the amplitudes should given approx 7% (it depends again if you take the
rates of the delayed power over the prompt power or over the total
power, apparently there is no consensus on this issue and each code
makes its own choice).

>
> I don't understand your 2nd question, perhaps Pablo does...

I do. I think that the "division" of nuSigmaF by keff should be
automatically done in the code (maybe through a keyword?) and not by
the user. Another option (which I never actually saw in a code) is to
use the alpha-associated reactor critical instead of the keff-
associated one.


> > > 1- I’ve started with a “critical” calculation for my 3x3
> > > core. This configuration should be critical (I’m using 2 group XS
> > > obtained from the same 3D model Serpent, where I include Sa,
> > > nuSf, Diff and Scattering Matrix with P1 in the diagonal and
> > > vacuum condition in the outside). When I analyze the results I
> > > get keff 1.00917951 (+900 pcm) which is quite good, but the flux
> > > shape is quite strange (some negative values in the top/bottom,
> > > some strange results in the central zone which I cannot figure
> > > out where came from). I attach 3 plots (the geometry, the
> > > boundaries and the flux g1). Any typical hint of what I’m doing
> > > wrong?

I agree this is mostly due to the mesh. Do the negative fluxes appear
at t=0 or after the transient. If it is the former, it shoudl be harder
to debug. If it is the latter, I bet the 2nd order elements introduce
oscillations in the time integration.

Pablo, I understand you used PETSc's TS object. Are there any up-winded
or spatially-smoother integrators that reduce oscillations? Can one
change the TS integrator from the input file?

BTW, I think it is time to merge the transient code into the main
branch.

> > > 2- Regarding the transient, my idea is just to adjust the
> > > absortion XS in the central FA as a function of time.
> > > Nevertheless I will have to depart from a critical initial state.
> > > How should I deal with that?

The usual way is to divide nuSigmaF by keff at t=0, but this should be
done directly from the code. Pablo, any thought?

> > > 3- Regarding the delayed power, I cannot figure out how to
> > > deal with that. Usually total power is considered and everything
> > > is just scaled to this initial value. But that’s probable not a
> > > very good approximation.

I did not get this point.

jeremy


César Pablo Camusso

unread,
Sep 5, 2018, 9:15:31 AM9/5/18
to was...@seamplex.com
Dears, I answer what I know:

Pablo, I understand you used PETSc's TS object. Are there any up-winded
or spatially-smoother integrators that reduce oscillations?
There are two oscillations which I found:
1) Spatially the precursors in 2 order elements. It hapens because I programmed it in a way in which some components of the null space of B corrupt the solution (manual of sleps 3.4.4). The advice to solve it is in the paper http://www.ams.org/journals/mcom/1997-66-218/S0025-5718-97-00844-2/S0025-5718-97-00844-2.pdf. It is only to calculate without the precursors in the eigenvalue problem and then calculate the transient.
I am working in it.

2) Oscilations because of the time step is big, the method is not implicit and that kind of thing.
----------------------------------------------
 Can one
change the TS integrator from the input file?
Yes, with:
MILONGA_SOLVER TS_TYPE beuler #backwards Euler
MILONGA_SOLVER TS_TYPE cn # Crank Nicolson
Not all the methods in http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/TS/TSType.html work.
----------------------------------------------

BTW, I think it is time to merge the transient code into the main
branch.

I agree, but the main drawback is that I deplicated some functions of milonga and added them a "_tra" behind. If I keep it, we will have double work to mantain the functions. I do not want to pollute the original version.
If you want, I can send the pull request today or a pull branch or something like that.
----------------------------------------------

The usual way is to divide nuSigmaF by keff at t=0, but this should be
done directly from the code. Pablo, any thought?

Yes.
At the beginning I prefered to touch as low as posible the original code.
I have added other issue with it.
--
You received this message because you are subscribed to the Google Groups "wasora" group.
To unsubscribe from this group and stop receiving emails from it, send an email to wasora+unsub...@seamplex.com.
To post to this group, send email to was...@seamplex.com.

jeremy theler

unread,
Sep 5, 2018, 7:38:56 PM9/5/18
to was...@seamplex.com
On Wed, 2018-09-05 at 13:15 +0000, 'César Pablo Camusso' via wasora
wrote:
> Dears, I answer what I know:
>
> Pablo, I understand you used PETSc's TS object. Are there any up-
> winded
> or spatially-smoother integrators that reduce oscillations?
> There are two oscillations which I found:
> 1) Spatially the precursors in 2 order elements. It hapens because I
> programmed it in a way in which some components of the null space of
> B corrupt the solution (manual of sleps 3.4.4). The advice to solve
> it is in the paper
> http://www.ams.org/journals/mcom/1997-66-218/S0025-5718-97-00844-2/S0025-5718-97-00844-2.pdf

Do the precursor equations involve spatial operators?
I am not sure, but I recall that only the neutron flux is affected by
the spatial derivatives.

Do you have one row of the matrices for each precursor n=1,...,N group
in each node k=1,...K adding N*K degrees of freedom to the problem
size?

> . It is only to calculate without the precursors in the eigenvalue
> problem and then calculate the transient.
> I am working in it.

This is what I would have done in the first place.
Actually, I would have solved the steady-state at t=0 as the base
milonga code (probably taking into account xenon poissoning setting
static_steps > 1) and then advance with TS just the neutron fluxes.
Afterward, I would explcitly compute the N*K precursor concentrations
using a first-order lag like the lag() function in
wasora/src/builtinfunctions.c:766 (actually I "stole" that trick with
the exponential from DYNETZ) over the local neutron flux (and the value
from the last step)

> 2) Oscilations because of the time step is big, the method is not
> implicit and that kind of thing.
> ----------------------------------------------
> Can one
> change the TS integrator from the input file?
> Yes, with:
> MILONGA_SOLVER TS_TYPE beuler #backwards Euler
> MILONGA_SOLVER TS_TYPE cn # Crank Nicolson
> Not all the methods in
> http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/TS/TSType.html
> work.
> ----------------------------------------------

isn't backwards euler implicit and thus inconditionally stable?

> BTW, I think it is time to merge the transient code into the main
> branch.
>
> I agree, but the main drawback is that I deplicated some functions of
> milonga and added them a "_tra" behind. If I keep it, we will have
> double work to mantain the functions. I do not want to pollute the
> original version.
> If you want, I can send the pull request today or a pull branch or
> something like that.
> ----------------------------------------------

add an "if t>0" or "if in_static==false" to the original function and
keep the original code in the "else" condition.

>
> The usual way is to divide nuSigmaF by keff at t=0, but this should
> be
> done directly from the code. Pablo, any thought?
>
> Yes.
> At the beginning I prefered to touch as low as posible the original
> code.
> I have added other issue with it.

ok

--
jeremy

César Pablo Camusso

unread,
Sep 6, 2018, 1:51:05 PM9/6/18
to was...@seamplex.com
Dear German,
I answer shortly:

Do the precursor equations involve spatial operators?
No.
------------------------------------------

Do you have one row of the matrices for each precursor n=1,...,N group
in each node k=1,...K adding N*K degrees of freedom to the problem
size?
Yes.
------------------------------------------
isn't backwards euler implicit and thus inconditionally stable?
Yes.



--
You received this message because you are subscribed to the Google Groups "wasora" group.
To unsubscribe from this group and stop receiving emails from it, send an email to wasora+unsub...@seamplex.com.
To post to this group, send email to was...@seamplex.com.
Visit this group at https://groups.google.com/a/seamplex.com/group/wasora/.

jeremy theler

unread,
Sep 6, 2018, 5:42:15 PM9/6/18
to was...@seamplex.com
On Thu, 2018-09-06 at 17:49 +0000, 'César Pablo Camusso' via wasora
wrote:
> Dear German,
> I answer shortly:
>
> Do the precursor equations involve spatial operators?
> No.
> ------------------------------------------
> Do you have one row of the matrices for each precursor n=1,...,N
> group
> in each node k=1,...K adding N*K degrees of freedom to the problem
> size?
> Yes.

I do not see any gain in having a giant stiffness matrix just to
advance simple exponentials in time. You are adding N*K rows (and thus
columns) to a matrix of size G*K (G is the number of groups) with no
real gain. Think that time and memory scale as the square of the matrix
size.

Can you advance just the flux part with TS and compute the precursor
concentration using first-order lags?

> ------------------------------------------
> isn't backwards euler implicit and thus inconditionally stable?
> Yes.

so why are there these oscilattions?



César Pablo Camusso

unread,
Sep 7, 2018, 8:31:31 AM9/7/18
to was...@seamplex.com
Dear German,
More or less I can answer some things.

I do not see any gain in having a giant stiffness matrix just to
advance simple exponentials in time. You are adding N*K rows (and thus
columns) to a matrix of size G*K (G is the number of groups) with no
real gain. Think that time and memory scale as the square of the matrix
size.

Can you advance just the flux part with TS and compute the precursor
concentration using first-order lags?

I agree. I did so because it was easier to program.
Both can be solved with TS because TS can solve the flux implicitly and the precursors explicitly.
As it is now, TS solve flux and precursors implicitly and takes too time. That is why I did not calculate 3D cases.
---------------------------------------------------------------------------------

so why are there these oscilattions?
I do not know why there were oscilations nor if they were in backwards Euler.
---------------------------------------------------------------------------------
This weekend I am going to look for old files because I wrote the reasons of some decisions.
I cannot answer all the question now.
However I was aware that it is first iteration and it was going to have several things to improve.





--
You received this message because you are subscribed to the Google Groups "wasora" group.
To unsubscribe from this group and stop receiving emails from it, send an email to wasora+unsub...@seamplex.com.
To post to this group, send email to was...@seamplex.com.
Visit this group at https://groups.google.com/a/seamplex.com/group/wasora/.

jeremy theler

unread,
Sep 7, 2018, 1:24:36 PM9/7/18
to was...@seamplex.com
On Fri, 2018-09-07 at 12:31 +0000, 'César Pablo Camusso' via wasora
wrote:
> Can you advance just the flux part with TS and compute the precursor
> concentration using first-order lags?
>
> I agree. I did so because it was easier to program.
> Both can be solved with TS because TS can solve the flux implicitly
> and the precursors explicitly.

But in this way you might end-up with different times for each part of
the problem, won't you?

> As it is now, TS solve flux and precursors implicitly and takes too
> time. That is why I did not calculate 3D cases.

Yes, I can imagine. Start with the precursors in the equilibrium
concentration and code the solution of just the flux part. Then advance
the precursors with explicit first-order lags like in wasora's
builtin_lag() function.

jeremy

Reply all
Reply to author
Forward
0 new messages