program main use SoybeanClass use ArrayClass use IOClass implicit none type(Soybean_) :: soy type(Environment_) :: env type(FEMDomain_),allocatable :: domains(:) integer(int32) :: hour real(real64),allocatable :: photosynthesis(:),carbon_concentration(:),& respiration(:),Photosynthate(:),internal_pressure(:),struct_vs_nonstruct(:),& void_ratio(:),DiffusionCoeff(:),volume(:),RHS(:),FixValue(:) integer(int32),allocatable :: target_element_list(:),FixBoundary(:) real(real64) :: dt type(CRS_) :: crs type(IO_) :: f ! initialize environment ! senario.csv: table data and time series ! senario.json: meta data ! all data should be contained. call env%init(json="Tutorial/obj/env.json") ! initialize soybean ! morphological data !soy%stem_division(1:3) = [2, 2, 20] !soy%peti_division(1:3) = [2, 2, 20] call soy%init(config="Tutorial/obj/soy.json") call soy%vtk("day"+zfill(0,4) ) !call soy%setPropertiesCarbonDiffusionCoefficient(default_value=0.00010d0) ! grow soybean dt = 1.d0 ! 1 sec do hour = 1, 100 ! (1-1) compute photosynthesis ratio (g/s) photosynthesis = soy%getPhotoSynthesis(env=env,dt=dt) ! unit: micro-gram call soy%vtk("photosyn_"+zfill(hour,4),scalar_field=photosynthesis,single_file=.true. ) volume = soy%getVolumePerElement() ! unit: micro-gram call soy%vtk("vol_"+zfill(hour,4),scalar_field=volume,single_file=.true. ) call soy%vtk("photosyn_vol_"+zfill(hour,4),scalar_field=photosynthesis/volume,single_file=.true. ) ! (1-2) search carbon requirement(sink capacity) (g/s) ! Fix carbo concentration (micro-g/m^3) at sink carbon_concentration = soy%getCarbon_concentration(env=env,& FixBoundary=FixBoundary,FixValue=FixValue) ! micro-gram/m^3/s call soy%vtk("carb_conc_"+zfill(hour,4),scalar_field=carbon_concentration,single_file=.true. ) ! get respiration (micro-g/m^3) respiration = soy%getRespiration(env=env) ! micro-gram/m^3/s call soy%vtk("carb_req_"+zfill(hour,4),scalar_field=respiration,single_file=.true. ) DiffusionCoeff = soy%getDiffusionCoefficient() call soy%vtk("diff_coef_"+zfill(hour,4),scalar_field=DiffusionCoeff,single_file=.true. ) ! (2) reaction-diffusion (g/s for each point) ! dC/dt = - D d^2/dx^2 c + P - R ! ソースの拡散係数(μg/m^3/m/s) DiffusionCoeff = dble(1.0e-15) ! Photosynthate = soy%getCarbonFlow( & dt = dt,& DiffusionCoeff = DiffusionCoeff,& photosynthesis=photosynthesis/volume,& ! reaction term (P) respiration = respiration,& ! reaction term (R) FixBoundary=[(i_i,i_i=1,10)],& FixValue=ones(10),& Photosynthate_n=soy%Photosynthate_n,& penalty = maxval(DiffusionCoeff)*(1.0),& RHS = RHS, & Matrix = CRS, & tol = dble(1.0e-1), & debug=.true. & ) ! initial condition soy%Photosynthate_n = Photosynthate call soy%vtk("biomass_"+zfill(hour,4),scalar_field=Photosynthate,single_file=.true. ) print *, size(RHS),maxval(RHS),minval(RHS) print *, soy%nn(),soy%ne() call soy%vtk("RHS_"+zfill(hour,4),scalar_field=RHS,single_file=.true. ) ! (3-1) decide ratio of structural carbon : non-structural carbon !struct_vs_nonstruct = soy%struct_vs_nonstruct( & ! Photosynthate=Photosynthate) ! (3-2) increment of biomass to internal pressure !internal_pressure = soy%internal_pressure( & ! Photosynthatet=Photosynthate,& ! struct_vs_nonstruct=struct_vs_nonstruct) ! (3-3) increment of biomass to void ratio !void_ratio = soy%internal_pressure( & ! Photosynthate=Photosynthate,& ! struct_vs_nonstruct=struct_vs_nonstruct) ! (4) increment of Young modulus and yield parameter to void ratio !call soy%update_Texture( & ! internal_pressure=internal_pressure,& ! struct_vs_nonstruct=struct_vs_nonstruct) ! (5) pressure-induced deformation !displacement = soy%static_deformation( & ! ground_level=0.0d0,& ! internal_pressure=internal_pressure) ! (6) deform soybean !call soy%deform(displacement=displacement) ! (7) remesh (create new object in apical) !call soy%remesh() ! (8) export morphology !call soy%vtk("day"+zfill(days,4) ) enddo end program