diff --git a/kmos/cli.py b/kmos/cli.py index 3badf2de..58f9bd9d 100644 --- a/kmos/cli.py +++ b/kmos/cli.py @@ -80,11 +80,16 @@ -s/--source-only Export source only and don't build binary - -b/--backend (local_smart|lat_int) + -b/--backend (local_smart|lat_int|otf) Choose backend. Default is "local_smart". lat_int is EXPERIMENTAL and not made for production, yet. + -w/--wasm + Export and compile for WebAssembly using flang-wasm. + This will generate .wasm and .js files instead of native binary. + Works with all backends (otf, local_smart, lat_int). + -d/--debug Turn on assertion statements in F90 code. (Only active in compile step) @@ -204,6 +209,15 @@ def get_options(args=None, get_parser=False): default=False, ) + parser.add_option( + "-w", + "--wasm", + dest="wasm", + action="store_true", + default=False, + help="Export and compile for WebAssembly using flang-wasm", + ) + # Detect available Fortran compiler # Note: numpy.distutils is deprecated and removed in Python >= 3.12 # Using direct detection instead @@ -327,12 +341,19 @@ def main(args=None): elif args[0] == "export": import kmos.types import kmos.io - from kmos.utils import build + from kmos.utils import build, build_wasm if len(args) < 2: parser.error("XML file and export path expected.") if len(args) < 3: - out_dir = "%s_%s" % (os.path.splitext(args[1])[0], options.backend) + # Include backend in directory name even for WASM builds + backend_suffix = "_%s" % options.backend + wasm_suffix = "_wasm" if options.wasm else "" + out_dir = "%s%s%s" % ( + os.path.splitext(args[1])[0], + backend_suffix, + wasm_suffix, + ) logger.info("No export path provided. Exporting to %s" % out_dir) args.append(out_dir) @@ -347,28 +368,120 @@ def main(args=None): kmos.io.export_source(project, export_dir, options=options) - if ( - (os.name == "posix" and os.uname()[0] in ["Linux", "Darwin"]) - or os.name == "nt" - ) and not options.source_only: + if not options.source_only: os.chdir(export_dir) - build(options) - for out in glob("kmc_*"): - if os.path.exists("../%s" % out): - if options.overwrite: - overwrite = "y" + + if options.wasm: + # Build for WebAssembly with enhanced output + import time + + logger.info("=" * 64) + logger.info("Building for WebAssembly (WASM)") + logger.info("=" * 64) + logger.info(f"Backend: {options.backend}") + logger.info(f"Export directory: {export_dir}") + logger.info("") + logger.info("This will:") + logger.info(" 1. Apply WASM-specific code modifications") + logger.info(" 2. Generate C bindings for browser interface") + logger.info(" 3. Compile Fortran to WASM using Docker") + logger.info("") + logger.info("Docker image: ghcr.io/r-wasm/flang-wasm:main") + logger.info("=" * 64) + logger.info("") + + try: + start_time = time.time() + result = build_wasm(options) + elapsed = time.time() - start_time + + # Get file sizes + js_size_kb = os.path.getsize("kmc_model.js") / 1024 + wasm_size_mb = os.path.getsize("kmc_model.wasm") / 1024 / 1024 + + logger.info("") + logger.info("=" * 64) + logger.info("WASM Build Successful!") + logger.info("=" * 64) + logger.info(f"Completed in {elapsed:.1f} seconds") + logger.info(f"JavaScript: kmc_model.js ({js_size_kb:.1f} KB)") + logger.info(f"WebAssembly: kmc_model.wasm ({wasm_size_mb:.2f} MB)") + logger.info("") + logger.info("Next steps:") + logger.info(" 1. Test the build in browser") + logger.info( + " 2. See kmos-web/test_kmc_wasm.html for example usage" + ) + logger.info("=" * 64) + logger.info("") + + # Move WASM files to parent directory + for out in glob("kmc_model.*"): + if out.endswith((".js", ".wasm")): + target = "../%s" % out + if os.path.exists(target): + if options.overwrite: + overwrite = "y" + else: + overwrite = input( + ("Should I overwrite existing %s ?[y/N] ") + % out + ).lower() + if overwrite.startswith("y"): + logger.info("Overwriting {out}".format(**locals())) + os.remove(target) + shutil.move(out, "..") + else: + logger.info("Skipping {out}".format(**locals())) + else: + shutil.move(out, "..") + + except RuntimeError as e: + logger.error("") + logger.error("=" * 64) + logger.error("WASM Build Failed") + logger.error("=" * 64) + logger.error(str(e)) + logger.error("") + logger.error("Troubleshooting:") + logger.error(" - Ensure Docker is running") + logger.error( + " - Try pulling the image: docker pull ghcr.io/r-wasm/flang-wasm:main" + ) + logger.error("") + logger.error("Alternative: Try source-only export:") + logger.error(f" kmos export {xml_file} -b {options.backend} -s") + logger.error("=" * 64) + raise + except IOError as e: + logger.error("") + logger.error("=" * 64) + logger.error("WASM Build Failed") + logger.error("=" * 64) + logger.error(str(e)) + logger.error("=" * 64) + raise + elif ( + os.name == "posix" and os.uname()[0] in ["Linux", "Darwin"] + ) or os.name == "nt": + # Build native binary + build(options) + for out in glob("kmc_*"): + if os.path.exists("../%s" % out): + if options.overwrite: + overwrite = "y" + else: + overwrite = input( + ("Should I overwrite existing %s ?[y/N] ") % out + ).lower() + if overwrite.startswith("y"): + logger.info("Overwriting {out}".format(**locals())) + os.remove("../%s" % out) + shutil.move(out, "..") + else: + logger.info("Skipping {out}".format(**locals())) else: - overwrite = input( - ("Should I overwrite existing %s ?[y/N] ") % out - ).lower() - if overwrite.startswith("y"): - logger.info("Overwriting {out}".format(**locals())) - os.remove("../%s" % out) shutil.move(out, "..") - else: - logger.info("Skipping {out}".format(**locals())) - else: - shutil.move(out, "..") elif args[0] == "settings-export": import kmos.io diff --git a/kmos/fortran_src/base.mpy b/kmos/fortran_src/base.mpy index 3ba57303..377d72d8 100644 --- a/kmos/fortran_src/base.mpy +++ b/kmos/fortran_src/base.mpy @@ -18,7 +18,7 @@ #@ ! along with kmos; if not, write to the Free Software #@ ! Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 #@ ! USA -#@ +#@ #@ !****h* kmos/base #@ ! FUNCTION #@ ! The base kMC module, which implements the kMC method on a :math:`d = 1` @@ -36,12 +36,12 @@ #@ use kind_values #@ !------ No implicit definition of variables ! #@ implicit none -#@ +#@ #@ !------ All variables, function and subroutine are private by default #@ private -#@ -#@ -#@ +#@ +#@ +#@ #@ ! The following subroutines and functions are made public #@ public :: add_proc, & #@ allocate_system, & @@ -79,11 +79,11 @@ #@ update_accum_rate, & #@ update_integ_rate, & #@ update_clocks -#@ -#@ +#@ +#@ #@ ! Public constants #@ integer(kind=iint) :: null_species = -1 -#@ +#@ #@ !---- Allocatable, module wide, variables #@ integer(kind=iint), dimension(:,:,:), allocatable, public :: avail_sites #@ !****v* base/avail_sites @@ -156,7 +156,7 @@ #@ ! Stores the total number of times each process has been executed #@ ! during one simulation. #@ !****** -#@ +#@ #@ real(kind=rdouble) :: kmc_time #@ !****v* base/kmc_time #@ ! FUNCTION @@ -182,9 +182,9 @@ #@ ! FUNCTION #@ ! The time increment of the current kMC step. #@ !****** -#@ -#@ -#@ +#@ +#@ +#@ #@ !--- Local copies of variables #@ integer(kind=iint) :: nr_of_proc #@ !****v* base/nr_of_proc @@ -203,11 +203,11 @@ #@ ! This name should not contain any characters that you don't want to #@ ! have in a filename either, i.e. only [A-Za-z0-9\_-]. #@ !****** -#@ +#@ #@ !**************** #@ contains #@ !**************** -#@ +#@ #@ subroutine del_proc(proc, site) #@ !****f* base/del_proc #@ ! FUNCTION @@ -229,27 +229,27 @@ #@ ! * ``site`` positive integer that encodes the site to be manipulated #@ !****** #@ integer(kind=iint), intent(in) :: proc, site -#@ +#@ #@ integer(kind=iint) :: memory_address -#@ +#@ #@ ! Make sure proc_nr is in the right range #@ ASSERT(proc.gt.0,"add_proc: proc has to be positive") #@ ASSERT(proc.le.nr_of_proc,"add_proc: proc has to be less or equal nr_of_proc.") #@ ! Make sure site is in the right range #@ ASSERT(site.gt.0,"add_proc: site has to be positive") #@ ASSERT(site.le.volume,"base/add_proc: site needs to be in volume") -#@ +#@ #@ ! assert consistency #@ ASSERT(avail_sites(proc, site, 2) .ne. 0 , "Error: tried to take ability from site that is not there!") -#@ +#@ #@ memory_address = avail_sites(proc, site, 2) #@ if(memory_address .lt. nr_of_sites(proc))then #@ ! check if we are deleting the last field -#@ +#@ #@ ! move last field to deleted field #@ avail_sites(proc, memory_address, 1) = avail_sites(proc, nr_of_sites(proc), 1) #@ avail_sites(proc, nr_of_sites(proc), 1) = 0 -#@ +#@ #@ ! change address of moved field #@ avail_sites(proc, avail_sites(proc, memory_address, 1), 2) = memory_address #@ else ! simply deleted last field @@ -257,14 +257,14 @@ #@ endif #@ ! delete address of deleted field #@ avail_sites(proc, site, 2) = 0 -#@ -#@ -#@ +#@ +#@ +#@ #@ ! decrement nr_of_sites(proc) #@ nr_of_sites(proc) = nr_of_sites(proc) - 1 #@ end subroutine del_proc -#@ -#@ +#@ +#@ #@ subroutine add_proc(proc, site) #@ !****f* base/add_proc #@ ! FUNCTION @@ -278,29 +278,29 @@ #@ ! * ``site`` positive integer number that represents the site to be manipulated #@ !****** #@ integer(kind=iint), intent(in) :: proc, site -#@ +#@ #@ ! Make sure proc_nr is in the right range #@ ASSERT(proc.gt.0,"base/add_proc: proc has to be positive") #@ ASSERT(proc.le.nr_of_proc,"base/add_proc: proc has to be less or equal nr_of_proc.") -#@ +#@ #@ ! Make sure site is in the right range #@ ASSERT(site.gt.0,"base/add_proc: site has to be positive") #@ ASSERT(site.le.volume,"base/add_proc: site needs to be in volume") -#@ +#@ #@ ! assert consistency #@ ASSERT(avail_sites(proc, site, 2) == 0, "base/add_proc Error: tried to add ability that is already there") -#@ +#@ #@ ! increment nr_of_sites(proc) #@ nr_of_sites(proc) = nr_of_sites(proc) + 1 -#@ +#@ #@ ! store site in nr_of_sites(proc)th slot #@ avail_sites(proc, nr_of_sites(proc), 1) = site -#@ +#@ #@ ! let address of added site point to nr_of_sites(proc)th slot #@ avail_sites(proc, site, 2) = nr_of_sites(proc) -#@ +#@ #@ end subroutine add_proc -#@ +#@ #@ pure function can_do(proc, site) #@ !****f* base/can_do #@ ! FUNCTION @@ -315,12 +315,12 @@ #@ !---------------I/O variables--------------- #@ logical :: can_do #@ integer(kind=iint), intent(in) :: proc, site -#@ +#@ #@ can_do = avail_sites(proc,site,2).ne.0 -#@ +#@ #@ end function can_do -#@ -#@ +#@ +#@ #@ subroutine reset_site(site, old_species) #@ !****f* base/reset_site #@ ! FUNCTION @@ -338,9 +338,9 @@ #@ integer(kind=iint), intent(in) :: site, old_species #@ !---------------internal variables--------------- #@ integer(kind=iint) :: proc, species -#@ +#@ #@ species = get_species(site) -#@ +#@ #@ ! Reset species if stated correctly #@ if(old_species.eq.species)then #@ call replace_species(site, species, null_species) @@ -349,17 +349,17 @@ #@ print *,'Expected',old_species,'but found',species,'on',site #@ stop #@ endif -#@ +#@ #@ ! Strip all available capabilities from this site #@ do proc = 1, nr_of_proc #@ if(can_do(proc, site))then #@ call del_proc(proc, site) #@ endif #@ enddo -#@ +#@ #@ end subroutine reset_site -#@ -#@ +#@ +#@ #@ subroutine reload_system(input_system_name, reloaded) #@ !****f* base/reload_system #@ ! FUNCTION @@ -375,22 +375,22 @@ #@ !---------------I/O variables--------------- #@ character(len=200), intent(in) :: input_system_name #@ logical, intent(out) :: reloaded -#@ +#@ #@ character(len=210) :: filename #@ character(len=20) :: label #@ character(len=10**6) :: buffer #@ integer :: io_state , line, pos, subindex, io -#@ +#@ #@ integer, parameter :: filehandler = 15 #@ logical :: file_exists -#@ -#@ +#@ +#@ #@ ! store system name in module variable #@ system_name = input_system_name #@ ! initialize input/output flag #@ io_state = 0 #@ line = 0 -#@ +#@ #@ write(filename,'(a,a)')TRIM(ADJUSTL(system_name)),'.reload' #@ inquire(file=trim(adjustl(filename)),exist=file_exists) #@ if(.not.file_exists)then @@ -400,7 +400,7 @@ #@ else #@ ! Open file #@ open(filehandler, file = filename) -#@ +#@ #@ ! First parse loop: parse scalar values and system size #@ ! to allocate appropriate arrays for 2nd loop #@ do while(io_state == 0) @@ -409,10 +409,10 @@ #@ if(io_state == 0) then #@ ! advance line number #@ line = line + 1 -#@ +#@ #@ ! Shuffle all non-space characters to the left. #@ buffer = adjustl(buffer) -#@ +#@ #@ ! Ignore comment lines #@ if(buffer(1:1) == '#')then #@ cycle @@ -423,7 +423,7 @@ #@ label = buffer(1:pos) #@ ! Everythings else remains in the buffer #@ buffer = buffer(pos+1:) -#@ +#@ #@ ! In the first go we are only interested in the variables that #@ ! determine the system's size (in memory) #@ select case(label) @@ -434,15 +434,15 @@ #@ endselect #@ endif #@ enddo -#@ +#@ #@ if(io_state>0)then #@ print *,"Some read error occured in the first reload loop, investigate!" #@ stop #@ endif -#@ -#@ +#@ +#@ #@ call allocate_system(nr_of_proc, volume, system_name) -#@ +#@ #@ ! Second loop: parse the "meat" of data #@ rewind(filehandler) #@ io_state = 0 @@ -451,9 +451,9 @@ #@ read(filehandler, '(a)', iostat=io_state) buffer #@ if(io_state == 0) then #@ line = line + 1 -#@ +#@ #@ buffer = adjustl(buffer) -#@ +#@ #@ ! Ignore comment lines #@ if(buffer(1:1) == '#')then #@ cycle @@ -486,10 +486,10 @@ #@ label = buffer(1:pos) #@ buffer = buffer(pos+1:) #@ buffer = adjustl(buffer) -#@ +#@ #@ read(label, * ,iostat = io_state) subindex #@ read(buffer , *, iostat = io) avail_sites(subindex,:,1) -#@ +#@ #@ case('avail_sites_back') #@ buffer = adjustl(buffer) #@ pos = scan(buffer, ' ') @@ -497,7 +497,7 @@ #@ buffer = buffer(pos+1:) #@ read(label, *, iostat = io) subindex #@ read(buffer , *, iostat = io) avail_sites(subindex,:,2) -#@ +#@ #@ endselect #@ endif #@ enddo @@ -505,15 +505,15 @@ #@ print *,"Some read error occured in the second reload loop, investigate!" #@ stop #@ endif -#@ +#@ #@ close(filehandler) -#@ +#@ #@ reloaded = .true. #@ endif -#@ +#@ #@ end subroutine reload_system -#@ -#@ +#@ +#@ #@ subroutine save_system() #@ !****f* base/save_system #@ ! FUNCTION @@ -537,9 +537,9 @@ #@ integer, parameter :: filehandler = 15 #@ character(len=210) :: filename #@ integer(kind=iint) :: i, io_state -#@ +#@ #@ character(len=10) :: dummy_string -#@ +#@ #@ write(filename,'(2a)',iostat=io_state) trim(adjustl(system_name)),'.reload' #@ open(filehandler, file=filename) #@ ! Write scalar fields @@ -550,18 +550,18 @@ #@ write(filehandler,'(a,i22)')' kmc_step ',kmc_step #@ write(filehandler,*)'nr_of_proc ',nr_of_proc #@ write(filehandler,*)'volume ',volume -#@ +#@ #@ ! Write array fields #@ write(dummy_string,'(i9)') nr_of_proc -#@ +#@ #@ write(filehandler,'(a)')"#Vector variables" #@ write(filehandler,'(a,'//trim(adjustl(dummy_string))//'i21)')'procstat ',procstat #@ write(filehandler,'(a,'//trim(adjustl(dummy_string))//'i9)')'nr_of_sites ',nr_of_sites #@ write(filehandler,'(a,'//trim(adjustl(dummy_string))//'es14.7)')'rates ',rates -#@ +#@ #@ write(dummy_string,'(i9)') volume #@ write(filehandler,'(a,'//trim(adjustl(dummy_string))//'i9)')'lattice ',lattice -#@ +#@ #@ ! Avail_sites need one more field than 'volume' because first one describes the row #@ write(dummy_string,'(i9)') volume+1 #@ do i = 1, nr_of_proc @@ -570,14 +570,14 @@ #@ do i = 1, nr_of_proc #@ write(filehandler,'(a,'//trim(adjustl(dummy_string))//'i9)')'avail_sites_back ',i,avail_sites(i,:,2) #@ enddo -#@ -#@ -#@ +#@ +#@ +#@ #@ close(filehandler) -#@ +#@ #@ end subroutine save_system -#@ -#@ +#@ +#@ #@ subroutine set_rate_const(proc_nr, rate) #@ !****f* base/set_rate_const #@ ! FUNCTION @@ -590,16 +590,16 @@ #@ !****** #@ integer(kind=iint), intent(in) :: proc_nr #@ real(kind=rdouble), intent(in) :: rate -#@ +#@ #@ ! Make sure proc_nr is in the right range #@ ASSERT(proc_nr.gt.0,"base/set_rate_const: proc_nr has to be positive") #@ ! * the field within the process, but the meaning differs as explained #@ ! under 'switch' #@ ASSERT(proc_nr.le.nr_of_proc,"base/set_rate_const: proc_nr less or equal nr_of_proc.") #@ rates(proc_nr) = rate -#@ +#@ #@ end subroutine set_rate_const -#@ +#@ #@ subroutine update_accum_rate() #@ !****f* base/update_accum_rate #@ ! FUNCTION @@ -609,19 +609,19 @@ #@ ! #@ ! ``none`` #@ !****** -#@ +#@ #@ integer(kind=iint) :: i -#@ +#@ #@ accum_rates(1)=nr_of_sites(1)*rates(1) #@ do i = 2, nr_of_proc #@ accum_rates(i)=accum_rates(i-1)+nr_of_sites(i)*rates(i) #@ enddo -#@ +#@ #@ ASSERT(accum_rates(nr_of_proc).gt.0.,"base/update_accum_rate found & #@ accum_rates(nr_of_proc)=0, so no process is available at all") -#@ +#@ #@ end subroutine update_accum_rate -#@ +#@ #@ !------ S. Matera 09/18/2012------ #@ subroutine update_integ_rate() #@ !****f* base/update_integ_rate @@ -632,19 +632,19 @@ #@ ! #@ ! ``none`` #@ !****** -#@ +#@ #@ integer(kind=iint) :: i -#@ -#@ +#@ +#@ #@ do i = 1, nr_of_proc #@ integ_rates(i)=integ_rates(i)+nr_of_sites(i)*rates(i)*kmc_time_step #@ enddo -#@ +#@ #@ ASSERT(accum_rates(nr_of_proc).gt.0.,"base/update_accum_rate: no process available") -#@ +#@ #@ end subroutine update_integ_rate #@ !------ S. Matera 09/18/2012------ -#@ +#@ #@ subroutine allocate_system(input_nr_of_proc, input_volume, input_system_name) #@ !****f* base/allocate_system #@ ! FUNCTION @@ -660,22 +660,22 @@ #@ character(len=200), intent(in) :: input_system_name #@ integer(kind=iint), intent(in) :: input_volume, input_nr_of_proc #@ logical :: system_allocated -#@ +#@ #@ system_allocated = .false. -#@ -#@ +#@ +#@ #@ ! Make sure we have at least one process #@ if(input_nr_of_proc.le.0)then #@ print *,"kmos/base/allocate_system: there needs to be at least one process in a kMC system" #@ stop #@ endif -#@ +#@ #@ ! Make sure we have at least one site #@ if(input_volume.le.0)then #@ print *,"kmos/base/allocate_system: there needs to be at least one site in the system" #@ stop #@ endif -#@ +#@ #@ ! Make sure we don't try to allocate twice #@ if(allocated(avail_sites))then #@ print *,"kmos/base/allocate_system: Tried to allocate avail_sites twice, please deallocate first" @@ -707,19 +707,19 @@ #@ print *,"kmos/base/allocate_system: Tried to allocate procstat twice, please deallocate first" #@ system_allocated = .true. #@ endif -#@ +#@ #@ if(.not. system_allocated)then #@ ! copy arguments to module variables #@ nr_of_proc = input_nr_of_proc #@ volume = input_volume #@ system_name = input_system_name -#@ +#@ #@ ! Set clocks and step counter to 0 #@ kmc_time = 0. #@ walltime = 0. #@ start_time = 0. #@ kmc_step = 0 -#@ +#@ #@ ! allocate data structures and initialize with 0 #@ allocate(avail_sites(nr_of_proc, volume, 2)) #@ avail_sites = 0 @@ -737,19 +737,19 @@ #@ !------ S. Matera 09/18/2012------ #@ allocate(procstat(nr_of_proc)) #@ procstat = 0 -#@ +#@ #@ endif -#@ -#@ +#@ +#@ #@ end subroutine allocate_system -#@ -#@ +#@ +#@ #@ subroutine is_allocated(result) #@ logical, intent(out) :: result #@ result = allocated(avail_sites) #@ end subroutine is_allocated -#@ -#@ +#@ +#@ #@ subroutine deallocate_system() #@ !****f* base/deallocate_system #@ ! FUNCTION @@ -797,10 +797,10 @@ #@ else #@ print *,"Warning: rates was not procstat, tried to deallocate." #@ endif -#@ +#@ #@ end subroutine deallocate_system -#@ -#@ +#@ +#@ #@ subroutine get_system_name(output_system_name) #@ !****f* base/get_system_name #@ ! FUNCTION @@ -815,8 +815,8 @@ #@ #@ output_system_name = system_name #@ end subroutine get_system_name -#@ -#@ +#@ +#@ #@ subroutine set_system_name(input_system_name) #@ !****f* base/set_system_name #@ ! FUNCTION @@ -828,12 +828,12 @@ #@ ! * ``system_name`` Readable string of type character(len=200). #@ !****** #@ character(len=200), intent(in) :: input_system_name -#@ +#@ #@ system_name = input_system_name -#@ +#@ #@ end subroutine set_system_name -#@ -#@ +#@ +#@ #@ subroutine set_kmc_time(new_kmc_time) #@ !****f* base/set_kmc_time #@ ! FUNCTION @@ -845,12 +845,12 @@ #@ !****** #@ !---------------I/O variables--------------- #@ real(kind=rdouble), intent(in) :: new_kmc_time -#@ +#@ #@ kmc_time = new_kmc_time -#@ +#@ #@ end subroutine set_kmc_time -#@ -#@ +#@ +#@ #@ subroutine get_kmc_time(return_kmc_time) #@ !****f* base/get_kmc_time #@ ! FUNCTION @@ -862,12 +862,12 @@ #@ !****** #@ !---------------I/O variables--------------- #@ real(kind=rdouble), intent(out) :: return_kmc_time -#@ +#@ #@ return_kmc_time = kmc_time -#@ +#@ #@ end subroutine get_kmc_time -#@ -#@ +#@ +#@ #@ subroutine get_kmc_time_step(return_kmc_time_step) #@ !****f* base/get_kmc_time_step #@ ! FUNCTION @@ -879,12 +879,12 @@ #@ !****** #@ !---------------I/O variables--------------- #@ real(kind=rdouble), intent(out) :: return_kmc_time_step -#@ +#@ #@ return_kmc_time_step = kmc_time_step -#@ +#@ #@ end subroutine get_kmc_time_step -#@ -#@ +#@ +#@ #@ subroutine get_procstat(proc, return_procstat) #@ !****f* base/get_procstat #@ ! FUNCTION @@ -898,12 +898,12 @@ #@ !---------------I/O variables--------------- #@ integer(kind=iint),intent(in) :: proc #@ integer(kind=ilong),intent(out) :: return_procstat -#@ +#@ #@ return_procstat = procstat(proc) -#@ +#@ #@ end subroutine get_procstat -#@ -#@ +#@ +#@ #@ subroutine get_nrofsites(proc, return_nrofsites) #@ !****f* base/get_nrofsites #@ ! FUNCTION @@ -917,11 +917,11 @@ #@ !****** #@ integer(kind=iint), intent(in) :: proc #@ integer(kind=iint), intent(out) :: return_nrofsites -#@ +#@ #@ return_nrofsites = nr_of_sites(proc) #@ end subroutine get_nrofsites -#@ -#@ +#@ +#@ #@ subroutine get_avail_site(proc_nr, field, switch, return_avail_site) #@ !****f* base/get_avail_site #@ ! FUNCTION @@ -936,12 +936,12 @@ #@ !---------------I/O variables--------------- #@ integer(kind=iint), intent(in) :: proc_nr, field, switch #@ integer(kind=iint), intent(out) :: return_avail_site -#@ +#@ #@ return_avail_site = avail_sites(proc_nr, field, switch) -#@ +#@ #@ end subroutine get_avail_site -#@ -#@ +#@ +#@ #@ subroutine get_accum_rate(proc_nr, return_accum_rate) #@ !****f* base/get_accum_rate #@ ! FUNCTION @@ -955,15 +955,15 @@ #@ !---------------I/O variables--------------- #@ integer(kind=iint), intent(in), optional :: proc_nr #@ real(kind=rdouble), intent(out) :: return_accum_rate -#@ +#@ #@ if(.not. present(proc_nr) .or. proc_nr.eq.0) then #@ return_accum_rate=accum_rates(nr_of_proc) #@ else #@ return_accum_rate=accum_rates(proc_nr) #@ endif -#@ +#@ #@ end subroutine get_accum_rate -#@ +#@ #@ !------ S. Matera 09/18/2012------ #@ subroutine get_integ_rate(proc_nr, return_integ_rate) #@ !****f* base/get_integ_rate @@ -978,16 +978,16 @@ #@ !---------------I/O variables--------------- #@ integer(kind=iint), intent(in), optional :: proc_nr #@ real(kind=rdouble), intent(out) :: return_integ_rate -#@ +#@ #@ if(.not. present(proc_nr) .or. proc_nr.eq.0) then #@ return_integ_rate=integ_rates(nr_of_proc) #@ else #@ return_integ_rate=integ_rates(proc_nr) #@ endif -#@ +#@ #@ end subroutine get_integ_rate #@ !------ S. Matera 09/18/2012------ -#@ +#@ #@ subroutine get_rate(proc_nr, return_rate) #@ !****f* base/get_rate #@ ! FUNCTION @@ -1001,12 +1001,12 @@ #@ !---------------I/O variables--------------- #@ integer(kind=iint), intent(in) :: proc_nr #@ real(kind=rdouble), intent(out) :: return_rate -#@ +#@ #@ return_rate=rates(proc_nr) -#@ +#@ #@ end subroutine get_rate -#@ -#@ +#@ +#@ #@ subroutine increment_procstat(proc) #@ !****f* base/increment_procstat #@ ! FUNCTION @@ -1017,12 +1017,12 @@ #@ ! * ``proc`` integer representing the process to be increment. #@ !****** #@ integer(kind=iint),intent(in) :: proc -#@ +#@ #@ procstat(proc) = procstat(proc) + 1 -#@ +#@ #@ end subroutine increment_procstat -#@ -#@ +#@ +#@ #@ subroutine get_walltime(return_walltime) #@ !****f* base/get_walltime #@ ! FUNCTION @@ -1034,11 +1034,11 @@ #@ !****** #@ !---------------I/O variables--------------- #@ real(kind=rsingle), intent(out) :: return_walltime -#@ +#@ #@ return_walltime = walltime -#@ +#@ #@ end subroutine get_walltime -#@ +#@ #@ subroutine get_volume(return_volume) #@ !****f* base/get_kmc_volume #@ ! FUNCTION @@ -1050,11 +1050,11 @@ #@ !****** #@ !---------------I/O variables--------------- #@ integer(kind=iint), intent(out) :: return_volume -#@ +#@ #@ return_volume = volume -#@ +#@ #@ end subroutine get_volume -#@ +#@ #@ subroutine get_kmc_step(return_kmc_step) #@ !****f* base/get_kmc_step #@ ! FUNCTION @@ -1066,12 +1066,12 @@ #@ !****** #@ !---------------I/O variables--------------- #@ integer(kind=ilong), intent(out) :: return_kmc_step -#@ +#@ #@ return_kmc_step = kmc_step -#@ +#@ #@ end subroutine get_kmc_step -#@ -#@ +#@ +#@ #@ subroutine determine_procsite(ran_proc, ran_site, proc, site) #@ !****f* base/determine_procsite #@ ! FUNCTION @@ -1092,34 +1092,34 @@ #@ real(kind=rsingle), intent(in) :: ran_proc, ran_site #@ integer(kind=iint), intent(out) :: proc, site #@ !---------------internal variables--------------- -#@ -#@ +#@ +#@ #@ ASSERT(ran_proc.ge.0,"base/determine_procsite: ran_proc has to be positive") #@ ASSERT(ran_proc.le.1,"base/determine_procsite: ran_proc has to be less or equal 1") #@ ASSERT(ran_site.ge.0,"base/determine_procsite: ran_site has to be positive") #@ ASSERT(ran_site.le.1,"base/determine_procsite: ran_site has to be less or equal 1") -#@ +#@ #@ ! ran_proc <- [0,1] so we multiply with larger value in accum_rates #@ call interval_search_real(accum_rates, ran_proc*accum_rates(nr_of_proc), proc) -#@ -#@ +#@ +#@ #@ ! the result shall be between 1 and nrofsite(proc) so we have to add 1 the #@ ! scaled random number. But if the random number is closer to 1 than machine #@ ! precision, e.g. 0.999999999, we would get nrofsits(proc)+1 so we have to #@ ! cap it with min(...) #@ site = avail_sites(proc, & #@ min(nr_of_sites(proc),int(1+ran_site*(nr_of_sites(proc)))),1) -#@ -#@ +#@ +#@ #@ ASSERT(nr_of_sites(proc).gt.0,"base/determine_procsite: chosen process is invalid & #@ because it has no sites available.") #@ ASSERT(site.gt.0,"kmos/base/determine_procsite: tries to return invalid site") #@ ASSERT(site.le.volume,"base/determine_procsite: tries to return site larger than volume") -#@ -#@ +#@ +#@ #@ end subroutine determine_procsite -#@ -#@ +#@ +#@ #@ subroutine update_clocks(ran_time) #@ !****f* base/update_clocks #@ ! FUNCTION @@ -1131,25 +1131,25 @@ #@ !****** #@ real(kind=rsingle), intent(in) :: ran_time #@ real(kind=rsingle) :: runtime -#@ -#@ +#@ +#@ #@ ! Make sure ran_time is in the right interval #@ ASSERT(ran_time.ge.0.,"base/update_clocks: ran_time variable has to be positive.") #@ ASSERT(ran_time.le.1.,"base/update_clocks: ran_time variable has to be less than 1.") -#@ +#@ #@ kmc_time_step = -log(ran_time)/accum_rates(nr_of_proc) #@ ! Make sure the difference is not so small, that it is rounded off #@ ! ASSERT(kmc_time+kmc_time_step>kmc_time,"base/update_clocks: precision of kmc_time is not sufficient") -#@ +#@ #@ call CPU_TIME(runtime) -#@ +#@ #@ ! Make sure we are not dividing by zero #@ ASSERT(accum_rates(nr_of_proc).gt.0,"base/update_clocks: total rate was found to be zero") #@ kmc_time = kmc_time + kmc_time_step -#@ +#@ #@ ! Increment kMC steps #@ kmc_step = kmc_step + 1 -#@ +#@ #@ ! Walltime is the time of this simulation run plus the walltime #@ ! when the simulation was reloaded, so walltime represents the total #@ ! walltime across reloads. @@ -1159,8 +1159,8 @@ #@ !call update_integ_rate() #@ !------ S. Matera 09/18/2012------ #@ end subroutine update_clocks -#@ -#@ +#@ +#@ #@ pure function get_species(site) #@ !****f* base/get_species #@ ! FUNCTION @@ -1173,17 +1173,17 @@ #@ !---------------I/O variables--------------- #@ integer(kind=iint) :: get_species #@ integer(kind=iint), intent(in) :: site -#@ +#@ #@ !! DEBUG #@ !print *, site #@ !ASSERT(site.ge.1,"kmos/base/get_species was asked for a zero or negative site") #@ !ASSERT(site.le.volume,"kmos/base/get_species was asked for a site outside the lattice") -#@ +#@ #@ get_species = lattice(site) -#@ +#@ #@ end function get_species -#@ -#@ +#@ +#@ #@ subroutine replace_species(site, old_species, new_species) #@ !****f* base/replace_species #@ ! FUNCTION @@ -1198,9 +1198,9 @@ #@ ! * ``new_species`` integer representing the species to be placed #@ !****** #@ integer(kind=iint), intent(in) :: site, old_species, new_species -#@ +#@ #@ ASSERT(site.le.volume,"kmos/base/replace_species was asked for a site outside the lattice") -#@ +#@ #@ ! Double-check that we actually remove the atom that we think is there #@ if(old_species.ne.lattice(site))then #@ print '(a)', "kmos/base/replace_species Tried to remove species from sites which is not there!" @@ -1222,15 +1222,15 @@ #@ print '(a,i2,a,i2,a,i2,a,i10,a,i10,a)', & #@ "model.post_mortem(err_code=(",old_species,", ",new_species, ", ", lattice(site), ", ", site, ", ", kmc_step, "))" #@ print '(a)', "model.show()" -#@ -#@ +#@ +#@ #@ stop #@ endif -#@ +#@ #@ lattice(site) = new_species #@ end subroutine replace_species -#@ -#@ +#@ +#@ #@ subroutine interval_search_real(arr, value, return_field) #@ !****f* base/interval_search_real #@ ! FUNCTION @@ -1266,14 +1266,14 @@ #@ integer(kind=iint), intent(out) :: return_field #@ !---------------internal variables--------------- #@ integer(kind=iint) :: left, mid, right -#@ +#@ #@ left = 1 #@ right = size(arr) -#@ +#@ #@ binarysearch: do #@ ! Determine the middle between left and right by bit shifting #@ mid = ISHFT(right+left, -1) -#@ +#@ #@ ! if left and right overlap, we are done for now #@ if(left.ge.right)then #@ exit binarysearch @@ -1287,7 +1287,7 @@ #@ ! So now, we have found a candidate fields, which is great. #@ ! We just have to make sure it fullfills all other requirements. #@ enddo binarysearch -#@ +#@ #@ ! If the value turns out to be zero, we do a linear search #@ ! to the right until we find a non-zero entry #@ if(arr(mid).eq.0.)then @@ -1309,8 +1309,8 @@ #@ endif #@ enddo nonzerosearch #@ endif -#@ -#@ +#@ +#@ #@ ! If we are here, the value is non-zero which is great, but we also #@ ! have to make sure we return the left-most field if one or more fields #@ ! have the same value. @@ -1325,18 +1325,18 @@ #@ ASSERT(arr(mid).gt.arr(mid-1),"interval_search_real did not return a leftmost") #@ endif #@ enddo leftmostsearch -#@ -#@ +#@ +#@ #@ ASSERT(mid>0,"Returned index has to be at least 1") #@ ASSERT(mid<=size(arr),"Returned index can be at most size(arr)") #@ ASSERT(arr(mid).gt.0.,"Value of returned field has to be greater then 0") #@ ASSERT(value.le.arr(mid),"interval_search_real has an internal error") -#@ +#@ #@ return_field = mid -#@ -#@ +#@ +#@ #@ end subroutine interval_search_real -#@ +#@ #@ subroutine assertion_fail(a, r) #@ !****f* base/assertion_fail #@ ! FUNCTION @@ -1358,22 +1358,22 @@ #@ ' walltime: '//trim(wt)// & #@ ' kmc_time:'//trim(kt) #@ stop -#@ +#@ #@ end subroutine assertion_fail -#@ +#@ #@ subroutine set_null_species(input_null_species) #@ integer(kind=iint), intent(in) :: input_null_species -#@ +#@ #@ null_species = input_null_species -#@ +#@ #@ end subroutine set_null_species -#@ +#@ #@ subroutine get_null_species(output_null_species) #@ integer(kind=iint), intent(out) :: output_null_species -#@ +#@ #@ output_null_species = null_species -#@ +#@ #@ end subroutine get_null_species -#@ -#@ +#@ +#@ #@ end module base diff --git a/kmos/fortran_src/kind_values.f90 b/kmos/fortran_src/kind_values.f90 index bf6a50cb..f34ec976 100644 --- a/kmos/fortran_src/kind_values.f90 +++ b/kmos/fortran_src/kind_values.f90 @@ -7,10 +7,10 @@ module kind_values !****** implicit none -integer, parameter :: rsingle = SELECTED_REAL_KIND(p=15, r=200) -integer, parameter :: rdouble = SELECTED_REAL_KIND(p=15, r=200) -integer, parameter :: ibyte = SELECTED_INT_KIND(p=2) -integer, parameter :: ishort = SELECTED_INT_KIND(p=4) -integer, parameter :: iint = SELECTED_INT_KIND(p=9) -integer, parameter :: ilong = SELECTED_INT_KIND(p=18) +integer, parameter :: rsingle = SELECTED_REAL_KIND(15, 200) +integer, parameter :: rdouble = SELECTED_REAL_KIND(15, 200) +integer, parameter :: ibyte = SELECTED_INT_KIND(2) +integer, parameter :: ishort = SELECTED_INT_KIND(4) +integer, parameter :: iint = SELECTED_INT_KIND(9) +integer, parameter :: ilong = SELECTED_INT_KIND(18) end module kind_values diff --git a/kmos/fortran_src/lattice.mpy b/kmos/fortran_src/lattice.mpy index 3a233184..1612bb44 100644 --- a/kmos/fortran_src/lattice.mpy +++ b/kmos/fortran_src/lattice.mpy @@ -256,7 +256,7 @@ elif data.meta.model_dimension == 1: #@ end do #@ end do #@ -#@ do check_nr=1, product(system_size)*spuck +#@ do check_nr=1, system_size(1)*system_size(2)*system_size(3)*spuck #@ if(.not.check_nr.eq.calculate_lattice2nr(calculate_nr2lattice(check_nr)))then #@ print *, "ERROR in Mapping:", check_nr #@ print *, "was mapped on", calculate_nr2lattice(check_nr) @@ -267,14 +267,14 @@ elif data.meta.model_dimension == 1: #@ if data.meta.debug > 1: #@ print *, " LATTICE/ALLOCATE_SYSTEM/MAPPING_OK" -#@ allocate(nr2lattice(1:product(system_size)*spuck,4)) +#@ allocate(nr2lattice(1:system_size(1)*system_size(2)*system_size(3)*spuck,4)) #@ allocate(lattice2nr(-system_size(1):2*system_size(1)-1, & #@ -system_size(2):2*system_size(2)-1, & #@ -system_size(3):2*system_size(3)-1, & #@ 1:spuck)) if data.meta.debug > 1: #@ print *, " LATTICE/ALLOCATE_SYSTEM/ALLOCATED_LOOKUP" -#@ do check_nr=1, product(system_size)*spuck +#@ do check_nr=1, system_size(1)*system_size(2)*system_size(3)*spuck #@ nr2lattice(check_nr, :) = calculate_nr2lattice(check_nr) #@ end do #@ do k = -system_size(3), 2*system_size(3)-1 diff --git a/kmos/utils/__init__.py b/kmos/utils/__init__.py index 2b24bf3f..a0332b37 100644 --- a/kmos/utils/__init__.py +++ b/kmos/utils/__init__.py @@ -21,7 +21,6 @@ import logging import re import os -from time import time from io import StringIO from kmos.utils.ordered_dict import OrderedDict as OrderedDict @@ -374,11 +373,11 @@ def real_kind(args): infile = open(infile) outfile = open(outfile, "w") int_pattern = re.compile( - (r"(?P.*)selected_int_kind" "\((?P.*)\)(?P.*)"), + r"(?P.*)selected_int_kind\((?P.*)\)(?P.*)", flags=re.IGNORECASE, ) real_pattern = re.compile( - (r"(?P.*)selected_real_kind" "\((?P.*)\)(?P.*)"), + r"(?P.*)selected_real_kind\((?P.*)\)(?P.*)", flags=re.IGNORECASE, ) @@ -522,120 +521,877 @@ def build(options): sys.argv = true_argv -def T_grid(T_min, T_max, n): - from numpy import linspace, array +def _run_docker_command(cmd, src_dir, operation_name, timeout=300): + """Helper to run Docker commands with consistent error handling. - """Return a list of n temperatures between - T_min and T_max such that the grid of T^(-1) - is evenly spaced. + Args: + cmd: Shell command to run inside container + src_dir: Source directory to mount + operation_name: Human-readable name for the operation (e.g., "Compilation") + timeout: Timeout in seconds (default: 300) + + Raises: + RuntimeError: If command fails or times out """ + import subprocess - T_min1 = T_min ** (-1.0) - T_max1 = T_max ** (-1.0) + docker_cmd = [ + "docker", + "run", + "--rm", + "--platform", + "linux/amd64", + "-v", + f"{src_dir}:/src", + "ghcr.io/r-wasm/flang-wasm:main", + "sh", + "-c", + cmd, + ] - grid = list(linspace(T_max1, T_min1, n)) - grid.reverse() - grid = [x ** (-1.0) for x in grid] + try: + result = subprocess.run( + docker_cmd, capture_output=True, text=True, check=True, timeout=timeout + ) + if result.stdout: + logger.debug(result.stdout) + if result.stderr: + logger.debug(result.stderr) + return result + + except subprocess.TimeoutExpired: + raise RuntimeError( + f"{operation_name} timed out after {timeout} seconds. " + "Try reducing model complexity or increasing timeout." + ) + except subprocess.CalledProcessError as e: + logger.error(f"{operation_name} failed with exit code {e.returncode}") + if e.stdout: + logger.error(f"stdout: {e.stdout}") + if e.stderr: + logger.error(f"stderr: {e.stderr}") + raise RuntimeError( + f"{operation_name} failed. Check that Docker is running and " + f"the flang-wasm image is available." + ) from e + except FileNotFoundError: + raise RuntimeError( + "Docker not found. Please install Docker to build WASM.\n" + "See: https://docs.docker.com/get-docker/" + ) - return array(grid) +def build_wasm(options): + """Build WebAssembly binary from complete set of source files. -def p_grid(p_min, p_max, n): - from numpy import logspace, log10 + Uses flang-wasm Docker container to compile Fortran to WebAssembly. + Applies necessary WASM-specific modifications to source files. - """Return a list of n pressures between - p_min and p_max such that the grid of log(p) - is evenly spaced. + Args: + options: Export options object + + Returns: + dict: Paths to generated files {'js': path, 'wasm': path} + + Raises: + RuntimeError: If Docker not available or compilation fails + IOError: If required source files missing """ - p_minlog = log10(p_min) - p_maxlog = log10(p_max) + import shutil + from os.path import isfile, abspath - grid = logspace(p_minlog, p_maxlog, n) + logger.info("Building WebAssembly binary...") - return grid + # 1. Validate Docker availability + if not shutil.which("docker"): + raise RuntimeError( + "Docker not found. Please install Docker to build WASM.\n" + "See: https://docs.docker.com/get-docker/" + ) + # Get the current directory (should contain Fortran sources) + src_dir = abspath(".") -def timeit(func): - """ - Generic timing decorator + # 2. Check for required source files + required_files = ["kind_values.f90", "base.f90", "lattice.f90", "proclist.f90"] + missing = [f for f in required_files if not isfile(f)] + if missing: + raise IOError( + f"Required source files missing: {', '.join(missing)}\n" + f"Make sure you're in a directory with exported Fortran sources." + ) - To stop time for function call f - just :: - from kmos.utils import timeit - @timeit - def f(): - ... + # 3. Apply WASM-specific modifications + logger.info("Applying WASM-specific modifications...") + modifications_count = apply_wasm_modifications(src_dir) + if modifications_count > 0: + logger.info(f"Applied {modifications_count} WASM modifications") + # 4. Create c_bindings.f90 if it doesn't exist + if not isfile("c_bindings.f90"): + logger.info("Creating c_bindings.f90...") + create_c_bindings() + else: + logger.info("Using existing c_bindings.f90") + + # 5. Compile with flang-wasm + logger.info("Compiling Fortran sources with flang-wasm...") + + # First, compile each Fortran file to object files + # Use /src (container path) instead of host path + compile_cmd = """ + cd /src && + /opt/flang/host/bin/flang -g -O0 -target wasm32-unknown-emscripten -c kind_values.f90 -o kind_values.o && + /opt/flang/host/bin/flang -g -O0 -target wasm32-unknown-emscripten -c base.f90 -I. -o base.o && + /opt/flang/host/bin/flang -g -O0 -target wasm32-unknown-emscripten -c lattice.f90 -I. -o lattice.o && + /opt/flang/host/bin/flang -g -O0 -target wasm32-unknown-emscripten -c proclist.f90 -I. -o proclist.o && + /opt/flang/host/bin/flang -g -O0 -target wasm32-unknown-emscripten -c c_bindings.f90 -I. -o c_bindings.o """ - def wrapper(*args, **kwargs): - time0 = time() - func(*args, **kwargs) - logger.info("Executing %s took %.3f s" % (func.__name__, time() - time0)) + _run_docker_command(compile_cmd, src_dir, "Compilation", timeout=1800) + + # 6. Link with emcc + logger.info("Linking with emscripten...") + logger.info("Note: WASM linking can take several minutes for large models...") + + # Use /src (container path) instead of host path + # Production build: O3 optimization, no debug, no assertions + link_cmd = """ + cd /src && + /opt/emsdk/upstream/emscripten/emcc -O3 -sASSERTIONS=0 -sSTACK_SIZE=10MB \ + -sINITIAL_MEMORY=64MB -sALLOW_MEMORY_GROWTH \ + -sEXPORTED_FUNCTIONS=_kmc_init,_kmc_do_step,_kmc_get_time,_kmc_get_species,_kmc_get_system_size,_kmc_set_rate,_kmc_get_nr_sites,_kmc_get_accum_rate,_kmc_get_rate,_kmc_update_accum_rate,_kmc_get_nr2lattice,_kmc_get_debug_last_proc,_kmc_get_debug_last_site,_malloc,_free \ + -sEXPORTED_RUNTIME_METHODS=ccall,cwrap,getValue,setValue \ + -L/opt/flang/wasm/lib -lFortranRuntime \ + kind_values.o base.o lattice.o proclist.o c_bindings.o \ + -o kmc_model.js + """ - return wrapper + _run_docker_command(link_cmd, src_dir, "Linking", timeout=1800) + # 7. Validate output files + js_file = abspath("kmc_model.js") + wasm_file = abspath("kmc_model.wasm") + + if not isfile(js_file): + raise RuntimeError("Build failed: kmc_model.js was not created") + if not isfile(wasm_file): + raise RuntimeError("Build failed: kmc_model.wasm was not created") + + # Check files are not empty + import os + + js_size = os.path.getsize(js_file) + wasm_size = os.path.getsize(wasm_file) + + if js_size == 0: + raise RuntimeError("Build failed: kmc_model.js is empty") + if wasm_size == 0: + raise RuntimeError("Build failed: kmc_model.wasm is empty") + + # 8. Log success with file sizes + js_size_kb = js_size / 1024 + wasm_size_mb = wasm_size / 1024 / 1024 + + logger.info("WebAssembly build completed successfully!") + logger.info(f"Generated kmc_model.js ({js_size_kb:.1f} KB)") + logger.info(f"Generated kmc_model.wasm ({wasm_size_mb:.2f} MB)") + + return {"js": js_file, "wasm": wasm_file} -def col_tuple2str(tup): - """Convenience function that turns a HTML type color - into a tuple of three float between 0 and 1 - """ - r, g, b = tup - b *= 255 - res = "#" - res += hex(int(255 * r))[-2:].replace("x", "0") - res += hex(int(255 * g))[-2:].replace("x", "0") - res += hex(int(255 * b))[-2:].replace("x", "0") - return res +def apply_wasm_modifications(src_dir): + """Apply WASM-specific modifications to Fortran source files. + Key modifications for WASM compatibility: + 1. proclist.f90: + - Replace random_number() with custom RNG (WASM doesn't support Fortran intrinsics) + - Replace array slice assignments with explicit element assignments + - Add debug variables for tracking + - Change seed_size to 1 (simple seed) + 2. base.f90: + - Remove CHARACTER assignments (causes _FortranAAssign crashes) + - Replace array operations with explicit DO loops + 3. lattice.f90: + - Convert calculate_nr2lattice from function to subroutine (WASM can't return arrays) + - Replace array operations with explicit assignments + - Comment out validation loops (trigger 'unreachable' errors) -def col_str2tuple(hex_string): - """Convenience function that turns a HTML type color - into a tuple of three float between 0 and 1 + Args: + src_dir: Directory containing Fortran source files + + Returns: + int: Number of modifications applied """ - import gtk + import re + from os.path import join, isfile - try: - color = gtk.gdk.Color(hex_string) - except ValueError: - raise UserWarning( - "GTK cannot decipher color string {hex_string}: {e}".format(**locals()) + modifications_count = 0 + + # ======================================================================== + # PART 1: Modify proclist.f90 + # ======================================================================== + proclist_file = join(src_dir, "proclist.f90") + if not isfile(proclist_file): + logger.debug( + f"proclist.f90 not found in {src_dir}, skipping WASM modifications" + ) + return 0 + + with open(proclist_file, "r") as f: + proclist_content = f.read() + + # Skip if already modified (idempotent) + if "! WASM: Custom RNG" in proclist_content: + logger.debug("proclist.f90 already contains WASM modifications, skipping") + return 0 + + logger.info("Applying WASM modifications to proclist.f90...") + + # 1.1: Add get_volume to imports from base module if not already there + # Look for increment_procstat in the use base statement + if ( + "get_volume" not in proclist_content.split("use lattice")[0] + ): # Only check base module imports + # Find and replace increment_procstat (last item in use base) with increment_procstat, get_volume + proclist_content = re.sub( + r"(use base,.*?increment_procstat)(\s*\n)", + r"\1, &\n get_volume\2", + proclist_content, + flags=re.DOTALL, + count=1, ) - return (color.red_float, color.green_float, color.blue_float) + # 1.2: Add debug variables after implicit none + implicit_none_pattern = r"(implicit none\s*\n)" + debug_vars = r"\1\n! Debug variables to store last selected process and site\ninteger(kind=iint), public :: debug_last_proc_nr = 0\ninteger(kind=iint), public :: debug_last_nr_site = 0\n" + proclist_content = re.sub( + implicit_none_pattern, debug_vars, proclist_content, count=1 + ) -def jmolcolor_in_hex(i): - """Return a given jmol color in hexadecimal representation.""" - from ase.data.colors import jmol_colors + # 1.3: Change seed_size from 12 to 1 + proclist_content = re.sub( + r"integer\(kind=iint\), public :: seed_size = \d+", + r"integer(kind=iint), public :: seed_size = 1", + proclist_content, + ) - color = [int(x) for x in 255 * jmol_colors[i]] - r, g, b = color - a = 255 - color = (r << 24) | (g << 16) | (b << 8) | a - return color + # 1.4: Add custom RNG state variable after seed_arr declaration + seed_arr_pattern = ( + r"(integer\(kind=iint\), public, dimension\(:\), allocatable :: seed_arr.*\n)" + ) + rng_state = r"\1\n! WASM: Custom RNG state (module-level variable)\ninteger(kind=iint) :: rng_state = 123456789_iint\n" + proclist_content = re.sub(seed_arr_pattern, rng_state, proclist_content, count=1) + + # 1.5: Add custom RNG implementation after "contains" + contains_pattern = r"(contains\s*\n)" + custom_rng = r"""\1 +! WASM: Custom RNG to replace random_number() intrinsic +! Simple Linear Congruential Generator (LCG) +! Uses the same algorithm as glibc: X_n+1 = (a * X_n + c) mod m +! where a = 1103515245, c = 12345, m = 2^31 +subroutine custom_random_seed(seed_val) + integer(kind=iint), intent(in) :: seed_val + rng_state = seed_val +end subroutine custom_random_seed + +subroutine custom_random_number(harvest) + real(kind=rsingle), intent(out) :: harvest + integer(kind=iint), parameter :: a = 1103515245_iint + integer(kind=iint), parameter :: c = 12345_iint + integer(kind=iint), parameter :: m = 2147483647_iint ! 2^31 - 1 (Mersenne prime, safer) + integer(kind=iint) :: temp + + ! Update state with careful modulo to avoid overflow issues + temp = mod(rng_state, m) + temp = mod(a * temp, m) + temp = mod(temp + c, m) + if (temp < 0) temp = temp + m ! Ensure positive + rng_state = temp + + ! Convert to [0, 1) range + harvest = real(rng_state, kind=rsingle) / real(m, kind=rsingle) +end subroutine custom_random_number +""" + proclist_content = re.sub(contains_pattern, custom_rng, proclist_content, count=1) -def evaluate_template(template, escape_python=False, **kwargs): - """Very simple template evaluation function using only exec and str.format() + # 1.6: Replace all random_number() calls with custom_random_number() + proclist_content = re.sub( + r"\bcall random_number\(", r"call custom_random_number(", proclist_content + ) - There are two flavors of the template language, depending on wether - the python parts or the template parts are escaped. + # 1.7: Replace random_seed initialization with custom_random_seed + # Find and replace the entire random seed initialization block + seed_init_pattern = r"(\s+)(! initialize random number generator\s*\n\s*allocate\(seed_arr\(seed_size\)\)\s*\n\s*seed = seed_in\s*\n\s*seed_arr = seed\s*\n\s*call random_seed\(seed_size\)\s*\n\s*call random_seed\(put=seed_arr\)\s*\n\s*deallocate\(seed_arr\))" + seed_init_replacement = r"\1! WASM: Use custom RNG instead of Fortran intrinsics due to ABI incompatibilities\n\1seed = seed_in\n\1call custom_random_seed(seed_in)" + proclist_content = re.sub( + seed_init_pattern, seed_init_replacement, proclist_content + ) - A template can use the full python syntax. Every line starts with '#@ ' - is interpreted as a template line. Please use proper indentation before - and note the space after '#@'. + # 1.8: Add debug variable storage in do_kmc_step + # Find the do_kmc_step subroutine and add debug storage before run_proc_nr + run_proc_pattern = r"(\s+)(call run_proc_nr\(proc_nr, nr_site\))" + debug_storage = r"\1! Store for debugging (used by exported functions)\n\1debug_last_proc_nr = proc_nr\n\1debug_last_nr_site = nr_site\n\n\1\2" + proclist_content = re.sub(run_proc_pattern, debug_storage, proclist_content) + + # 1.9: Replace array slice assignments with explicit element assignments + slice_pattern = r"(\s+)lsite = nr2lattice\(nr_site, :\)" + explicit_replacement = r"\1! lsite = lattice_site, (vs. scalar site)\n\1! WASM: Explicit element assignment instead of array slice\n\n\1lsite(1) = nr2lattice(nr_site, 1)\n\n\1lsite(2) = nr2lattice(nr_site, 2)\n\n\1lsite(3) = nr2lattice(nr_site, 3)\n\n\1lsite(4) = nr2lattice(nr_site, 4)\n" + matches = re.findall(slice_pattern, proclist_content) + if len(matches) > 0: + proclist_content = re.sub(slice_pattern, explicit_replacement, proclist_content) + modifications_count += len(matches) + + # Write modified proclist.f90 + with open(proclist_file, "w") as f: + f.write(proclist_content) + + logger.info(f"Applied {modifications_count + 7} modifications to proclist.f90") + + # ======================================================================== + # PART 2: Modify base.f90 + # ======================================================================== + base_file = join(src_dir, "base.f90") + if isfile(base_file): + with open(base_file, "r") as f: + base_content = f.read() + + # Skip if already modified + if ( + "! WASM: Leave uninitialized" in base_content + or "! WASM: Skip CHARACTER assignment" in base_content + ): + logger.debug("base.f90 already contains WASM modifications") + else: + logger.info("Applying WASM modifications to base.f90...") - The template lines are converted to TEMPLATE_LINE.format(locals()) - and thefore every variable in the template line should be escape - with {}. + # 2.1: Leave system_name uninitialized (avoid _FortranAAssign crash) + base_content = re.sub( + r"character\(len=200\) :: system_name\s*\n", + r"character(len=200) :: system_name ! WASM: Leave uninitialized to avoid _FortranAAssign crash\n", + base_content, + ) + + # 2.2: Comment out system_name assignment in allocate_system + base_content = re.sub( + r"(\s+)(system_name = input_system_name\s*\n)", + r"\1! system_name = input_system_name ! WASM: Skip CHARACTER assignment - causes ABI crash in _FortranAAssign\n", + base_content, + ) + + # 2.3: Add WASM comment after variable declaration in update_accum_rate + # Find update_accum_rate subroutine and add comment after integer(kind=iint) :: i + base_content = re.sub( + r"(subroutine update_accum_rate.*?\n\s+integer\(kind=iint\) :: i\s*\n)", + r"\1\n ! WASM: Can't use print statements with INTEGER(8) due to ABI mismatch\n ! Use kmc_get_nr_sites() from JavaScript to debug instead\n", + base_content, + flags=re.DOTALL, + ) + + # 2.4: Replace array operations with explicit DO loops in allocate_system + # Add loop variables (i already exists, just add i_proc, i_vol, i_dim) + base_content = re.sub( + r"(subroutine allocate_system.*?integer\(kind=iint\), intent\(in\) :: input_volume, input_nr_of_proc\s*\n\s+logical :: system_allocated\s*\n)", + r"\1 integer(kind=iint) :: i_proc, i_vol, i_dim ! WASM: Loop variables for explicit init\n", + base_content, + flags=re.DOTALL, + ) - A valid template could be + # Replace avail_sites = 0 with explicit loops + base_content = re.sub( + r"(\s+)allocate\(avail_sites\(nr_of_proc, volume, 2\)\)\s*\n\s+avail_sites = 0", + r"\1allocate(avail_sites(nr_of_proc, volume, 2))\n" + r"\1do i_dim = 1, 2\n" + r"\1 do i_vol = 1, volume\n" + r"\1 do i_proc = 1, nr_of_proc\n" + r"\1 avail_sites(i_proc, i_vol, i_dim) = 0\n" + r"\1 end do\n" + r"\1 end do\n" + r"\1end do", + base_content, + ) + + # Replace lattice = null_species with explicit loop + base_content = re.sub( + r"(\s+)allocate\(lattice\(volume\)\)\s*\n\s+lattice = null_species", + r"\1allocate(lattice(volume))\n" + r"\1do i_vol = 1, volume\n" + r"\1 lattice(i_vol) = null_species\n" + r"\1end do", + base_content, + ) + + # Replace nr_of_sites = 0 with explicit loop + base_content = re.sub( + r"(\s+)allocate\(nr_of_sites\(nr_of_proc\)\)\s*\n\s+nr_of_sites = 0", + r"\1allocate(nr_of_sites(nr_of_proc))\n" + r"\1do i_proc = 1, nr_of_proc\n" + r"\1 nr_of_sites(i_proc) = 0\n" + r"\1end do", + base_content, + ) + + # Replace rates = 0 with explicit loop + base_content = re.sub( + r"(\s+)allocate\(rates\(nr_of_proc\)\)\s*\n\s+rates = 0", + r"\1allocate(rates(nr_of_proc))\n" + r"\1do i_proc = 1, nr_of_proc\n" + r"\1 rates(i_proc) = 0\n" + r"\1end do", + base_content, + ) + + # Replace accum_rates = 0 with explicit loop + base_content = re.sub( + r"(\s+)allocate\(accum_rates\(nr_of_proc\)\)\s*\n\s+accum_rates = 0", + r"\1allocate(accum_rates(nr_of_proc))\n" + r"\1do i_proc = 1, nr_of_proc\n" + r"\1 accum_rates(i_proc) = 0\n" + r"\1end do", + base_content, + ) + + # Replace integ_rates = 0 with explicit loop + base_content = re.sub( + r"(\s+)allocate\(integ_rates\(nr_of_proc\)\)\s*\n\s+integ_rates = 0", + r"\1allocate(integ_rates(nr_of_proc))\n" + r"\1do i_proc = 1, nr_of_proc\n" + r"\1 integ_rates(i_proc) = 0\n" + r"\1end do", + base_content, + ) + + # Write modified base.f90 + with open(base_file, "w") as f: + f.write(base_content) + + modifications_count += 8 + logger.info("Applied 8 modifications to base.f90") + + # ======================================================================== + # PART 3: Modify lattice.f90 + # ======================================================================== + lattice_file = join(src_dir, "lattice.f90") + if isfile(lattice_file): + with open(lattice_file, "r") as f: + lattice_content = f.read() + + # Skip if already modified + if "subroutine calculate_nr2lattice" in lattice_content: + logger.debug("lattice.f90 already has calculate_nr2lattice as subroutine") + else: + logger.debug("Applying WASM modifications to lattice.f90...") + + # 1. Convert function to subroutine + # First, extract the function to modify it + func_start = lattice_content.find("pure function calculate_nr2lattice(nr)") + func_end = lattice_content.find("end function calculate_nr2lattice") + len( + "end function calculate_nr2lattice" + ) - for i in range: - #@ Hello World {i} + if func_start != -1 and func_end != -1: + before = lattice_content[:func_start] + func_body = lattice_content[func_start:func_end] + after = lattice_content[func_end:] + # Replace function declaration with subroutine + func_body = func_body.replace( + "pure function calculate_nr2lattice(nr)", + "subroutine calculate_nr2lattice(nr, coords)", + ) + func_body = func_body.replace( + " integer(kind=iint), dimension(4) :: calculate_nr2lattice", + " integer(kind=iint), dimension(4), intent(out) :: coords", + ) + # Replace calculate_nr2lattice( with coords( in assignment statements only + # (not in the subroutine declaration or end statement) + lines = func_body.split("\n") + new_lines = [] + for line in lines: + if "subroutine" not in line.lower() and "coords(" not in line: + line = line.replace("calculate_nr2lattice(", "coords(") + new_lines.append(line) + func_body = "\n".join(new_lines) + + func_body = func_body.replace( + "end function calculate_nr2lattice", + "end subroutine calculate_nr2lattice", + ) + + lattice_content = before + func_body + after + else: + logger.warning( + "Could not find calculate_nr2lattice function to convert" + ) + + # 2. Add temp_coords variable + lattice_content = lattice_content.replace( + " integer(kind=iint) :: volume\n", + " integer(kind=iint) :: volume\n integer(kind=iint), dimension(4) :: temp_coords ! WASM\n", + ) + + # 2.5: Replace system_size array constructor with explicit assignments + # Avoid product() and array constructors due to ABI mismatch + lattice_content = re.sub( + r"(\s+)system_size = \(/input_system_size\(1\), input_system_size\(2\), 1/\)\s*\n\s+volume = system_size\(1\)\*system_size\(2\)\*system_size\(3\)\*spuck", + r"\1system_size(1) = input_system_size(1)\n\1system_size(2) = input_system_size(2)\n\1system_size(3) = 1\n\1! Avoid product() due to ABI mismatch with WASM runtime\n\1volume = system_size(1)*system_size(2)*system_size(3)*spuck", + lattice_content, + ) + + # 3. Comment out first validation loop + # Find the block starting with "! Let's check" and ending with "end do" (4 levels deep) + lines = lattice_content.split("\n") + new_lines = [] + in_validation1 = False + validation1_depth = 0 + + for line in lines: + if "! Let's check if the works correctly" in line: + in_validation1 = True + new_lines.append( + " ! Validation checks disabled for WASM - these trigger 'unreachable' errors" + ) + new_lines.append(" ! Let's check if the works correctly, first") + continue + + if in_validation1: + if ( + "do k = 0" in line + or "do j = 0" in line + or "do i = 0" in line + or ("do nr = 1" in line and "do nr = 1, spuck" in line) + ): + validation1_depth += 1 + new_lines.append(" ! " + line.lstrip()) + elif "end do" in line and validation1_depth > 0: + new_lines.append(" ! " + line.lstrip()) + validation1_depth -= 1 + if validation1_depth == 0: + in_validation1 = False + else: + new_lines.append(" ! " + line.lstrip()) + else: + new_lines.append(line) + + lattice_content = "\n".join(new_lines) + + # 4. Comment out second validation loop + lines = lattice_content.split("\n") + new_lines = [] + in_validation2 = False + + for line in lines: + if ( + "do check_nr=1," in line + and "calculate_lattice2nr(calculate_nr2lattice" in lattice_content + ): + # Check if this is the validation loop (next line has if statement) + idx = lines.index(line) + if ( + idx + 1 < len(lines) + and "if(.not.check_nr.eq.calculate_lattice2nr" in lines[idx + 1] + ): + in_validation2 = True + new_lines.append(" ! " + line.lstrip()) + continue + + if in_validation2: + new_lines.append(" ! " + line.lstrip()) + if "end do" in line: + in_validation2 = False + else: + new_lines.append(line) + + lattice_content = "\n".join(new_lines) + + # 5. Replace array slice assignment (may be commented out in source) + # First try to match uncommented version + if ( + "do check_nr=1, system_size(1)*system_size(2)*system_size(3)*spuck" + in lattice_content + and "nr2lattice(check_nr, :) = calculate_nr2lattice(check_nr)" + in lattice_content + ): + lattice_content = re.sub( + r" do check_nr=1, system_size\(1\)\*system_size\(2\)\*system_size\(3\)\*spuck\n nr2lattice\(check_nr, :\) = calculate_nr2lattice\(check_nr\)\n end do", + """ ! WASM: Use temp variable and explicit element assignment instead of array slice + do check_nr=1, system_size(1)*system_size(2)*system_size(3)*spuck + call calculate_nr2lattice(check_nr, temp_coords) + nr2lattice(check_nr, 1) = temp_coords(1) + nr2lattice(check_nr, 2) = temp_coords(2) + nr2lattice(check_nr, 3) = temp_coords(3) + nr2lattice(check_nr, 4) = temp_coords(4) + end do""", + lattice_content, + ) + logger.debug( + "Replaced nr2lattice array slice with explicit assignments" + ) + # Try to match commented-out version (from kmos templates) + # Pattern with consistent 4-space indentation: + # ! do check_nr=1, ... + # ! nr2lattice(check_nr, :) = ... + # ! end do + has_check_nr = ( + "! do check_nr=1, system_size(1)*system_size(2)*system_size(3)*spuck" + in lattice_content + ) + has_nr2lattice = ( + "! nr2lattice(check_nr, :) = calculate_nr2lattice(check_nr)" + in lattice_content + ) + + if has_check_nr and has_nr2lattice: + # Match the specific 3-line commented block for nr2lattice initialization + before_sub = lattice_content + lattice_content = re.sub( + r" ! do check_nr=1, system_size\(1\)\*system_size\(2\)\*system_size\(3\)\*spuck\n ! nr2lattice\(check_nr, :\) = calculate_nr2lattice\(check_nr\)\n ! end do", + """ ! WASM: Use temp variable and explicit element assignment instead of array slice + do check_nr=1, system_size(1)*system_size(2)*system_size(3)*spuck + call calculate_nr2lattice(check_nr, temp_coords) + nr2lattice(check_nr, 1) = temp_coords(1) + nr2lattice(check_nr, 2) = temp_coords(2) + nr2lattice(check_nr, 3) = temp_coords(3) + nr2lattice(check_nr, 4) = temp_coords(4) + end do""", + lattice_content, + ) + if before_sub != lattice_content: + logger.info( + "Uncommented and fixed nr2lattice initialization loop (was commented out in template)" + ) + else: + logger.warning( + "Found commented nr2lattice pattern but regex substitution failed" + ) + else: + logger.warning( + f"Could not find nr2lattice initialization loop (has_check_nr={has_check_nr}, has_nr2lattice={has_nr2lattice})" + ) + + with open(lattice_file, "w") as f: + f.write(lattice_content) + + modifications_count += 1 + logger.debug("Applied WASM modifications to lattice.f90") + + return modifications_count + + +def create_c_bindings(): + """Create c_bindings.f90 for WebAssembly exports. + + Generates Fortran bindings with C calling convention for JavaScript/WASM interface. + + Exported functions: + - kmc_init: Initialize KMC system + - kmc_do_step: Perform one KMC step + - kmc_get_time: Get current simulation time + - kmc_get_species: Get species at a site + - kmc_get_system_size: Get total number of sites + - kmc_set_rate: Set rate constant for a process + - kmc_get_nr_sites: Get number of sites + - kmc_get_accum_rate: Get accumulated rate + - kmc_get_rate: Get rate for a process + - kmc_update_accum_rate: Update accumulated rate + - kmc_get_nr2lattice: Get lattice coordinates from site number + - kmc_get_debug_last_proc: Get last executed process (debug) + - kmc_get_debug_last_site: Get last executed site (debug) + + Returns: + str: Absolute path to created c_bindings.f90 file + """ + from os.path import abspath + + # Use the WORKING c_bindings template from AB_model_wasm + c_bindings_code = """! C bindings for kMC functions to enable JavaScript interoperability +! This demonstrates how to make Fortran functions callable from WASM/JavaScript + +module c_bindings + use iso_c_binding + use kind_values + use base, only: get_kmc_time, set_rate_const, get_nrofsites, get_accum_rate, get_rate, update_accum_rate + use lattice, only: allocate_system, system_size, default, get_species, nr2lattice + use proclist, only: init, do_kmc_step + implicit none + +contains + + ! Initialize the kMC system + subroutine c_init(nx, ny, nz) bind(C, name="kmc_init") + use lattice, only: allocate_system + use proclist, only: initialize_state, nr_of_proc + integer(c_int), intent(in), value :: nx, ny, nz + integer(kind=iint) :: system_size_f(2) + integer(kind=iint) :: seed_in + character(len=200), parameter :: system_name = "kmc_system" ! WASM: Use parameter to avoid _FortranAAssign + + ! Convert C integers to Fortran integers (2D model, ignore nz) + system_size_f(1) = int(nx, kind=iint) + system_size_f(2) = int(ny, kind=iint) + + seed_in = 42 ! Random seed + + ! Call allocate_system and initialize_state directly to avoid optional parameter issues + call allocate_system(nr_of_proc, system_size_f, system_name) + call initialize_state(default, seed_in) + + end subroutine c_init + + ! Perform one kMC step + subroutine c_do_step() bind(C, name="kmc_do_step") + call do_kmc_step() + end subroutine c_do_step + + ! Get the current kMC time + function c_get_time() bind(C, name="kmc_get_time") result(time) + real(c_double) :: time + real(kind=rdouble) :: kmc_time_val + call get_kmc_time(kmc_time_val) + time = real(kmc_time_val, kind=c_double) + end function c_get_time + + ! Get species at a site (simplified - assumes default layer) + function c_get_species(x, y, z) bind(C, name="kmc_get_species") result(species) + integer(c_int), intent(in), value :: x, y, z + integer(c_int) :: species + integer(kind=iint) :: site(4) + + site(1) = int(x, kind=iint) + site(2) = int(y, kind=iint) + site(3) = int(z, kind=iint) + site(4) = default ! layer + + species = int(get_species(site), kind=c_int) + end function c_get_species + + ! Get system size + subroutine c_get_system_size(nx, ny, nz) bind(C, name="kmc_get_system_size") + integer(c_int), intent(out) :: nx, ny, nz + nx = int(system_size(1), kind=c_int) + ny = int(system_size(2), kind=c_int) + nz = int(system_size(3), kind=c_int) + end subroutine c_get_system_size + + ! Set rate constant for a process + subroutine c_set_rate_const(proc_nr, rate) bind(C, name="kmc_set_rate") + integer(c_int), intent(in), value :: proc_nr + real(c_double), intent(in), value :: rate + call set_rate_const(int(proc_nr, kind=iint), real(rate, kind=rdouble)) + end subroutine c_set_rate_const + + ! Get number of available sites for a process (for debugging) + function c_get_nr_sites(proc_nr) bind(C, name="kmc_get_nr_sites") result(nr_sites) + integer(c_int), intent(in), value :: proc_nr + integer(c_int) :: nr_sites + integer(kind=iint) :: nr_sites_f + call get_nrofsites(int(proc_nr, kind=iint), nr_sites_f) + nr_sites = int(nr_sites_f, kind=c_int) + end function c_get_nr_sites + + ! Get accumulated rate for a process (for debugging) + function c_get_accum_rate(proc_nr) bind(C, name="kmc_get_accum_rate") result(accum_rate) + integer(c_int), intent(in), value :: proc_nr + real(c_double) :: accum_rate + real(kind=rdouble) :: accum_rate_f + call get_accum_rate(int(proc_nr, kind=iint), accum_rate_f) + accum_rate = real(accum_rate_f, kind=c_double) + end function c_get_accum_rate + + ! Get rate for a process (for debugging) + function c_get_rate(proc_nr) bind(C, name="kmc_get_rate") result(rate) + integer(c_int), intent(in), value :: proc_nr + real(c_double) :: rate + real(kind=rdouble) :: rate_f + call get_rate(int(proc_nr, kind=iint), rate_f) + rate = real(rate_f, kind=c_double) + end function c_get_rate + + ! Manually trigger update_accum_rate (for debugging) + subroutine c_update_accum_rate() bind(C, name="kmc_update_accum_rate") + call update_accum_rate() + end subroutine c_update_accum_rate + + ! Debug: Get nr2lattice mapping for a given site number + ! Returns the coordinate at index coord_idx (1=x, 2=y, 3=z, 4=layer) + function c_get_nr2lattice(nr_site, coord_idx) bind(C, name="kmc_get_nr2lattice") result(coord_val) + integer(c_int), intent(in), value :: nr_site, coord_idx + integer(c_int) :: coord_val + integer(kind=iint) :: nr_site_f, coord_idx_f + + nr_site_f = int(nr_site, kind=iint) + coord_idx_f = int(coord_idx, kind=iint) + + coord_val = int(nr2lattice(nr_site_f, coord_idx_f), kind=c_int) + end function c_get_nr2lattice + + ! Debug: Get last selected process number from determine_procsite + ! Note: Debug variables not available in auto-generated code, return 0 + function c_get_debug_last_proc() bind(C, name="kmc_get_debug_last_proc") result(proc_nr) + integer(c_int) :: proc_nr + proc_nr = 0 + end function c_get_debug_last_proc + + ! Debug: Get last selected site number from determine_procsite + ! Note: Debug variables not available in auto-generated code, return 0 + function c_get_debug_last_site() bind(C, name="kmc_get_debug_last_site") result(nr_site) + integer(c_int) :: nr_site + nr_site = 0 + end function c_get_debug_last_site + +end module c_bindings +""" + + bindings_file = "c_bindings.f90" + + with open(bindings_file, "w") as f: + f.write(c_bindings_code) + + logger.debug(f"Created {bindings_file} with 13 exported C-callable functions") + + return abspath(bindings_file) + + +def T_grid(T_min, T_max, n): + from numpy import linspace, array + + """Return a list of n temperatures between + T_min and T_max such that the grid of T^(-1) + is evenly spaced. """ + + T_min1 = T_min ** (-1.0) + T_max1 = T_max ** (-1.0) + + grid = list(linspace(T_max1, T_min1, n)) + grid.reverse() + grid = [x ** (-1.0) for x in grid] + + return array(grid) + + +def p_grid(p_min, p_max, n): + from numpy import logspace, log10 + + """Return a list of n pressures between + p_min and p_max such that the grid of log(p) + is evenly spaced. + """ + + return logspace(log10(p_min), log10(p_max), n) + + +def evaluate_template(template, **kwargs): + """Generates code from a template with some simple python preprocessing. + + Preprocessor lines are started with #@ and end with @#. Preprocessor + lines can use all variables from the calling scope. Result lines + can also use local variables defined in preprocessor lines. + """ + + escape_python = kwargs.pop("escape_python", False) + # Create a namespace dict for exec() - Python 3 requires this for variable modification namespace = dict(kwargs) namespace["result"] = "" @@ -649,7 +1405,7 @@ def evaluate_template(template, escape_python=False, **kwargs): python_lines = "" matched = False for line in lines: - if re.match("^\s*%s ?" % PREFIX, line): + if re.match(r"^\s*%s ?" % PREFIX, line): python_lines += line.lstrip()[3:] matched = True else: @@ -663,9 +1419,9 @@ def evaluate_template(template, escape_python=False, **kwargs): # second turn literary lines into write statements python_lines = "" for line in lines: - if re.match("^\s*%s " % PREFIX, line): + if re.match(r"^\s*%s " % PREFIX, line): python_lines += line.lstrip()[3:] - elif re.match("^\s*%s$" % PREFIX, line): + elif re.match(r"^\s*%s$" % PREFIX, line): python_lines += '%sresult += "\\n"\n' % ( " " * (len(line) - len(line.lstrip())) ) @@ -685,7 +1441,7 @@ def evaluate_template(template, escape_python=False, **kwargs): python_lines = "" matched = False for line in lines: - if re.match("\s*%s ?" % PREFIX, line): + if re.match(r"\s*%s ?" % PREFIX, line): python_lines += "%spass %s" % ( " " * (len(line) - len(line.lstrip())), line.lstrip(), @@ -701,12 +1457,12 @@ def evaluate_template(template, escape_python=False, **kwargs): # second turn literary lines into write statements python_lines = "" for line in lines: - if re.match("\s*%s " % PREFIX, line): + if re.match(r"\s*%s " % PREFIX, line): python_lines += '%sresult += ("""%s""".format(**dict(locals())))\n' % ( " " * (len(line) - len(line.lstrip())), line.lstrip()[3:], ) - elif re.match("\s*%s" % PREFIX, line): + elif re.match(r"\s*%s" % PREFIX, line): python_lines += '%sresult += "\\n"\n' % ( " " * (len(line) - len(line.lstrip())) ) diff --git a/pyproject.toml b/pyproject.toml index 47058de4..941c7c16 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -44,6 +44,7 @@ dependencies = [ "numpy>=2.0.2", "meson>=1.2.0", "meson-python>=0.15.0", + "pandas>=2.3.3", ] [project.urls] @@ -84,6 +85,20 @@ kmos = [ [tool.pytest.ini_options] # Tests now run deterministically without requiring PYTHONHASHSEED=0 # Example: PYTHONPATH=. uv run pytest tests/ +# Disable setuptools' vendored typeguard plugin to avoid conflicts +addopts = [ + "-p", "no:typeguard", +] + +# Ignore test_types directory - these are scripts that require pre-built Fortran modules +# and cause collection errors (floating-point exceptions when importing compiled modules) +norecursedirs = ["test_types"] + +markers = [ + "docker: tests requiring Docker - VERY SLOW, 5-10 min per test (deselect with '-m \"not docker\"')", + "slow: slow-running tests (typically >1 minute)", + "wasm: WebAssembly-related tests (fast unit tests are NOT marked docker)", +] [tool.bumpversion] current_version = "0.4.2" diff --git a/tests/export_test/reference_export/lattice.f90 b/tests/export_test/reference_export/lattice.f90 index c3a55578..5ff214b4 100644 --- a/tests/export_test/reference_export/lattice.f90 +++ b/tests/export_test/reference_export/lattice.f90 @@ -231,7 +231,7 @@ subroutine allocate_system(nr_of_proc, input_system_size, system_name) end do end do - do check_nr=1, product(system_size)*spuck + do check_nr=1, system_size(1)*system_size(2)*system_size(3)*spuck if(.not.check_nr.eq.calculate_lattice2nr(calculate_nr2lattice(check_nr)))then print *, "ERROR in Mapping:", check_nr print *, "was mapped on", calculate_nr2lattice(check_nr) @@ -240,12 +240,12 @@ subroutine allocate_system(nr_of_proc, input_system_size, system_name) endif end do - allocate(nr2lattice(1:product(system_size)*spuck,4)) + allocate(nr2lattice(1:system_size(1)*system_size(2)*system_size(3)*spuck,4)) allocate(lattice2nr(-system_size(1):2*system_size(1)-1, & -system_size(2):2*system_size(2)-1, & -system_size(3):2*system_size(3)-1, & 1:spuck)) - do check_nr=1, product(system_size)*spuck + do check_nr=1, system_size(1)*system_size(2)*system_size(3)*spuck nr2lattice(check_nr, :) = calculate_nr2lattice(check_nr) end do do k = -system_size(3), 2*system_size(3)-1 diff --git a/tests/export_test/reference_export_intZGB_otf/lattice.f90 b/tests/export_test/reference_export_intZGB_otf/lattice.f90 index 37c27f99..9705d6f7 100644 --- a/tests/export_test/reference_export_intZGB_otf/lattice.f90 +++ b/tests/export_test/reference_export_intZGB_otf/lattice.f90 @@ -231,7 +231,7 @@ subroutine allocate_system(nr_of_proc, input_system_size, system_name) end do end do - do check_nr=1, product(system_size)*spuck + do check_nr=1, system_size(1)*system_size(2)*system_size(3)*spuck if(.not.check_nr.eq.calculate_lattice2nr(calculate_nr2lattice(check_nr)))then print *, "ERROR in Mapping:", check_nr print *, "was mapped on", calculate_nr2lattice(check_nr) @@ -240,12 +240,12 @@ subroutine allocate_system(nr_of_proc, input_system_size, system_name) endif end do - allocate(nr2lattice(1:product(system_size)*spuck,4)) + allocate(nr2lattice(1:system_size(1)*system_size(2)*system_size(3)*spuck,4)) allocate(lattice2nr(-system_size(1):2*system_size(1)-1, & -system_size(2):2*system_size(2)-1, & -system_size(3):2*system_size(3)-1, & 1:spuck)) - do check_nr=1, product(system_size)*spuck + do check_nr=1, system_size(1)*system_size(2)*system_size(3)*spuck nr2lattice(check_nr, :) = calculate_nr2lattice(check_nr) end do do k = -system_size(3), 2*system_size(3)-1 diff --git a/tests/export_test/reference_export_lat_int/lattice.f90 b/tests/export_test/reference_export_lat_int/lattice.f90 index c3a55578..5ff214b4 100644 --- a/tests/export_test/reference_export_lat_int/lattice.f90 +++ b/tests/export_test/reference_export_lat_int/lattice.f90 @@ -231,7 +231,7 @@ subroutine allocate_system(nr_of_proc, input_system_size, system_name) end do end do - do check_nr=1, product(system_size)*spuck + do check_nr=1, system_size(1)*system_size(2)*system_size(3)*spuck if(.not.check_nr.eq.calculate_lattice2nr(calculate_nr2lattice(check_nr)))then print *, "ERROR in Mapping:", check_nr print *, "was mapped on", calculate_nr2lattice(check_nr) @@ -240,12 +240,12 @@ subroutine allocate_system(nr_of_proc, input_system_size, system_name) endif end do - allocate(nr2lattice(1:product(system_size)*spuck,4)) + allocate(nr2lattice(1:system_size(1)*system_size(2)*system_size(3)*spuck,4)) allocate(lattice2nr(-system_size(1):2*system_size(1)-1, & -system_size(2):2*system_size(2)-1, & -system_size(3):2*system_size(3)-1, & 1:spuck)) - do check_nr=1, product(system_size)*spuck + do check_nr=1, system_size(1)*system_size(2)*system_size(3)*spuck nr2lattice(check_nr, :) = calculate_nr2lattice(check_nr) end do do k = -system_size(3), 2*system_size(3)-1 diff --git a/tests/export_test/reference_export_otf/lattice.f90 b/tests/export_test/reference_export_otf/lattice.f90 index 42ea0bfe..d9e58ef6 100644 --- a/tests/export_test/reference_export_otf/lattice.f90 +++ b/tests/export_test/reference_export_otf/lattice.f90 @@ -232,7 +232,7 @@ subroutine allocate_system(nr_of_proc, input_system_size, system_name) end do end do - do check_nr=1, product(system_size)*spuck + do check_nr=1, system_size(1)*system_size(2)*system_size(3)*spuck if(.not.check_nr.eq.calculate_lattice2nr(calculate_nr2lattice(check_nr)))then print *, "ERROR in Mapping:", check_nr print *, "was mapped on", calculate_nr2lattice(check_nr) @@ -241,12 +241,12 @@ subroutine allocate_system(nr_of_proc, input_system_size, system_name) endif end do - allocate(nr2lattice(1:product(system_size)*spuck,4)) + allocate(nr2lattice(1:system_size(1)*system_size(2)*system_size(3)*spuck,4)) allocate(lattice2nr(-system_size(1):2*system_size(1)-1, & -system_size(2):2*system_size(2)-1, & -system_size(3):2*system_size(3)-1, & 1:spuck)) - do check_nr=1, product(system_size)*spuck + do check_nr=1, system_size(1)*system_size(2)*system_size(3)*spuck nr2lattice(check_nr, :) = calculate_nr2lattice(check_nr) end do do k = -system_size(3), 2*system_size(3)-1 diff --git a/tests/export_test/reference_pdopd_lat_int/lattice.f90 b/tests/export_test/reference_pdopd_lat_int/lattice.f90 index 84fd1547..d7599e23 100644 --- a/tests/export_test/reference_pdopd_lat_int/lattice.f90 +++ b/tests/export_test/reference_pdopd_lat_int/lattice.f90 @@ -255,7 +255,7 @@ subroutine allocate_system(nr_of_proc, input_system_size, system_name) end do end do - do check_nr=1, product(system_size)*spuck + do check_nr=1, system_size(1)*system_size(2)*system_size(3)*spuck if(.not.check_nr.eq.calculate_lattice2nr(calculate_nr2lattice(check_nr)))then print *, "ERROR in Mapping:", check_nr print *, "was mapped on", calculate_nr2lattice(check_nr) @@ -264,12 +264,12 @@ subroutine allocate_system(nr_of_proc, input_system_size, system_name) endif end do - allocate(nr2lattice(1:product(system_size)*spuck,4)) + allocate(nr2lattice(1:system_size(1)*system_size(2)*system_size(3)*spuck,4)) allocate(lattice2nr(-system_size(1):2*system_size(1)-1, & -system_size(2):2*system_size(2)-1, & -system_size(3):2*system_size(3)-1, & 1:spuck)) - do check_nr=1, product(system_size)*spuck + do check_nr=1, system_size(1)*system_size(2)*system_size(3)*spuck nr2lattice(check_nr, :) = calculate_nr2lattice(check_nr) end do do k = -system_size(3), 2*system_size(3)-1 diff --git a/tests/export_test/reference_pdopd_local_smart/lattice.f90 b/tests/export_test/reference_pdopd_local_smart/lattice.f90 index 84fd1547..d7599e23 100644 --- a/tests/export_test/reference_pdopd_local_smart/lattice.f90 +++ b/tests/export_test/reference_pdopd_local_smart/lattice.f90 @@ -255,7 +255,7 @@ subroutine allocate_system(nr_of_proc, input_system_size, system_name) end do end do - do check_nr=1, product(system_size)*spuck + do check_nr=1, system_size(1)*system_size(2)*system_size(3)*spuck if(.not.check_nr.eq.calculate_lattice2nr(calculate_nr2lattice(check_nr)))then print *, "ERROR in Mapping:", check_nr print *, "was mapped on", calculate_nr2lattice(check_nr) @@ -264,12 +264,12 @@ subroutine allocate_system(nr_of_proc, input_system_size, system_name) endif end do - allocate(nr2lattice(1:product(system_size)*spuck,4)) + allocate(nr2lattice(1:system_size(1)*system_size(2)*system_size(3)*spuck,4)) allocate(lattice2nr(-system_size(1):2*system_size(1)-1, & -system_size(2):2*system_size(2)-1, & -system_size(3):2*system_size(3)-1, & 1:spuck)) - do check_nr=1, product(system_size)*spuck + do check_nr=1, system_size(1)*system_size(2)*system_size(3)*spuck nr2lattice(check_nr, :) = calculate_nr2lattice(check_nr) end do do k = -system_size(3), 2*system_size(3)-1 diff --git a/tests/export_test/test_export/kind_values.f90 b/tests/export_test/test_export/kind_values.f90 index bf6a50cb..f34ec976 100644 --- a/tests/export_test/test_export/kind_values.f90 +++ b/tests/export_test/test_export/kind_values.f90 @@ -7,10 +7,10 @@ module kind_values !****** implicit none -integer, parameter :: rsingle = SELECTED_REAL_KIND(p=15, r=200) -integer, parameter :: rdouble = SELECTED_REAL_KIND(p=15, r=200) -integer, parameter :: ibyte = SELECTED_INT_KIND(p=2) -integer, parameter :: ishort = SELECTED_INT_KIND(p=4) -integer, parameter :: iint = SELECTED_INT_KIND(p=9) -integer, parameter :: ilong = SELECTED_INT_KIND(p=18) +integer, parameter :: rsingle = SELECTED_REAL_KIND(15, 200) +integer, parameter :: rdouble = SELECTED_REAL_KIND(15, 200) +integer, parameter :: ibyte = SELECTED_INT_KIND(2) +integer, parameter :: ishort = SELECTED_INT_KIND(4) +integer, parameter :: iint = SELECTED_INT_KIND(9) +integer, parameter :: ilong = SELECTED_INT_KIND(18) end module kind_values diff --git a/tests/export_test/test_export/lattice.f90 b/tests/export_test/test_export/lattice.f90 index c3a55578..5ff214b4 100644 --- a/tests/export_test/test_export/lattice.f90 +++ b/tests/export_test/test_export/lattice.f90 @@ -231,7 +231,7 @@ subroutine allocate_system(nr_of_proc, input_system_size, system_name) end do end do - do check_nr=1, product(system_size)*spuck + do check_nr=1, system_size(1)*system_size(2)*system_size(3)*spuck if(.not.check_nr.eq.calculate_lattice2nr(calculate_nr2lattice(check_nr)))then print *, "ERROR in Mapping:", check_nr print *, "was mapped on", calculate_nr2lattice(check_nr) @@ -240,12 +240,12 @@ subroutine allocate_system(nr_of_proc, input_system_size, system_name) endif end do - allocate(nr2lattice(1:product(system_size)*spuck,4)) + allocate(nr2lattice(1:system_size(1)*system_size(2)*system_size(3)*spuck,4)) allocate(lattice2nr(-system_size(1):2*system_size(1)-1, & -system_size(2):2*system_size(2)-1, & -system_size(3):2*system_size(3)-1, & 1:spuck)) - do check_nr=1, product(system_size)*spuck + do check_nr=1, system_size(1)*system_size(2)*system_size(3)*spuck nr2lattice(check_nr, :) = calculate_nr2lattice(check_nr) end do do k = -system_size(3), 2*system_size(3)-1 diff --git a/tests/export_test/test_export_lat_int/kind_values.f90 b/tests/export_test/test_export_lat_int/kind_values.f90 index bf6a50cb..f34ec976 100644 --- a/tests/export_test/test_export_lat_int/kind_values.f90 +++ b/tests/export_test/test_export_lat_int/kind_values.f90 @@ -7,10 +7,10 @@ module kind_values !****** implicit none -integer, parameter :: rsingle = SELECTED_REAL_KIND(p=15, r=200) -integer, parameter :: rdouble = SELECTED_REAL_KIND(p=15, r=200) -integer, parameter :: ibyte = SELECTED_INT_KIND(p=2) -integer, parameter :: ishort = SELECTED_INT_KIND(p=4) -integer, parameter :: iint = SELECTED_INT_KIND(p=9) -integer, parameter :: ilong = SELECTED_INT_KIND(p=18) +integer, parameter :: rsingle = SELECTED_REAL_KIND(15, 200) +integer, parameter :: rdouble = SELECTED_REAL_KIND(15, 200) +integer, parameter :: ibyte = SELECTED_INT_KIND(2) +integer, parameter :: ishort = SELECTED_INT_KIND(4) +integer, parameter :: iint = SELECTED_INT_KIND(9) +integer, parameter :: ilong = SELECTED_INT_KIND(18) end module kind_values diff --git a/tests/export_test/test_export_lat_int/lattice.f90 b/tests/export_test/test_export_lat_int/lattice.f90 index c3a55578..5ff214b4 100644 --- a/tests/export_test/test_export_lat_int/lattice.f90 +++ b/tests/export_test/test_export_lat_int/lattice.f90 @@ -231,7 +231,7 @@ subroutine allocate_system(nr_of_proc, input_system_size, system_name) end do end do - do check_nr=1, product(system_size)*spuck + do check_nr=1, system_size(1)*system_size(2)*system_size(3)*spuck if(.not.check_nr.eq.calculate_lattice2nr(calculate_nr2lattice(check_nr)))then print *, "ERROR in Mapping:", check_nr print *, "was mapped on", calculate_nr2lattice(check_nr) @@ -240,12 +240,12 @@ subroutine allocate_system(nr_of_proc, input_system_size, system_name) endif end do - allocate(nr2lattice(1:product(system_size)*spuck,4)) + allocate(nr2lattice(1:system_size(1)*system_size(2)*system_size(3)*spuck,4)) allocate(lattice2nr(-system_size(1):2*system_size(1)-1, & -system_size(2):2*system_size(2)-1, & -system_size(3):2*system_size(3)-1, & 1:spuck)) - do check_nr=1, product(system_size)*spuck + do check_nr=1, system_size(1)*system_size(2)*system_size(3)*spuck nr2lattice(check_nr, :) = calculate_nr2lattice(check_nr) end do do k = -system_size(3), 2*system_size(3)-1 diff --git a/tests/export_test/test_pdopd_lat_int/kind_values.f90 b/tests/export_test/test_pdopd_lat_int/kind_values.f90 index bf6a50cb..f34ec976 100644 --- a/tests/export_test/test_pdopd_lat_int/kind_values.f90 +++ b/tests/export_test/test_pdopd_lat_int/kind_values.f90 @@ -7,10 +7,10 @@ module kind_values !****** implicit none -integer, parameter :: rsingle = SELECTED_REAL_KIND(p=15, r=200) -integer, parameter :: rdouble = SELECTED_REAL_KIND(p=15, r=200) -integer, parameter :: ibyte = SELECTED_INT_KIND(p=2) -integer, parameter :: ishort = SELECTED_INT_KIND(p=4) -integer, parameter :: iint = SELECTED_INT_KIND(p=9) -integer, parameter :: ilong = SELECTED_INT_KIND(p=18) +integer, parameter :: rsingle = SELECTED_REAL_KIND(15, 200) +integer, parameter :: rdouble = SELECTED_REAL_KIND(15, 200) +integer, parameter :: ibyte = SELECTED_INT_KIND(2) +integer, parameter :: ishort = SELECTED_INT_KIND(4) +integer, parameter :: iint = SELECTED_INT_KIND(9) +integer, parameter :: ilong = SELECTED_INT_KIND(18) end module kind_values diff --git a/tests/export_test/test_pdopd_lat_int/lattice.f90 b/tests/export_test/test_pdopd_lat_int/lattice.f90 index 84fd1547..d7599e23 100644 --- a/tests/export_test/test_pdopd_lat_int/lattice.f90 +++ b/tests/export_test/test_pdopd_lat_int/lattice.f90 @@ -255,7 +255,7 @@ subroutine allocate_system(nr_of_proc, input_system_size, system_name) end do end do - do check_nr=1, product(system_size)*spuck + do check_nr=1, system_size(1)*system_size(2)*system_size(3)*spuck if(.not.check_nr.eq.calculate_lattice2nr(calculate_nr2lattice(check_nr)))then print *, "ERROR in Mapping:", check_nr print *, "was mapped on", calculate_nr2lattice(check_nr) @@ -264,12 +264,12 @@ subroutine allocate_system(nr_of_proc, input_system_size, system_name) endif end do - allocate(nr2lattice(1:product(system_size)*spuck,4)) + allocate(nr2lattice(1:system_size(1)*system_size(2)*system_size(3)*spuck,4)) allocate(lattice2nr(-system_size(1):2*system_size(1)-1, & -system_size(2):2*system_size(2)-1, & -system_size(3):2*system_size(3)-1, & 1:spuck)) - do check_nr=1, product(system_size)*spuck + do check_nr=1, system_size(1)*system_size(2)*system_size(3)*spuck nr2lattice(check_nr, :) = calculate_nr2lattice(check_nr) end do do k = -system_size(3), 2*system_size(3)-1 diff --git a/tests/export_test/test_pdopd_local_smart/kind_values.f90 b/tests/export_test/test_pdopd_local_smart/kind_values.f90 index bf6a50cb..f34ec976 100644 --- a/tests/export_test/test_pdopd_local_smart/kind_values.f90 +++ b/tests/export_test/test_pdopd_local_smart/kind_values.f90 @@ -7,10 +7,10 @@ module kind_values !****** implicit none -integer, parameter :: rsingle = SELECTED_REAL_KIND(p=15, r=200) -integer, parameter :: rdouble = SELECTED_REAL_KIND(p=15, r=200) -integer, parameter :: ibyte = SELECTED_INT_KIND(p=2) -integer, parameter :: ishort = SELECTED_INT_KIND(p=4) -integer, parameter :: iint = SELECTED_INT_KIND(p=9) -integer, parameter :: ilong = SELECTED_INT_KIND(p=18) +integer, parameter :: rsingle = SELECTED_REAL_KIND(15, 200) +integer, parameter :: rdouble = SELECTED_REAL_KIND(15, 200) +integer, parameter :: ibyte = SELECTED_INT_KIND(2) +integer, parameter :: ishort = SELECTED_INT_KIND(4) +integer, parameter :: iint = SELECTED_INT_KIND(9) +integer, parameter :: ilong = SELECTED_INT_KIND(18) end module kind_values diff --git a/tests/export_test/test_pdopd_local_smart/lattice.f90 b/tests/export_test/test_pdopd_local_smart/lattice.f90 index 84fd1547..d7599e23 100644 --- a/tests/export_test/test_pdopd_local_smart/lattice.f90 +++ b/tests/export_test/test_pdopd_local_smart/lattice.f90 @@ -255,7 +255,7 @@ subroutine allocate_system(nr_of_proc, input_system_size, system_name) end do end do - do check_nr=1, product(system_size)*spuck + do check_nr=1, system_size(1)*system_size(2)*system_size(3)*spuck if(.not.check_nr.eq.calculate_lattice2nr(calculate_nr2lattice(check_nr)))then print *, "ERROR in Mapping:", check_nr print *, "was mapped on", calculate_nr2lattice(check_nr) @@ -264,12 +264,12 @@ subroutine allocate_system(nr_of_proc, input_system_size, system_name) endif end do - allocate(nr2lattice(1:product(system_size)*spuck,4)) + allocate(nr2lattice(1:system_size(1)*system_size(2)*system_size(3)*spuck,4)) allocate(lattice2nr(-system_size(1):2*system_size(1)-1, & -system_size(2):2*system_size(2)-1, & -system_size(3):2*system_size(3)-1, & 1:spuck)) - do check_nr=1, product(system_size)*spuck + do check_nr=1, system_size(1)*system_size(2)*system_size(3)*spuck nr2lattice(check_nr, :) = calculate_nr2lattice(check_nr) end do do k = -system_size(3), 2*system_size(3)-1 diff --git a/tests/test_cli_generated_model/kind_values.f90 b/tests/test_cli_generated_model/kind_values.f90 index bf6a50cb..f34ec976 100644 --- a/tests/test_cli_generated_model/kind_values.f90 +++ b/tests/test_cli_generated_model/kind_values.f90 @@ -7,10 +7,10 @@ module kind_values !****** implicit none -integer, parameter :: rsingle = SELECTED_REAL_KIND(p=15, r=200) -integer, parameter :: rdouble = SELECTED_REAL_KIND(p=15, r=200) -integer, parameter :: ibyte = SELECTED_INT_KIND(p=2) -integer, parameter :: ishort = SELECTED_INT_KIND(p=4) -integer, parameter :: iint = SELECTED_INT_KIND(p=9) -integer, parameter :: ilong = SELECTED_INT_KIND(p=18) +integer, parameter :: rsingle = SELECTED_REAL_KIND(15, 200) +integer, parameter :: rdouble = SELECTED_REAL_KIND(15, 200) +integer, parameter :: ibyte = SELECTED_INT_KIND(2) +integer, parameter :: ishort = SELECTED_INT_KIND(4) +integer, parameter :: iint = SELECTED_INT_KIND(9) +integer, parameter :: ilong = SELECTED_INT_KIND(18) end module kind_values diff --git a/tests/test_cli_generated_model/lattice.f90 b/tests/test_cli_generated_model/lattice.f90 index 1f544529..0d4435fb 100644 --- a/tests/test_cli_generated_model/lattice.f90 +++ b/tests/test_cli_generated_model/lattice.f90 @@ -221,7 +221,7 @@ subroutine allocate_system(nr_of_proc, input_system_size, system_name) end do end do - do check_nr=1, product(system_size)*spuck + do check_nr=1, system_size(1)*system_size(2)*system_size(3)*spuck if(.not.check_nr.eq.calculate_lattice2nr(calculate_nr2lattice(check_nr)))then print *, "ERROR in Mapping:", check_nr print *, "was mapped on", calculate_nr2lattice(check_nr) @@ -230,12 +230,12 @@ subroutine allocate_system(nr_of_proc, input_system_size, system_name) endif end do - allocate(nr2lattice(1:product(system_size)*spuck,4)) + allocate(nr2lattice(1:system_size(1)*system_size(2)*system_size(3)*spuck,4)) allocate(lattice2nr(-system_size(1):2*system_size(1)-1, & -system_size(2):2*system_size(2)-1, & -system_size(3):2*system_size(3)-1, & 1:spuck)) - do check_nr=1, product(system_size)*spuck + do check_nr=1, system_size(1)*system_size(2)*system_size(3)*spuck nr2lattice(check_nr, :) = calculate_nr2lattice(check_nr) end do do k = -system_size(3), 2*system_size(3)-1 diff --git a/tests/test_wasm.py b/tests/test_wasm.py new file mode 100644 index 00000000..864bc849 --- /dev/null +++ b/tests/test_wasm.py @@ -0,0 +1,556 @@ +""" +Comprehensive unit tests for WASM functionality. + +Tests cover: +- apply_wasm_modifications() - Source code modifications for WASM compatibility +- create_c_bindings() - C bindings generation +- build_wasm() - Full WASM compilation pipeline (requires Docker) +- Backend compatibility - All three backends (local_smart, lat_int, otf) +""" + +import pytest +import os +from pathlib import Path +from unittest.mock import patch, MagicMock + +# Import functions to test +from kmos.utils import apply_wasm_modifications, create_c_bindings, build_wasm + + +class TestApplyWasmModifications: + """Test WASM-specific source code modifications.""" + + def test_array_slice_replacement(self, tmp_path): + """Test that array slice assignments are replaced with explicit assignments.""" + # Create proclist.f90 with array slice pattern + proclist = tmp_path / "proclist.f90" + proclist_content = """ +subroutine do_something() + integer :: lsite(4) + integer :: nr_site + + ! This line should be modified + lsite = nr2lattice(nr_site, :) + +end subroutine +""" + proclist.write_text(proclist_content) + + # Apply modifications + apply_wasm_modifications(str(tmp_path)) + + # Verify replacement + content = proclist.read_text() + assert "! WASM: Explicit element assignment instead of array slice" in content + assert "lsite(1) = nr2lattice(nr_site, 1)" in content + assert "lsite(2) = nr2lattice(nr_site, 2)" in content + assert "lsite(3) = nr2lattice(nr_site, 3)" in content + assert "lsite(4) = nr2lattice(nr_site, 4)" in content + # Original line should be gone + assert "lsite = nr2lattice(nr_site, :)" not in content.replace("! WASM:", "") + + def test_multiple_array_slices(self, tmp_path): + """Test that multiple array slice patterns are all replaced.""" + proclist = tmp_path / "proclist.f90" + proclist_content = """ +subroutine proc1() + lsite = nr2lattice(nr_site, :) +end subroutine + +subroutine proc2() + lsite = nr2lattice(nr_site, :) +end subroutine +""" + proclist.write_text(proclist_content) + + apply_wasm_modifications(str(tmp_path)) + + content = proclist.read_text() + # Both occurrences should be replaced + assert content.count("lsite(1) = nr2lattice(nr_site, 1)") == 2 + assert content.count("lsite(2) = nr2lattice(nr_site, 2)") == 2 + + def test_idempotency(self, tmp_path): + """Test that modifications are idempotent (running twice doesn't double-modify).""" + proclist = tmp_path / "proclist.f90" + proclist_content = """ + lsite = nr2lattice(nr_site, :) +""" + proclist.write_text(proclist_content) + + # Apply once + apply_wasm_modifications(str(tmp_path)) + content_after_first = (tmp_path / "proclist.f90").read_text() + + # Apply again + apply_wasm_modifications(str(tmp_path)) + content_after_second = (tmp_path / "proclist.f90").read_text() + + # Content should be identical after both applications + assert content_after_first == content_after_second + # Should only have one WASM comment + assert content_after_second.count("! WASM:") == 1 + + def test_no_modification_when_already_modified(self, tmp_path): + """Test that files already containing WASM modifications are not modified again.""" + proclist = tmp_path / "proclist.f90" + proclist_content = """ + ! WASM: Explicit element assignment instead of array slice + lsite(1) = nr2lattice(nr_site, 1) + lsite(2) = nr2lattice(nr_site, 2) + lsite(3) = nr2lattice(nr_site, 3) + lsite(4) = nr2lattice(nr_site, 4) +""" + proclist.write_text(proclist_content) + + # Get content before + content_before = proclist.read_text() + + # Apply modifications + apply_wasm_modifications(str(tmp_path)) + + # Content should be unchanged + content_after = proclist.read_text() + assert content_before == content_after + + def test_no_pattern_found(self, tmp_path): + """Test behavior when no array slice patterns exist.""" + proclist = tmp_path / "proclist.f90" + proclist_content = """ +subroutine simple() + integer :: x + x = 5 +end subroutine +""" + proclist.write_text(proclist_content) + + # Should not raise error + apply_wasm_modifications(str(tmp_path)) + + # Content should be unchanged + content = proclist.read_text() + assert content == proclist_content + assert "! WASM:" not in content + + def test_missing_proclist_file(self, tmp_path): + """Test behavior when proclist.f90 doesn't exist.""" + # Should not raise error + apply_wasm_modifications(str(tmp_path)) + + def test_empty_proclist_file(self, tmp_path): + """Test behavior with empty proclist.f90.""" + proclist = tmp_path / "proclist.f90" + proclist.write_text("") + + # Should not raise error + apply_wasm_modifications(str(tmp_path)) + + # File should still be empty + assert proclist.read_text() == "" + + def test_preserves_indentation(self, tmp_path): + """Test that indentation is preserved in replacements.""" + proclist = tmp_path / "proclist.f90" + proclist_content = """ +subroutine test1() + lsite = nr2lattice(nr_site, :) +end subroutine + +subroutine test2() + lsite = nr2lattice(nr_site, :) +end subroutine + +subroutine test3() + lsite = nr2lattice(nr_site, :) +end subroutine +""" + proclist.write_text(proclist_content) + + apply_wasm_modifications(str(tmp_path)) + + content = proclist.read_text() + # Check that different indentation levels are preserved + assert " lsite(1) = nr2lattice(nr_site, 1)" in content + assert " lsite(1) = nr2lattice(nr_site, 1)" in content + assert " lsite(1) = nr2lattice(nr_site, 1)" in content + + +class TestCreateCBindings: + """Test C bindings generation.""" + + def test_bindings_file_created(self, tmp_path): + """Test that c_bindings.f90 is created.""" + os.chdir(tmp_path) + + create_c_bindings() + + bindings_file = tmp_path / "c_bindings.f90" + assert bindings_file.exists() + assert bindings_file.stat().st_size > 0 + + def test_bindings_content_structure(self, tmp_path): + """Test that c_bindings.f90 has the correct structure.""" + os.chdir(tmp_path) + + create_c_bindings() + + content = (tmp_path / "c_bindings.f90").read_text() + + # Check module declaration + assert "module c_bindings" in content + assert "use iso_c_binding" in content + assert "use kind_values" in content + assert "use base" in content + assert "use lattice" in content + assert "use proclist" in content + + # Check end module + assert "end module c_bindings" in content + + def test_required_functions_present(self, tmp_path): + """Test that all required C-callable functions are present.""" + os.chdir(tmp_path) + + create_c_bindings() + + content = (tmp_path / "c_bindings.f90").read_text() + + # List of required functions for WASM interface + required_functions = [ + "kmc_init", + "kmc_do_step", + "kmc_get_time", + "kmc_get_species", + "kmc_get_system_size", + "kmc_set_rate", + "kmc_get_nr_sites", + "kmc_get_accum_rate", + "kmc_get_rate", + "kmc_update_accum_rate", + "kmc_get_nr2lattice", + "kmc_get_debug_last_proc", + "kmc_get_debug_last_site", + ] + + for func_name in required_functions: + assert f'bind(C, name="{func_name}")' in content + + def test_overwrites_existing_bindings(self, tmp_path): + """Test that creating bindings overwrites existing c_bindings.f90.""" + os.chdir(tmp_path) + + # Create initial bindings + bindings_file = tmp_path / "c_bindings.f90" + bindings_file.write_text("old content") + + # Create new bindings + create_c_bindings() + + # Content should be new, not "old content" + content = bindings_file.read_text() + assert "old content" not in content + assert "module c_bindings" in content + + def test_bindings_fortran_syntax(self, tmp_path): + """Test that generated bindings have valid Fortran syntax.""" + os.chdir(tmp_path) + + create_c_bindings() + + content = (tmp_path / "c_bindings.f90").read_text() + + # Check for common Fortran syntax elements + assert "implicit none" in content + assert "contains" in content + assert "subroutine" in content or "function" in content + + +class TestBuildWasmValidation: + """Test WASM build validation (without Docker).""" + + def test_docker_not_available_error(self, tmp_path): + """Test that missing Docker gives clear error.""" + os.chdir(tmp_path) + + # Create minimal required files + for filename in ["kind_values.f90", "base.f90", "lattice.f90", "proclist.f90"]: + (tmp_path / filename).write_text("! minimal") + + options = MagicMock() + + # Mock shutil.which to return None (Docker not found) + with patch("shutil.which", return_value=None): + with pytest.raises(RuntimeError, match="Docker not found"): + build_wasm(options) + + def test_missing_files_error(self, tmp_path): + """Test that missing source files give clear error.""" + os.chdir(tmp_path) + + # Don't create any files + options = MagicMock() + + # Mock Docker as available + with patch("shutil.which", return_value="/usr/bin/docker"): + with pytest.raises(IOError, match="Required source files missing"): + build_wasm(options) + + +@pytest.mark.docker +@pytest.mark.slow +class TestBuildWasm: + """Test WASM compilation (requires Docker, very slow ~5 min per test).""" + + def create_minimal_fortran_sources(self, directory): + """Helper to create minimal Fortran source files for testing.""" + # Create kind_values.f90 + (directory / "kind_values.f90").write_text(""" +module kind_values + implicit none + integer, parameter :: iint = selected_int_kind(9) + integer, parameter :: rsingle = selected_real_kind(6, 37) + integer, parameter :: rdouble = selected_real_kind(15, 307) +end module kind_values +""") + + # Create base.f90 + (directory / "base.f90").write_text(""" +module base + use kind_values + implicit none + real(kind=rdouble), private :: kmc_time = 0.0d0 +contains + subroutine get_kmc_time(time) + real(kind=rdouble), intent(out) :: time + time = kmc_time + end subroutine + + subroutine get_kmc_step(step) + integer(kind=iint), intent(out) :: step + step = 0 + end subroutine +end module base +""") + + # Create lattice.f90 + (directory / "lattice.f90").write_text(""" +module lattice + use kind_values + implicit none + integer(kind=iint), dimension(3), public :: system_size = [10, 10, 1] + integer(kind=iint), public :: spuck = 1 + integer(kind=iint), dimension(:,:), allocatable :: nr2lattice +contains + subroutine allocate_system(n, system_name, layer, seed, no_banner) + integer(kind=iint), dimension(2), intent(in) :: n + character(len=*), intent(in) :: system_name + integer(kind=iint), intent(in) :: layer, seed + logical, intent(in) :: no_banner + system_size(1:2) = n + allocate(nr2lattice(100, 4)) + end subroutine + + function lattice2nr(lattice) + integer(kind=iint), dimension(4), intent(in) :: lattice + integer(kind=iint) :: lattice2nr + lattice2nr = 1 + end function + + function get_species(lattice) + integer(kind=iint), dimension(4), intent(in) :: lattice + integer(kind=iint) :: get_species + get_species = 0 + end function +end module lattice +""") + + # Create proclist.f90 + (directory / "proclist.f90").write_text(""" +module proclist + use kind_values + use base + use lattice + implicit none + integer(kind=iint), public :: nr_of_proc = 1 + integer(kind=iint), dimension(:), allocatable, public :: avail_sites +contains + subroutine init(system_size, system_name, layer, seed, no_banner) + integer(kind=iint), dimension(2), intent(in) :: system_size + character(len=*), intent(in) :: system_name + integer(kind=iint), intent(in) :: layer, seed + logical, intent(in) :: no_banner + call allocate_system(system_size, system_name, layer, seed, no_banner) + end subroutine + + subroutine do_kmc_step() + ! Minimal implementation + end subroutine + + subroutine set_rate_const(proc_nr, rate) + integer(kind=iint), intent(in) :: proc_nr + real(kind=rsingle), intent(in) :: rate + end subroutine +end module proclist +""") + + def test_build_success_with_valid_sources(self, tmp_path): + """Test successful WASM build with valid Fortran sources. + + WARNING: This test actually compiles WASM with Docker and takes ~5-10 minutes. + Run with: pytest -m docker + """ + os.chdir(tmp_path) + self.create_minimal_fortran_sources(tmp_path) + + # Create a mock options object + options = MagicMock() + + # This will actually try to run Docker + try: + build_wasm(options) + + # Check that output files were created + assert (tmp_path / "kmc_model.js").exists() + assert (tmp_path / "kmc_model.wasm").exists() + + # Check file sizes are reasonable (not empty, not too small) + js_size = (tmp_path / "kmc_model.js").stat().st_size + wasm_size = (tmp_path / "kmc_model.wasm").stat().st_size + + assert js_size > 1000, "JS file too small" + assert wasm_size > 1000, "WASM file too small" + + except Exception as e: + if "docker" in str(e).lower(): + pytest.skip(f"Docker not available: {e}") + else: + raise + + +@pytest.mark.wasm +@pytest.mark.parametrize("backend", ["local_smart", "lat_int", "otf"]) +class TestWasmBackendCompatibility: + """Test WASM export with all backends.""" + + def test_export_backend_creates_bindings(self, backend, tmp_path): + """Test that WASM export creates c_bindings.f90 for all backends.""" + # This test would typically use kmos export command + # For now, we test that create_c_bindings works regardless of backend + os.chdir(tmp_path) + + create_c_bindings() + + bindings_file = tmp_path / "c_bindings.f90" + assert bindings_file.exists() + + # Bindings content should be the same for all backends + content = bindings_file.read_text() + assert "module c_bindings" in content + + def test_modifications_work_for_backend(self, backend, tmp_path): + """Test that WASM modifications work regardless of backend.""" + # Create a proclist.f90 with array slice pattern + proclist = tmp_path / "proclist.f90" + proclist.write_text(""" + lsite = nr2lattice(nr_site, :) +""") + + # Apply modifications (backend-agnostic) + apply_wasm_modifications(str(tmp_path)) + + content = proclist.read_text() + assert "! WASM:" in content + assert "lsite(1) = nr2lattice(nr_site, 1)" in content + + +class TestWasmHelperFunctions: + """Test helper functions and edge cases.""" + + def test_apply_modifications_preserves_file_permissions(self, tmp_path): + """Test that file permissions are preserved after modifications.""" + proclist = tmp_path / "proclist.f90" + proclist.write_text(" lsite = nr2lattice(nr_site, :)") + + # Set specific permissions + os.chmod(proclist, 0o644) + original_mode = os.stat(proclist).st_mode + + apply_wasm_modifications(str(tmp_path)) + + # Permissions should be preserved + new_mode = os.stat(proclist).st_mode + assert original_mode == new_mode + + def test_create_bindings_in_current_directory(self, tmp_path): + """Test that bindings are created in current working directory.""" + # Change to tmp directory + original_dir = os.getcwd() + try: + # Clean up any existing c_bindings.f90 from previous tests + bindings_in_original = Path(original_dir) / "c_bindings.f90" + if bindings_in_original.exists(): + bindings_in_original.unlink() + + os.chdir(tmp_path) + + create_c_bindings() + + # File should be in tmp_path + assert (tmp_path / "c_bindings.f90").exists() + + finally: + os.chdir(original_dir) + + +@pytest.mark.docker +class TestWasmDockerIntegration: + """Test Docker-specific functionality.""" + + @pytest.mark.slow + def test_docker_image_availability(self): + """Test that the flang-wasm Docker image can be pulled.""" + import subprocess + + try: + result = subprocess.run( + ["docker", "pull", "ghcr.io/r-wasm/flang-wasm:main"], + capture_output=True, + timeout=300, # 5 minutes + ) + assert result.returncode == 0, "Failed to pull flang-wasm image" + except FileNotFoundError: + pytest.skip("Docker not installed") + except subprocess.TimeoutExpired: + pytest.fail("Docker pull timed out") + + def test_docker_command_construction(self, tmp_path): + """Test that Docker commands are constructed correctly.""" + # This is more of a smoke test to ensure build_wasm doesn't crash + # when constructing Docker commands + + os.chdir(tmp_path) + + # Create minimal sources + for filename in ["kind_values.f90", "base.f90", "lattice.f90", "proclist.f90"]: + (tmp_path / filename).write_text("! minimal") + + # Mock subprocess to check command construction + with patch("subprocess.run") as mock_run: + mock_run.return_value = MagicMock(returncode=0, stdout="", stderr="") + + options = MagicMock() + + try: + build_wasm(options) + except Exception: + pass # We're just testing command construction + + # Check that subprocess.run was called with docker commands + assert mock_run.called + first_call_args = mock_run.call_args_list[0][0][0] + assert first_call_args[0] == "docker" + assert "ghcr.io/r-wasm/flang-wasm:main" in first_call_args + + +if __name__ == "__main__": + pytest.main([__file__, "-v"])