Imposing scalar at inlet by user subroutine

Questions and remarks about code_saturne usage
Forum rules
Please read the forum usage recommendations before posting.
Post Reply
Tsubasa
Posts: 175
Joined: Fri Nov 20, 2020 4:09 am

Imposing scalar at inlet by user subroutine

Post by Tsubasa »

Hello,

Now I want to impose a value for scalar transport at inlet face by cs_user_boundary_conditions.f90, not GUI.

Can we impose a value of scalar transport at inlet by cs_user_boundary_conditions.f90?
At the same time I am imposing pressure value following the code.
However, errors happend when computing.
Just wIth the code for pressure, it works, but with the code of scalar it does not work.

Can we impose the specific value of scalar transpot at inlet by usersubroutine?
Once I succeed, I will edit the code to impose time-dependent value for scalar transport.

Code: Select all

call getfbr('inlet', nlelt, lstelt)
do ilelt = 1, nlelt
   ifac = lstelt(ilelt)
   itypfb(ifac) = ientre

!Dirichlet B.C for pressure:
   icodcl(ifac,ipr) = 1
   rcodcl(ifac,ipr,1) = 101325.633

!Dirichlet B.C for scalar:
icodcl(ifac,isca(1))=1
rcodcl(ifac,isca(1),1)=10 ! desired value 

!Neumann B.C for velocity:
   icodcl(ifac,iu)= 3
   icodcl(ifac,iv)= 3
   icodcl(ifac,iw)= 3
enddo
In the documentation, I found that "isca(i) is the index of the scalar i".
What is "i"? and Can I find "i" in setup log?
Would someone help me?

Tsubasa
Yvan Fournier
Posts: 4070
Joined: Mon Feb 20, 2012 3:25 pm

Re: Imposing scalar at inlet by user subroutine

Post by Yvan Fournier »

Hello,

I see we still have examples using isca(1) directly, which is unfortunate, as this directly uses a scalar's index (1) instead of accessing it by name (which is safer/clearer/more independent of scalrs creation order).

instread of :

Code: Select all

isca(1)
Use something like:

Code: Select all

integer :: f_id, keysca, iscal
...
call field_get_id(scalar_name, f_id)
call field_get_key_id("scalar_id", keysca)
call field_get_key_int(f_id, keysca, iscal)
Where iscal replaces isca(1) and scalar_name is the actual name given to the scalar.

Best regards,

Yvan
Tsubasa
Posts: 175
Joined: Fri Nov 20, 2020 4:09 am

Re: Imposing scalar at inlet by user subroutine

Post by Tsubasa »

Hello Yvan,

Thank you for the advice.

You mean if I assign "scalar1" as a name of scalar in GUI like in picture,
I should write like this for the first line, right?

Code: Select all

integer :: f_id, keysca, iscal
...
all field_get_id(scalar1, f_id)
call field_get_key_id("scalar_id", keysca)
call field_get_key_int(f_id, keysca, iscal)
When I put this, the message showed

Code: Select all

/home/tsubasa/Desktop/test/course/RESU/20211210-1300/src/cs_meg_initialization.c: In function ‘cs_meg_initialization’:
/home/tsubasa/Desktop/test/course/RESU/20211210-1300/src/cs_meg_initialization.c:47:17: warning: unused variable ‘c_id’ [-Wunused-variable]
       cs_lnum_t c_id = zone->elt_ids[e_id];
                 ^~~~
/home/tsubasa/Desktop/test/course/RESU/20211210-1300/src/cs_user_boundary_conditions.f90:485:25:

 call field_get_id(scalar1, f_id)
                         1
Error: Symbol ‘scalar1’ at (1) has no IMPLICIT type
scalar1.png
For this, shoudl I write something like

Code: Select all

field_get_name(...)
or
just should I define the "scalar1" at the begining of code?

Sorry this may be easy quesition for you, but I am confused by difference of program manner between C, Fortran and pyrhon and so on. I am relatively new for most of languages.

Best regards,
Hamada
Yvan Fournier
Posts: 4070
Joined: Mon Feb 20, 2012 3:25 pm

Re: Imposing scalar at inlet by user subroutine

Post by Yvan Fournier »

Hello,

Since "scalar1" is the name of your field, in the Fortran or C code, you need to put it in quotes:

call field_get_id("scalar1", f_id)

Regards,

Yvan
Tsubasa
Posts: 175
Joined: Fri Nov 20, 2020 4:09 am

Re: Imposing scalar at inlet by user subroutine

Post by Tsubasa »

Hello Yvan,

Thank you for the advice.
I remenbered it.

I suceeded to impose time history of scalar value for species at inlet.
I put them here, and if you find some mistakes, please let me know.

Code: Select all

integer i
real i_scalar, i_direchlet_scalar
real, save :: i_narray(50)         ! should be change : total calculation step for scalar calculation
integer :: total_nstep = 50        ! should be change : total calculation step for scalar

!!Dirichlet B.C for scalar transport:
if(ntcabs ==101)then ! This number should be (restart time step + 1)
  open (17, file='/home/tsubasa/Desktop/test/course/SRC/inlet_scalar_time_history.csv', status='old')
  read (17, '()') 
  do i = 1, total_nstep
    read (17, *) i_scalar
    i_narray(i) = i_scalar  
  end do
  close (17)
else
endif
call field_get_id("scalar1", f_id)
call field_get_key_id("scalar_id", keysca)
call field_get_key_int(f_id, keysca, iscal)
   icodcl(ifac,isca)=1
   rcodcl(ifac,isca,1)=i_narray(ntcabs-100) !You should subtract the number of previous calculation

Best regards,
Tsubasa
Yvan Fournier
Posts: 4070
Joined: Mon Feb 20, 2012 3:25 pm

Re: Imposing scalar at inlet by user subroutine

Post by Yvan Fournier »

Hello,

You can improve you code a bit by replacing "100" with "ntpabs", which is the number of the last time step of the prvevious computation when you do a restart (and 0 for a non-restarted case). This way, it will still work even if you restart after a different number of time steps.

Best regards,

Yvan
Tsubasa
Posts: 175
Joined: Fri Nov 20, 2020 4:09 am

Re: Imposing scalar at inlet by user subroutine

Post by Tsubasa »

Hello Yvan,

Thank you for that.
That was what I wanted to know.

Best regards,
Tsubasa

Code: Select all

!!Dirichlet B.C for scalar transport:
if(ntcabs ==ntpabs+1)then
  open (17, file='/home/tsubasa/Desktop/test/course/SRC/inlet_scalar_time_history.csv', status='old')
  read (17, '()')
  do i = 1, total_nstep
    read (17, *) i_scalar
    i_narray(i) = i_scalar 
  end do
  close (17)
else
endif
call field_get_id("scalar1", f_id)
call field_get_key_id("scalar_id", keysca)
call field_get_key_int(f_id, keysca, iscal)
   icodcl(ifac,isca)=1
   rcodcl(ifac,isca,1)=i_narray(ntcabs-ntpabs)

enddo
Post Reply