! ############################# module ReactorClass use MathClass implicit none integer(int32), parameter :: PF_MAX_REACTION_NUM = 20 integer(int32), parameter :: PF_MAX_INPUT_NUM = 5 integer(int32), parameter :: PF_MAX_OUTPUT_NUM = 5 type :: Reactor_ ! substances type(Real64Ptr_), allocatable :: substances(:) type(String_), allocatable :: substance_names(:) integer :: n_substances = 0 ! reactions integer :: n_reaction = 0 integer(int32):: input_list(PF_MAX_REACTION_NUM, PF_MAX_INPUT_NUM) = 0 integer(int32):: output_list(PF_MAX_REACTION_NUM, PF_MAX_OUTPUT_NUM) = 0 real(real64) :: reaction_order(PF_MAX_REACTION_NUM, PF_MAX_INPUT_NUM + PF_MAX_OUTPUT_NUM) = 0 real(real64) :: mol_rate(PF_MAX_REACTION_NUM, PF_MAX_INPUT_NUM + PF_MAX_OUTPUT_NUM) = 0 real(real64) :: constants(PF_MAX_REACTION_NUM) = 0 contains procedure :: Init => InitReactor procedure :: put => putReactor procedure :: define => defineReactor procedure :: SearchSubstanceID => SearchSubstanceIDReactor procedure :: run => runReactor end type contains ! --------------------------- subroutine InitReactor(this, max_num_substances) class(Reactor_), intent(inout) :: this integer(int32), optional, intent(in) :: max_num_substances integer(int32) :: n if (present(max_num_substances)) then n = max_num_substances else n = 100 end if this%n_substances = 0 if (allocated(this%substances)) then deallocate (this%substances) end if if (allocated(this%substance_names)) then deallocate (this%substance_names) end if allocate (this%substances(n)) allocate (this%substance_names(n)) end subroutine InitReactor ! --------------------------- ! --------------------------- subroutine putReactor(this, substance_name, substance) class(Reactor_), intent(inout) :: this character(*), intent(in) :: substance_name real(real64), target, intent(in) :: substance integer(int32) :: i, num num = this%n_substances + 1 if (num > size(this%substances)) then print *, "ERROR :: putReactor num>size(this%substances ) " else this%substances(num)%ptr => substance this%substance_names(num) = substance_name end if this%n_substances = this%n_substances + 1 end subroutine putReactor ! --------------------------- pure function SearchSubstanceIDReactor(this, name) result(ret) class(Reactor_), intent(in) :: this character(*), intent(in) :: name integer(int32) :: i, ret ret = 0 do i = 1, this%n_substances if (this%substance_names(i)%all == name) then ret = i return end if end do end function ! --------------------------- subroutine defineReactor(this, input_list, output_list, constant, mol_rate, reaction_order) class(Reactor_), intent(inout) :: this type(String_), intent(in) :: input_list(:) type(String_), intent(in) :: output_list(:) real(real64), intent(in) :: constant ! reacton constant k real(real64), optional, intent(in) :: reaction_order(:), mol_rate(:) real(real64), allocatable :: mrate(:), rrate(:) integer(int32) :: i, j, n this%n_reaction = this%n_reaction + 1 if (present(mol_rate)) then mrate = mol_rate else n = size(input_list) + size(output_list) allocate (mrate(n)) mrate(:) = 1.0d0 end if if (present(reaction_order)) then rrate = reaction_order else n = size(input_list) + size(output_list) allocate (rrate(n)) rrate(:) = 1.0d0 end if if (size(input_list) > PF_MAX_INPUT_NUM .or. size(output_list) > PF_MAX_OUTPUT_NUM) then print *, "ERROR :: defineReactor size(mrate) > PF_MAX_INPUT_NUM+PF_MAX_OUTPUT_NUM" stop end if ! register reaction ! INPUT do i = 1, size(input_list) ! search substance id this%input_list(this%n_reaction, i) = this%SearchSubstanceID(input_list(i)%all) end do do i = 1, size(output_list) ! search substance id this%output_list(this%n_reaction, i) = this%SearchSubstanceID(output_list(i)%all) end do this%reaction_order(this%n_reaction, 1:size(mrate)) = mrate(1:size(mrate)) this%mol_rate(this%n_reaction, 1:size(mrate)) = mrate(1:size(mrate)) this%constants(this%n_reaction) = constant end subroutine ! --------------------------- subroutine runReactor(this, dt) class(Reactor_), intent(inout) :: this real(real64), intent(in) :: dt integer(int32) :: i, j, num_input, num_output real(real64) :: A_B_ratio, xi, xi_n, coef, A0, B0, nA, nB real(real64) :: k1, k2, k3, k4 ! run simulation if (this%n_reaction == 0) then print *, "ERROR :: runReactor this%n_reaction==0" stop elseif (this%n_reaction == 1) then ! single reaction equation num_input = 0 do i = 1, size(this%input_list, 2) if (this%input_list(1, i) /= 0) then num_input = num_input + 1 end if end do num_output = 0 do i = 1, size(this%output_list, 2) if (this%output_list(1, i) /= 0) then num_output = num_output + 1 end if end do if (num_input == 1 .and. num_output == 1) then ! x*A -> y*B reaction equation print *, "[ok] Reactor Started >> " A_B_ratio = this%reaction_order(1, 2)/this%reaction_order(1, 1) this%substances(this%input_list(1, 1))%ptr = this%substances(this%input_list(1, 1))%ptr*exp(this%constants(1)*dt) this%substances(this%output_list(1, 1))%ptr = this%substances(this%output_list(1, 1))%ptr*exp(A_B_ratio*this%constants(1)*dt) elseif (num_input == 2 .and. num_output == 1) then print *, "[ok] Reactor Started >> " ! x*A + y*B -> z*C reaction equation ! A + y/x*B -> z/x*C reaction equation ! by Runge-Kutta mehod ! Reacted A = A0 - xi ! d xi/dt = k*(A0-xi)**(nA) * (B0 - xi)**(nB) ! 反応比率も必要！！ ! 反応速度式の次数は反応比率によらないが， ! 反応量は反応比率によるため． coef = this%constants(1) A0 = this%substances(this%input_list(1, 1))%ptr B0 = this%substances(this%input_list(1, 2))%ptr nA = this%reaction_order(1, 1) nB = this%reaction_order(1, 2) xi_n = 0.0d0 k1 = coef*((A0 - xi_n)**(nA))*((B0 - xi_n)**(nB)) k2 = coef*((A0 - (xi_n + dt*0.50d0*k1))**(nA)) & *((B0 - (xi_n + dt*0.50d0*k1))**(nB)) k3 = coef*((A0 - (xi_n + dt*0.50d0*k2))**(nA)) & *((B0 - (xi_n + dt*0.50d0*k2))**(nB)) k4 = coef*((A0 - (xi_n + dt*k3))**(nA)) & *((B0 - (xi_n + dt*k3))**(nB)) xi = xi_n + 1.0d0/6.0d0*dt*(k1 + k2 + k3 + k4) this%substances(this%input_list(1, 1))%ptr = A0 - xi this%substances(this%input_list(1, 2))%ptr = & B0 - dble(this%mol_rate(1, 2)/this%mol_rate(1, 1))*xi this%substances(this%output_list(1, 1))%ptr = & this%substances(this%output_list(1, 1))%ptr & + dble(this%mol_rate(1, 3)/this%mol_rate(1, 1))*xi else print *, "ERROR :: runReactor num_input/=1 .or. num_output/=1, now being implemented." stop end if else print *, "ERROR :: runReactor this%n_reaction>=2, now being implemented." stop end if end subroutine end module ReactorClass ! #############################