Author Topic: Lambda funksjon fra prosjekt 1  (Read 213 times)

Anonymous

  • Guest
Lambda funksjon fra prosjekt 1
« on: 31. March 2021, 01:21:05 »
Jeg får litt problemer med lambdafunksjonen fra prosjekt 1. De forskjellige lambdaene er avhengig av \( T^9 \), mens vi i denne oppgaven kun bruker \( T \). Dette har jeg fikset ved at input temperaturen blir ganget med \( 10^9 \). Da får jeg derimot et problem for de laveste temperaturene. De gir utrolig høye verdier for lambda. Dette gir overflow og runtime warnings. Noen som vet hva som er greia her?

Helle Bakke

  • Global Moderator
  • Member
  • Posts: 75
Re: Lambda funksjon fra prosjekt 1
« Reply #1 on: 31. March 2021, 10:35:48 »
\(T^9\) er avhenging av \(T\), så dersom du sender temperaturen inn i lambda-funksjonen skal det kun være nødvendig å gjøre konverteringen fra \(T\) til \(T^9\). Hvis dette fungerte i prosjekt 1 burde du se over hvordan du gjør konverteringen, og sjekke at det stemmer overens med det du hadde i prosjekt 1 :)

Anonymous

  • Guest
Re: Lambda funksjon fra prosjekt 1
« Reply #2 on: 31. March 2021, 10:48:39 »
Skrev at jeg ganget med \( 10^9 \), mente selvfølgelig at jeg deler. I prosjekt 1 gjør jeg ikke konverteringen i det hele tatt, jeg holder meg til \( T9 \) gjennom hele.

Har dette:

Code: [Select]
def lamfunc(reac,T_conv):
    T9 = T_conv/1e9
    if reac == 'pp':
        return kdelta*conv*(4.01e-15*T9**(-2/3)*np.exp(-3.380*T9**(-1/3))*(1+0.123*T9**(1/3)+1.09*T9**(2/3)+0.938*T9))
    elif reac == '33':
        return kdelta*conv*6.04e10*T9**(-2/3)*np.exp(-12.276*T9**(-1/3))*(1+0.034*T9**(1/3)-0.522*T9**(2/3)-0.124*T9+0.353*T9**(4/3)+0.213*T9**(5/3))
    elif reac == '34':
        return conv*5.61e6*(T9/(1+4.95e-2*T9))**(5/6)*T9**(-3/2)*np.exp(-12.826*(T9/(1+4.95e-2*T9))**(-1/3))
    elif reac == 'e7':
        return conv*(1.34e-10*T9**(-1/2)*(1-0.537*T9**(1/3)+3.86*T9**(2/3)+0.0027*T9**(-1)*np.exp(2.515e-3*T9**(-1))))
    elif reac == '17*':
        return conv*(1.096*10**9*T9**(-2/3)*np.exp(-8.472*T9**(-1/3))-4.830*10**8*(T9/(1+0.759*T9))**(5/6)*T9**(-3/2)*np.exp(-8.472*(T9/(1+0.759*T9))**(-1/3))+1.06*10**10*T9**(-3/2)*np.exp(-30.442*T9**(-1)))
    elif reac == '17':
        return conv*(3.11e5*T9**(-2/3)*np.exp(-10.262*T9**(-1/3))+2.53e3*T9**(-3/2)*np.exp(-7.306*T9**(-1)))
    elif reac == 'p14':
        return conv*(4.9e7*T9**(-2/3) * np.exp(-15.228*T9**(-1/3)-0.092*T9**2) * (1 + 0.027*T9**(1/3) - 0.778*T9**(2/3) - 0.149*T9 + 0.261*T9**(4/3) + 0.127*T9**(5/3)) + 2.37e3*T9**(-3/2) * np.exp(-3.011*T9**(-1)) + 2.19e4*np.exp(-12.53*T9**(-1)))

Så hvis jeg sender inn \( T_{conv} = 8000K \), så gjør den beregningene med \( T9 = 8.0 \cdot 10^{-6} \), hvilket fører til veldig høye lambda-verdier. Hva burde jeg gjøre annerledes?

Helle Bakke

  • Global Moderator
  • Member
  • Posts: 75
Re: Lambda funksjon fra prosjekt 1
« Reply #3 on: 31. March 2021, 11:12:19 »
Pass på at du får med de variasjonene av \(T_9\) som står i tabellen i kapittel 3 (\(T_{9*}\)), jeg ser ikke de umiddelbart i kodesnutten din :)

Anonymous

  • Guest
Re: Lambda funksjon fra prosjekt 1
« Reply #4 on: 31. March 2021, 11:18:59 »
It is included. I use the relation between T9 and T9* as given in the table. The T9*'s are exchanged with this expression directly.

Frederik Clemmensen

  • Global Moderator
  • Member
  • Posts: 127
Re: Lambda funksjon fra prosjekt 1
« Reply #5 on: 31. March 2021, 11:40:00 »
Hvis din kode regner alt med T9, bør du vel give den T9 som input fra starten, ikke kun i utregningen af lambda. Kanske der mangler en omregning et andet sted.

Anonymous

  • Guest
Re: Lambda funksjon fra prosjekt 1
« Reply #6 on: 31. March 2021, 11:47:29 »
So you suggest that I should use T9 all the way through project 2? Won't that mess with the results when everything else is in SI-units?

Frederik Clemmensen

  • Global Moderator
  • Member
  • Posts: 127
Re: Lambda funksjon fra prosjekt 1
« Reply #7 on: 31. March 2021, 12:01:42 »
No, I meant in your project 1. I assumed you used the lamda function from your project1.py (or whatever you called it) directly, but if you moved it to your project2.py file, it should of course be the only place to make that conversion.

Anonymous

  • Guest
Re: Lambda funksjon fra prosjekt 1
« Reply #8 on: 31. March 2021, 12:15:17 »
Yeah, that's what I figured, but the lambda-function does not seem very keen on accepting the temperatures of project 2, even though they are converted to T9 within the function.

Anonymous

  • Guest
Re: Lambda funksjon fra prosjekt 1
« Reply #9 on: 31. March 2021, 13:05:51 »
Okay, I think I solved that problem, there was a problem within a different function. However, I still did not get any plots, as there still is something wrong. Is it possible to get any help with the code over zoom or something like that? I know we should focus on the report at this point, and I have written a big part of it, but I feel like I am so close to at least get some plots so that I have something to write about.

Frederik Clemmensen

  • Global Moderator
  • Member
  • Posts: 127
Re: Lambda funksjon fra prosjekt 1
« Reply #10 on: 31. March 2021, 14:24:57 »
We will not be doing another zoom session, but explain your problem here and we will try to help :)

Anonymous

  • Guest
Re: Lambda funksjon fra prosjekt 1
« Reply #11 on: 31. March 2021, 14:41:23 »
Okay, so I believe my problem is with the dynamic step thingy. I have:

Code: [Select]
for i in range(N-1):
        dm_r = np.abs(p * rad[i] / DEr(rad[i],rho1[i]))
        dm_P = np.abs(p * P[i] / DEP(rad[i],m[i]))
        dm_L = np.abs(p * L1[i] / epsilon(T1[i],rho1[i]))
        dm_T = np.abs(p * T1[i] / DET(L1[i],T1[i],rho1[i],rad[i],P[i],m[i]))
        print(dm_r,dm_P,dm_L,dm_T)
        #dm_T_min = min(dm_T)
        dm_wrong = np.array([dm_r,dm_P,dm_L,dm_T])
        dm = min(dm_wrong.ravel()[np.flatnonzero(dm_wrong)])
        print(dm)

Where the DEr, DEP, epsilon and DET solves the differential equations (1,2,3,5) given in the assignment. They are correct, but for instance r = r0, DEr returns a very small value, as r^2 is in the denominator. Similarly for DEP, r^4 is in the denominator. This makes the first dm-value huge (3.45 e18), and that is after choosing the smallest one. Thus, I get overflow and runtime warnings right away.

Any clues to what I am doing wrong?

Frederik Clemmensen

  • Global Moderator
  • Member
  • Posts: 127
Re: Lambda funksjon fra prosjekt 1
« Reply #12 on: 31. March 2021, 14:56:33 »
That seems reasonable, dm may be large ~\(  10^{18} \), but since \( M_0=1.989 \times10^{30} \) it is still a small change in mass. That should be ok.

Anonymous

  • Guest
Re: Lambda funksjon fra prosjekt 1
« Reply #13 on: 31. March 2021, 15:14:35 »
Okay, then maybe there is a problem with my kappa.

I have:

Code: [Select]
func = interpolate.interp2d(R,T,k)

def infunc(Tg,Rg):
    if Tg < T[0] or Tg > T[-1] or Rg < R[0] or Rg > R[-1]:
        print('Warning, out of bounds')
    return(func(Rg,Tg))

def k(rho,T):
    return (0.1*10**infunc(np.log(T),-np.log(rho/10/((T*1e-6)**3)))[0])

And I immediately get a warning in the function. I have tried the k function with T0 and rho0, which returns k = 1.69e-2 with no warning. Any ideas?

Frederik Clemmensen

  • Global Moderator
  • Member
  • Posts: 127
Re: Lambda funksjon fra prosjekt 1
« Reply #14 on: 31. March 2021, 15:20:39 »
You have a minus in front of your \( \log_{10}R \), when you input it into infunc(). And maybe you already take this into account, but remember that the rho in the table is in cgs.

Anonymous

  • Guest
Re: Lambda funksjon fra prosjekt 1
« Reply #15 on: 31. March 2021, 15:39:42 »
Yeah, okay. If I remove the minus sign, the infunc function returns zero. That's not right... I fixed the conversion now, but it still has the same issue.

Code: [Select]
func = interpolate.interp2d(R,T,k)

def infunc(Tg,Rg):
    if Tg < T[0] or Tg > T[-1] or Rg < R[0] or Rg > R[-1]:
        print('Warning, out of bounds')
    return(func(Rg,Tg))

def k(rho,T):
    return (0.1*10**infunc(np.log(T),-np.log(rho*1e3/((T*1e-6)**3)))[0])

printing for rho0 and T0 gives:

Warning, out of bounds
0.024766274906383903

Frederik Clemmensen

  • Global Moderator
  • Member
  • Posts: 127
Re: Lambda funksjon fra prosjekt 1
« Reply #16 on: 31. March 2021, 15:43:42 »
Oh, also use np.log10, not np.log. That is the natural logarithm.

Frederik Clemmensen

  • Global Moderator
  • Member
  • Posts: 127
Re: Lambda funksjon fra prosjekt 1
« Reply #17 on: 31. March 2021, 15:47:50 »
And \( \rho_{cgs} = \rho_{SI}\times 10^{-3} \).

Anonymous

  • Guest
Re: Lambda funksjon fra prosjekt 1
« Reply #18 on: 31. March 2021, 16:00:09 »
Yup! Well spotted. I still get the same warnings, though.

Warning, out of bounds
Warning, out of bounds
Warning, out of bounds
Warning, out of bounds
Warning, out of bounds
Warning, out of bounds
Warning, out of bounds
RuntimeWarning: divide by zero encountered in double_scalars
  scaling_He = rd['pp']/(2*rd['33']+rd['34'])
RuntimeWarning: overflow encountered in double_scalars
  scaling_Li = rd['e7']/rd['17*']
RuntimeWarning: divide by zero encountered in double_scalars
  T1[i+1] = T1 + 3*k(T1,rho1)*L1/(256*pi**2*sig**rad**4*T1**3)*dm
RuntimeWarning: invalid value encountered in double_scalars
  return mu*m_u*((Press(rho,T)-a/3*T**3)/kb)

I get the 'warning, out of bounds' after running the solver, but I have also told it to print(i) for every iteration, and it doesn't, which makes me think it does not go through it correctly a single time. After these messages, it yields the error message: LinAlgError: Array must not contain infs or NaNs.

Thoughts?

Anonymous

  • Guest
Re: Lambda funksjon fra prosjekt 1
« Reply #19 on: 31. March 2021, 16:01:01 »
But you were right with the kappa thing. I no longer get errors for printing with T0 and rho0.

Frederik Clemmensen

  • Global Moderator
  • Member
  • Posts: 127
Re: Lambda funksjon fra prosjekt 1
« Reply #20 on: 31. March 2021, 16:35:37 »
Good!

Well, if you your program tells you you are dividing by 0, use print statements to check which variable is 0 and why.

Anonymous

  • Guest
Re: Lambda funksjon fra prosjekt 1
« Reply #21 on: 31. March 2021, 16:53:35 »
Code: [Select]
T1[i+1] = T1[i] + 3*k(rho1[i],T1[i])*L1[i]/(256*pi**2*sig**rad[i]**4*T1[i]**3)*dm
        print('here',T[i+1])
        rho1[i+1] = rhofunc(P[i+1],T1[i+1],rho1[i])
        print(rho1[i+1])

Tells me that T[i+1] is 3.8. Which I think maybe should be 10^(3.8) ??
It also returns rho[i+1] = nan

Anonymous

  • Guest
Re: Lambda funksjon fra prosjekt 1
« Reply #22 on: 31. March 2021, 16:54:49 »
Code: [Select]
T1[i+1] = T1[i] + 3*k(rho1[i],T1[i])*L1[i]/(256*pi**2*sig**rad[i]**4*T1[i]**3)*dm
        print('here',T[i+1])
        rho1[i+1] = rhofunc(P[i+1],T1[i+1],rho1[i])
        print(rho1[i+1])

Tells me that T[i+1] is 3.8. Which I think maybe should be 10^(3.8) ??
It also returns rho[i+1] = nan

Anonymous

  • Guest
Re: Lambda funksjon fra prosjekt 1
« Reply #23 on: 31. March 2021, 16:55:19 »
Frickin smileys

Frederik Clemmensen

  • Global Moderator
  • Member
  • Posts: 127
Re: Lambda funksjon fra prosjekt 1
« Reply #24 on: 31. March 2021, 18:20:07 »
I don't know why you would get the result in 10^XX.

I think in general you would benefit from defining more variables instead of calculating everything at once. It makes it much harder to read and easier to make mistakes.

Anonymous

  • Guest
Re: Lambda funksjon fra prosjekt 1
« Reply #25 on: 31. March 2021, 19:13:52 »
Okay, so I am starting to get the hang of this crappy code. I must therefore ask, is it possible that the density, being 0.00019... for rho[0], then changes to rho[1] = 1.164 after the first integration?

Anonymous

  • Guest
Re: Lambda funksjon fra prosjekt 1
« Reply #26 on: 31. March 2021, 19:55:00 »
Okay, so now the code works for N = 1, but it seems that on the second integration, I get a new error. I have been trying to fix it for some time now, without luck. Here is the error message:

~\documents\ast3310\project 2\energy_transport.py in solver()
    359         dm_P = np.abs(p * P / DEP(rad,m))
    360         dm_L = np.abs(p * L1 / epsilon(T1,rho1))
--> 361         dm_T = np.abs(p * T1 / DET(L1,T1,rho1,rad,P,m))
    362         print(dm_r,dm_P,dm_L,dm_T)
    363         #dm_T_min = min(dm_T)

~\documents\ast3310\project 2\energy_transport.py in DET(L, T, rho, r, P, m)
    313 def DET(L,T,rho,r,P,m):
    314     if nabla_stable(L,T,rho,r,m) > nabla_ad:
--> 315         return nabla_star(rho,T,r,m,L)*T/P*DEP(r,m)
    316     else:
    317         return -(3*k(rho,T)*L)/(256*pi**2*sig*r**4*T**3)

TypeError: can't multiply sequence by non-int of type 'numpy.float64'


All help greatly appreciated :))

Frederik Clemmensen

  • Global Moderator
  • Member
  • Posts: 127
Re: Lambda funksjon fra prosjekt 1
« Reply #27 on: 31. March 2021, 20:39:11 »
That error message would suggest you may have a list that is not a numpy array, that you are trying to multiply with a floating number. If that is the case, converting the list to a numpy array should do the trick. Otherwise, it could be that what you expect to be a single number is actually a list.

Anonymous

  • Guest
Re: Lambda funksjon fra prosjekt 1
« Reply #28 on: 31. March 2021, 20:56:57 »
Yeah, I figured the same, but it is strange that the integration works fine for the first iteration, but gets an error message for the second round. Oh well, it is too late now anyways... Thanks, though :)