diffusion_benchmark_EbO_Scene_penalty.f90 Source File

program main
use FEMDomainClass
use DiffusionEquationClass
use SceneClass
use FactoryClass
implicit none

type(DiffusionEq_) :: solver
type(FEMDomain_) :: cube(4)
real(real64),allocatable :: c(:)
real(real64),allocatable :: DiffusionCoeff(:)
real(real64),allocatable :: Reaction(:)
real(real64),allocatable :: FixValue(:)
real(real64),allocatable :: C_n(:),x(:)
integer(int32),allocatable :: FixBoundary(:)
integer(int32) :: itr,offset,max_itr
real(real64) :: Length,epsi,dt,penalty,d,AL_lambda, AL_tol, AL_alpha,bc_val
type(IO_) :: f
logical :: passed
type(Scene_) :: scene
type(Factory_) :: factory

type(COO_) :: coo
type(CRS_) :: crs

bc_val = -1.0d0
! diffusion_benchmark_EbO_Scene.f90

Length=1.00d0
epsi=0.010d0
max_itr = 1000
dt = 0.030d0 
!dt = 0.0010d0 ! for Forward-Euler
!penalty = 10.0d0**(-4) ! for Forward-Euler
penalty = 10.0d0**(2) 
d       = 1.0e-2

AL_lambda= 0.00d0
AL_tol   = dble(5.0e-9)
AL_alpha =  3.00d0

! penaltyを上げるとgapが閉じずに発散する
cube = factory%cube([30,5,5],n=4)
call scene%add(cube)

call scene%resize(object=[1,2,3,4], &
    x=(Length+3*epsi)/4.0d0 ,y=Length/10.0d0,z=Length/10.0d0) 
    

call scene%move(object=[1], x=0.0d0)
call scene%move(object=[2], x=cube(1)%xmax()-epsi)
call scene%move(object=[3], x=cube(2)%xmax()-epsi)
call scene%move(object=[4], x=cube(3)%xmax()-epsi)

call scene%vtk("cube")
call scene%overset()

FixBoundary = scene%select_point_ID([1],x_max=scene%xmin([1]) ) &
    // scene%select_point_ID([4],x_min=scene%xmax([4]) )

FixValue = scene%full([1],values=bc_val,x_max=scene%xmin([1]) ) &
    // scene%full([4],values=bc_val,x_min=scene%xmax([4]) )

DiffusionCoeff = d*eyes(scene%nn() )
C_n = scene%full(func=sin_function)
call scene%vtk("init",scalar=C_n)


call f%open("init.csv","w")
x = scene%x()
do i_i=1,size(scene%x())
    write(f%fh,*) x(i_i), sin_function(x=[x(i_i)])  
enddo
call f%close()

call f%open("analytical.csv","w")
x = scene%x()
do i_i=1,size(scene%x())
    write(f%fh,*) x(i_i), analytical_solution(  &   
    x=[x(i_i), 0.0d0], &
    params=[DiffusionCoeff(1)]  )  
    
enddo
call f%close()


c = c_n

call solver%check_stability_condition(dt=dt,dx=Length/80,coefficient=dble(1.0e-4),passed=passed)
    if(.not.passed)stop

do itr=1,max_itr

    Reaction = zeros( scene%nn() )

    solver%solver%debug=True
    solver%solver%itrmax=100000
    
!    AL
!    c = solver%getDiffusionField(&
!        FEMDomains=cube, &
!        DiffusionCoeff=DiffusionCoeff, &
!        Reaction=Reaction, &
!        Penalty=penalty ,&
!        FixBoundary=FixBoundary, &
!        FixValue=FixValue   ,    &
!        C_n=C_n,&
!        AL_lambda=AL_lambda,&
!        AL_tol=AL_tol,&
!        AL_alpha=AL_alpha,&
!        dt=dt &
!        )

    c = solver%getDiffusionField(&
        FEMDomains=cube, &
        DiffusionCoeff=DiffusionCoeff, &
        Reaction=Reaction, &
        Penalty=penalty ,&
        FixBoundary=FixBoundary, &
        FixValue=FixValue   ,    &
        C_n=C_n,&
        dt=dt &
        )
    C_n = c


    call f%open("result.csv","w")
    call f%write(scene%x(),c(:) )
    call f%close()

    call f%cp("result.csv","result"+str(penalty)+zfill(itr,4)+".csv")


    call f%open("analytical.csv","w")
    x = scene%x()
    do i_i=1,size(scene%x())
        call f%write(x(i_i), analytical_solution(  &   
            x=[x(i_i), dt*itr], &
            params=[DiffusionCoeff(1)]  )  )
    enddo
    call f%close()
    
enddo




contains


function sin_function(x,params) result(ret)
    real(real64),intent(in) :: x(:)
    real(real64),optional,intent(in) :: params(:)
    real(real64) :: ret
    type(Math_) :: math

    ret = sin(math%pi*x(1))+bc_val

end function

function analytical_solution(x,params) result(ret)
    real(real64),intent(in) :: x(:)
    real(real64),optional,intent(in) :: params(:)
    real(real64) :: ret
    type(Math_) :: math

    ! x1 = x
    ! x2 = t
    ret = sin(math%PI*x(1) )&
                    *exp(-params(1)*math%PI*math%PI*x(2) )+bc_val

end function


end program