More question for transients

12 views
Skip to first unread message

Ferraro, Diego Ernesto (INR)

unread,
Sep 17, 2018, 10:11:22 AM9/17/18
to wasora

Hi guys!

As promised, I continued with my intention of model the transients and I have new doubts J !

I first made several static (criticality) calculations extracting the CR for my 3x3 3D minicore case. When comparing my results with the same from Serpent everything looks OK, so I suppose that there are no big errors in the geometry, XS and so.

 

So encouraged by this, I’ve tried to jump to the next step (basically to put a SigmaA_1 and SigmaA_2 as function of T with the XS from above), using milonga transient capabilities (from https://bitbucket.org/iepale/examples/src).  After some errors I could figure out, I get to one I cannot understand:

 

../../milotra/milo-tra/milonga 3x3.mil

error: gsl error #4 'second index out of range' in ../gsl/gsl_matrix_double.h

error: unspecified error

 

Any ideas? Something I’m missing? I suppose that some of the cards are just wrong, but I cannot figure out which one. I’m attaching the inputs just in case.

Thanks in advance,

Diego

 

 

3x3.geo
3x3.mil

César Pablo Camusso

unread,
Sep 17, 2018, 1:07:04 PM9/17/18
to was...@seamplex.com
Dear Diego,
I saw this line:
PHYSICAL_ENTITY NAME  external BC black
Try using:
PHYSICAL_ENTITY NAME  external BC null
Because I do not remember the black option. I have already done it and the new error mesage is:
error: Velocity <=0 for material 'FA_unrodded' in group 1

It is because you have to set the velocities in all the materials with  vel_1 and vel_2. I added this in your file as example.
You had a mistake in the SigA1(t) function, I corrected it.
I also had to add the keff calculation in static because I have not corrected it yet. However it is an issue and it is writen.
I added other variable in the print because power is the setpoint; whereas current_power is the current power. It also can change.

I am working slowly in correcting several things that we talked some weeks ago.
I also changed the geo because it uses too memory.
Check it, it runs:
pablo@debian:~/Descargas$ ~/petsc/wasora-suite/milo-tra/milonga 3x3.mil
# keff    t    current_power
1.0695419499    0.0000000000    140900000.0000000000
1.0000000432    0.0000000000    140900000.0000000000
1.0000000432    0.0100000000    139138644.5529071093
1.0000000432    0.0200000000    139254002.2049786747
1.0000000432    0.0300000000    139060207.6668723822
1.0000000432    0.0400000000    139258278.4846861660
1.0000000432    0.0500000000    139093852.3822170496
1.0000000432    0.0600000000    139271307.4797407389
1.0000000432    0.0700000000    139126149.3733670115
1.0000000432    0.0800000000    139285046.7580571175
1.0000000432    0.0900000000    139156750.8754683733
1.0000000432    0.1000000000    139299330.6523751915
1.0000000432    0.1100000000    139185827.2395524681
1.0000000432    0.1200000000    139314049.5939984322
1.0000000432    0.1300000000    139213530.8341576159
1.0000000432    0.1400000000    139329072.5328560174
1.0000000432    0.1500000000    139239993.1129422486
1.0000000432    0.1600000000    139344384.7020989060
1.0000000432    0.1700000000    139265375.2067396939
.
.
.
.









--
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/4b73dc591e5345d2bfe77f490623f573%40kit-msx-25.kit.edu.
For more options, visit https://groups.google.com/a/seamplex.com/d/optout.
3x3changed.zip

César Pablo Camusso

unread,
Sep 17, 2018, 2:32:08 PM9/17/18
to was...@seamplex.com
I had to add other change because it did not work:

FUNCTION SigA1(t) INTERPOLATION linear DATA {
0.0        1.60947E-03
0.200    0.0/40.0*1.05658E-02+(40-0)/40.0*1.60947E-03
0.325    5.0/40.0*1.05658E-02+(40-5)/40.0*1.60947E-03
0.450   10.0/40.0*1.05658E-02+(40-10)/40.0*1.60947E-03
0.575   15.0/40.0*1.05658E-02+(40-15)/40.0*1.60947E-03
0.700   20.0/40.0*1.05658E-02+(40-20)/40.0*1.60947E-03
0.825   25.0/40.0*1.05658E-02+(40-25)/40.0*1.60947E-03
0.950   30.0/40.0*1.05658E-02+(40-30)/40.0*1.60947E-03
1.075   35.0/40.0*1.05658E-02+(40-35)/40.0*1.60947E-03
1.200   40.0/40.0*1.05658E-02+(40-40)/40.0*1.60947E-03
1.325   35.0/40.0*1.05658E-02+(40-35)/40.0*1.60947E-03
1.450   30.0/40.0*1.05658E-02+(40-30)/40.0*1.60947E-03
1.575   25.0/40.0*1.05658E-02+(40-25)/40.0*1.60947E-03
1.700   20.0/40.0*1.05658E-02+(40-20)/40.0*1.60947E-03
1.825   15.0/40.0*1.05658E-02+(40-15)/40.0*1.60947E-03
1.950   10.0/40.0*1.05658E-02+(40-10)/40.0*1.60947E-03
2.075     5.0/40.0*1.05658E-02+(40-5)/40.0*1.60947E-03
2.200     0.0/40.0*1.05658E-02+(40-0)/40.0*1.60947E-03
}


FUNCTION SigA2(t) INTERPOLATION linear DATA {
0.0         4.67323E-02
0.200    0.0/40.0*1.21032E-01+(40-0)/40.0*4.67323E-02
0.325    5.0/40.0*1.21032E-01+(40-5)/40.0*4.67323E-02
0.450   10.0/40.0*1.21032E-01+(40-10)/40.0*4.67323E-02
0.575   15.0/40.0*1.21032E-01+(40-15)/40.0*4.67323E-02
0.700   20.0/40.0*1.21032E-01+(40-20)/40.0*4.67323E-02
0.825   25.0/40.0*1.21032E-01+(40-25)/40.0*4.67323E-02
0.950   30.0/40.0*1.21032E-01+(40-30)/40.0*4.67323E-02
1.075   35.0/40.0*1.21032E-01+(40-35)/40.0*4.67323E-02
1.200   40.0/40.0*1.21032E-01+(40-40)/40.0*4.67323E-02
1.325   35.0/40.0*1.21032E-01+(40-35)/40.0*4.67323E-02
1.450   30.0/40.0*1.21032E-01+(40-30)/40.0*4.67323E-02
1.575   25.0/40.0*1.21032E-01+(40-25)/40.0*4.67323E-02
1.700   20.0/40.0*1.21032E-01+(40-20)/40.0*4.67323E-02
1.825   15.0/40.0*1.21032E-01+(40-15)/40.0*4.67323E-02
1.950   10.0/40.0*1.21032E-01+(40-10)/40.0*4.67323E-02
2.075    5.0/40.0*1.21032E-01+(40-5)/40.0*4.67323E-02
2.200    0.0/40.0*1.21032E-01+(40-0)/40.0*4.67323E-02
}




El lunes, 17 de septiembre de 2018 11:12:25 a. m. GMT-3, Ferraro, Diego Ernesto (INR) <diego....@kit.edu> escribió:


--

jeremy theler

unread,
Sep 17, 2018, 4:36:23 PM9/17/18
to was...@seamplex.com
On Mon, 2018-09-17 at 14:11 +0000, Ferraro, Diego Ernesto (INR) wrote:
> Hi guys!
> As promised, I continued with my intention of model the transients
> and I have new doubts J !
> I first made several static (criticality) calculations extracting the
> CR for my 3x3 3D minicore case. When comparing my results with the
> same from Serpent everything looks OK, so I suppose that there are no
> big errors in the geometry, XS and so.

great job for the results, bad job for the presentation...
a plot in a bitmap format? you would be fired by now if you worked at
seamplex :-)

>
> So encouraged by this, I’ve tried to jump to the next step (basically
> to put a SigmaA_1 and SigmaA_2 as function of T with the XS from
> above), using milonga transient capabilities (from
> https://bitbucket.org/iepale/examples/src). After some errors I
> could figure out, I get to one I cannot understand:
>
> ../../milotra/milo-tra/milonga 3x3.mil
> error: gsl error #4 'second index out of range' in
> ../gsl/gsl_matrix_double.h
> error: unspecified error

this is an error from GSL that wasora could not catch so it passed the
acutal error to the user... I know this is cryptic, but it is a last
resource...

so I do not know what triggered that (probably something in Pablo's
code) but here are my contributions

1. some time ago, FUNCTIONSs could be point-defined with DATA involving
expresions in variables, but then I had to re-write some stuff and
could not handle this anymore (does anybody remember why?). So DATA can
involve expressions but only involing numbers, not variables. In any
case, mind the spaces. The parser looks for blank spaces separating
tokens so spaces tabs and newlines are counted as token separators
(because you are within a line continuation block with the bracket '{'
after DATA) so the expressions cannot contain blanks. At the parser
level, your input starts to shift from independent to dependent
variable values because there are some spaces in the expressions and
that can be a source of GSL error.

2. I would re-write the time-dependent functions this way (it works and
it is easier to understand):

FUNCTION theta(t) DATA {
0.0 0
0.200 0/40
0.325 5/40
0.450 10/40
0.575 15/40
0.700 20/40
0.825 25/40
0.950 30/40
1.075 35/40
1.200 40/40
1.325 35/40
1.450 30/40
1.575 25/40
1.700 20/40
1.825 15/40
1.950 10/40
2.075 5/40
2.200 0/40
}

SigA1(t) := theta(t)*SigmaA1out + (1-theta(t))*SigmaA1in
SigA2(t) := theta(t)*SigmaA2out + (1-theta(t))*SigmaA2in


In this case, I get negative values though but probably I got the
values wrong. See the attached sigmas.was
Maybe that is also triggering another error in milonga.

> Any ideas? Something I’m missing? I suppose that some of the cards

This term dates back to 1950s where input was written as punched cards.
Please do not use it for software designed and written after 1990.

> are just wrong, but I cannot figure out which one. I’m attaching the
> inputs just in case.
> Thanks in advance,
> Diego
>
>
sigmas.was

jeremy theler

unread,
Sep 17, 2018, 4:45:40 PM9/17/18
to was...@seamplex.com
On Mon, 2018-09-17 at 17:07 +0000, 'César Pablo Camusso' via wasora
wrote:
> Dear Diego,
> I saw this line:
> PHYSICAL_ENTITY NAME external BC black
> Try using:
> PHYSICAL_ENTITY NAME external BC null
> Because I do not remember the black option. I have already done it

There is no "black" option. Milonga does not complain (probably we can
make it complain), it just uses the default BC which is vacuum I think.

> and the new error mesage is:
> error: Velocity <=0 for material 'FA_unrodded' in group 1
>
> It is because you have to set the velocities in all the materials
> with vel_1 and vel_2. I added this in your file as example.

can you set a default value? or compute the velocities from the energy
at which the group changes?
remember that the printf format in PRINT can be changed for each
printed expression, so for a clearer output use something like

PRINT %.4f (keff–1)/keff %.2f t %.5e current_power HEADER

>
>
>
>
>
>
>
>
> El lunes, 17 de septiembre de 2018 11:12:25 a. m. GMT-3, Ferraro,
> Diego Ernesto (INR) <diego....@kit.edu> escribió:
>
>
> Hi guys!
> As promised, I continued with my intention of model the transients
> and I have new doubts J !
> I first made several static (criticality) calculations extracting the
> CR for my 3x3 3D minicore case. When comparing my results with the
> same from Serpent everything looks OK, so I suppose that there are no
> big errors in the geometry, XS and so.
>
>

Ferraro, Diego Ernesto (INR)

unread,
Sep 18, 2018, 9:51:26 AM9/18/18
to was...@seamplex.com

Thanks guys!

A couple of comments:

1-      Regarding requirements: If I maintain the original ~3cm size of the mesh, I get ~10Gb of RAM. Also it takes several minutes per step.

2-      I used the approach:

CONST SigmaA1out SigmaA2out SigmaA1in SigmaA2in nuSigmaF1in  nuSigmaF2in  nuSigmaF1out  nuSigmaF2out

SigmaA1out = 1.05658E-02

SigmaA2out = 1.21032E-01

SigmaA1in = 1.40071E-02

SigmaA2in = 1.49093E-01

nuSigmaF1in = 8.62877E-03

nuSigmaF2in = 1.88954E-01

nuSigmaF1out = 8.77546E-03

nuSigmaF2out = 1.84418E-01

 

 

FUNCTION theta(t) DATA {

0.0     0

 0.200       0/40 

 0.325   5/40

0.450   10/40

0.575   15/40

0.700   20/40

0.825   25/40

0.950   30/40

1.075   35/40

1.200   40/40

1.325   35/40

1.450   30/40

1.575   25/40

1.700   20/40

1.825   15/40

1.950   10/40

2.075   5/40

2.200       0/40

}

 SigA1(t) := theta(t)*SigmaA1out + (1-theta(t))*SigmaA1in

SigA2(t) := theta(t)*SigmaA2out + (1-theta(t))*SigmaA2in

nuSigF1(t) := theta(t)*nuSigmaF1out + (1-theta(t))*nuSigmaF1in

nuSigF2(t) := theta(t)*nuSigmaF2out + (1-theta(t))*nuSigmaF2in

 

3-      But for some reason, when I do PRINT %.2f t %.10f keff  %.5e current_power %.4f SigA1(t) SigA2(t) nuSigF1(t) nuSigF2(t) HEADER  the keff are not correct:

0.00 1.0000000634    1.40900e+08     0.0140     0.1491     0.0086     0.1890

0.10 1.0000000634    1.39977e+08     0.0140     0.1491     0.0086     0.1890

0.20 1.0000000634    1.39295e+08     0.0140     0.1491     0.0086     0.1890

0.30 1.0000000634    1.39028e+08     0.0137     0.1463     0.0086     0.1885

0.40 1.0000000634    1.39197e+08     0.0133     0.1435     0.0087     0.1880

0.50 1.0000000634    1.39636e+08     0.0130     0.1407     0.0087     0.1876

0.60 1.0000000634    1.40482e+08     0.0126     0.1379     0.0087     0.1871

0.70 1.0000000634    1.41610e+08     0.0123     0.1351     0.0087     0.1867

0.80 1.0000000634    1.43182e+08     0.0119     0.1323     0.0087     0.1862

0.90 1.0000000634    1.45107e+08     0.0116     0.1295     0.0087     0.1858

1.00 1.0000000634    1.47582e+08     0.0113     0.1266     0.0087     0.1853

1.10 1.0000000634    1.50556e+08     0.0109     0.1238     0.0088     0.1849

1.20 1.0000000634    1.54282e+08     0.0106     0.1210     0.0088     0.1844

1.30 1.0000000634    1.56783e+08     0.0109     0.1238     0.0088     0.1849

1.40 1.0000000634    1.58451e+08     0.0113     0.1266     0.0087     0.1853

1.50 1.0000000634    1.59297e+08     0.0116     0.1295     0.0087     0.1858

1.60 1.0000000634    1.59652e+08     0.0119     0.1323     0.0087     0.1862

1.70 1.0000000634    1.59449e+08     0.0123     0.1351     0.0087     0.1867

1.80 1.0000000634    1.58989e+08     0.0126     0.1379     0.0087     0.1871

1.90 1.0000000634    1.58155e+08     0.0130     0.1407     0.0087     0.1876

4-      The results are somehow underestimating the power peak. I’ll check with better discretization, but it looks like it’s another thing. Up to now I get the png attached. I suppose that the last part (after end of movement at 2.2 sec) the  FUNCTION theta(t) DATA is extrapolating the table (that’s easy to correct). But nevertheless I cannot figure out why the peak is so small. Ideas?

Thanks,

Diego

transient.png

César Pablo Camusso

unread,
Sep 18, 2018, 10:11:27 AM9/18/18
to was...@seamplex.com
Hello.

2) I think you should do it:

FUNCTION theta(t) DATA {

0.0     0

 0.200       0/40 

 0.325   5/40

0.450   10/40

0.575   15/40

0.700   20/40

0.825   25/40

0.950   30/40

1.075   35/40

1.200   40/40

1.325   35/40

1.450   30/40

1.575   25/40

1.700   20/40

1.825   15/40

1.950   10/40

2.075   5/40

2.200       0/40

3.0       0

}
So theta ends in 0 and remains 0.

3) Only the first two keff import because it is not anomore calculated. It is only used to get the critical initial condition.

4)In point 2 I added other point because when you extrapolate the returned value is taken from the last rect line; therefore I added a 0 in order to the extrapolated value is 0. For example in your case, whitout the last point, after the 2.2 s you get negative values.

I ran the files that you sent us yesterday and I got that the power decreases and then it increases. It is right because SigA1(t) and SigA2(t) increases and then the power decreases. I attached all (unless the msh file).


3x3.geo
3x3.mil
pow.pdf
SigA.pdf

jeremy theler

unread,
Sep 18, 2018, 1:46:00 PM9/18/18
to was...@seamplex.com
On Tue, 2018-09-18 at 13:51 +0000, Ferraro, Diego Ernesto (INR) wrote:
> Thanks guys!
> A couple of comments:
> 1- Regarding requirements: If I maintain the original ~3cm size
> of the mesh, I get ~10Gb of RAM. Also it takes several minutes per
> step.

are your geometry and materials symmetric? if yes, use 1/8th of the
core, that should significantly reduce the computational effort


> 3- But for some reason, when I do PRINT %.2f t %.10f keff %.5e
> current_power %.4f SigA1(t) SigA2(t) nuSigF1(t) nuSigF2(t) HEADER
> the keff are not correct:

the keff you are printing is the eigen-value of the first static
calculation
I do not think Pablo is computing the kinetic reactivity at each step
of the transient calculation

> 4- The results are somehow underestimating the power peak. I’ll
> check with better discretization, but it looks like it’s another
> thing. Up to now I get the png attached. I suppose that the last part
> (after end of movement at 2.2 sec) the FUNCTION theta(t) DATA is
> extrapolating the table (that’s easy to correct). But nevertheless I
> cannot figure out why the peak is so small. Ideas?
>

no idea about that, but what I would do first is to have milonga to run
a quasi-static case, i.e. the same thing but instead of transient just
to compute the keff (i.e. the static reactivity) at each time step. If
you use the same input with the original milonga it should work this
way (btw, I changed some stuff and added the msh4 parser)

then, with the static reactivity as a function of time, I would solve
the point kinetics equations with wasora using rho(t) to see what is
going on. See

https://bitbucket.org/seamplex/wasora/src/master/examples/reactor.was



Reply all
Reply to author
Forward
0 new messages