EDA workflow

[metdyn]

Description

This is a first step to make EDA workflow run within SWELL framework. These 32 ensemble members from @rtodling run independently on each node. Both GSI B and all observation type from 3dvar-atmos are used here. The single node run for each member is made possible with the current setting: C90 resolution using thinned observations (rejecting 80% of obs at random).
The memory usage on each node is roughly 72% (There are still details to work out and improve).

The overall workflow:

The EDA-3D implementation here is a jedi block-EDA version. It is a very mechanical version of the EDA method.
A conventional 3D-Var minimization is invoked for each member state and the resulting ensemble anlysis states are final products of 3D-EDA. The observation perturbation term is added to the cost function when the member is a perturbed state.

1. The driver for this mechnism is in the .cylc file:
   

   swell/src/swell/suites/eda/flow.cylc

   Lines 79 to 81 in 35e846c

   	{% for i in range( 1, models[model_component]['ensemble_num_members'] + 1 ) %}
   	SP2-{{model_component}} => RunJediEdaExecutable_mem{{i}}-{{model_component}} => SP3-{{model_component}}
   	{% endfor %}

   
   

   swell/src/swell/suites/eda/flow.cylc

   Lines 170 to 179 in 35e846c

   	{% for i in range( 1, models[model_component]['ensemble_num_members'] + 1 ) %}
   	[[RunJediEdaExecutable_mem{{i}}-{{model_component}}]]
   	script = "swell task RunJediEdaExecutable $config -d $datetime -m {{model_component}} -imem {{i}}"
   	platform = {{platform}}
   	execution time limit = {{scheduling["RunJediEdaExecutable"]["execution_time_limit"]}}
   	[[[directives]]]
   	{%- for key, value in scheduling["RunJediEdaExecutable"]["directives"][model_component].items() %}
   	--{{key}} = {{value}}
   	{%- endfor %}
   	{% endfor %}

Each job submission completes one cost function minization for one member:
"swell task RunJediEdaExecutable $config -d $datetime -m {{model_component}} -imem {{i}}"

2. The observation thinning can be disabled or tuned by setting obs_thinning_rej_fraction from zero (no thinning) to 1 (reject all).
   

   swell/src/swell/suites/eda/suite_config.py

   Line 107 in 35e846c

   qd.obs_thinning_rej_fraction(0.8),

Thinning can also be turned off by commenting out this line:

swell/src/swell/suites/eda/flow.cylc

Line 76 in 35e846c

RenderJediObservations-{{model_component}} => RunJediObsfiltersExecutable-{{model_component}} => SP2-{{model_component}}

3. Each EDA member run uses its own copy of observation files to prevent HDF5 file lock when multiple jobs are accessing a single file. Hence cleaning up the directory for each 3D-Var run is important:
   

   swell/src/swell/suites/eda/flow.cylc

   Lines 88 to 89 in 35e846c

   	# CleanupEda directory
   	SP4-{{model_component}} => CleanEdaFiles-{{model_component}}

The diagnostic tools:

For diagnostic purpose, we have designed handy codelets to compute mean/variance and diffstate (for increment, e.g.) in swell.

A. To compute the mean and variance of a set of states including bkg and analysis,
user can define the file specs for a list containing multiple entries where each can contain
its own input file, output file, and grid_type.

There are two ways to create the configure file.

Method-1. modify the EDA suite config file:

swell/src/swell/suites/eda/suite_config.py

Lines 108 to 119 in 35e846c

	qd.ensmeanvariance_spec([
	{"state": "bkg",
	"fn_input": "ebkg/mem%mem%/geos.mem%mem%.%yyyy%mm%dd_%hh%MM%ssz.nc4",
	"fn_output_mean": "geos.prior.mean",
	"fn_output_variance": "geos.prior.variance",
	"grid_type": ['cs', 'latlon']},
	{"state": "analysis",
	"fn_input": "analysis/mem%mem%/eda.ana.mem%mem%.%yyyy%mm%dd_%hh%MM%ssz.nc4",
	"fn_output_mean": "eda.ana.mean",
	"fn_output_variance": "eda.ana.variance",
	"grid_type": ['cs', 'latlon']},
	]),

Method-2. create a section in override.yaml

models:
  geos_atmosphere:
    ensmeanvariance_spec:
    - state: bkg
      fn_input: ebkg/mem%mem%/geos.mem%mem%.%yyyy%mm%dd_%hh%MM%ssz.nc4
      fn_output_mean: geos.prior.mean
      fn_output_variance: geos.prior.variance
      grid_type: [ 'latlon', 'cs' ]
    - state: analysis
      fn_input: analysis/mem%mem%/eda.ana.mem%mem%.%yyyy%mm%dd_%hh%MM%ssz.nc4
      fn_output_mean: eda.ana.mean
      fn_output_variance: eda.ana.variance
      grid_type: [ 'latlon', 'cs' ]

and run: swell create eda_atmos -o override.yaml to override the config default

Finally to compute the mean and variance, run
swell task RunJediEnsembleMeanVariance PATHTO/experiment.yaml -d $date -m geos_atmosphere

B. To compute the difference between two states, similar approaches can be used.

Method-1. modify the EDA suite config file:

swell/src/swell/suites/eda/suite_config.py

Lines 120 to 128 in 35e846c

	qd.diffstates_spec({
	"state1":
	{"fn_input": "geos.prior.mean.%yyyy%mm%dd_%hh%MM%ssz.nc4"},
	"state2":
	{"fn_input": "eda.ana.mean.%yyyy%mm%dd_%hh%MM%ssz.nc4"},
	"state_diff":
	{"fn_output": "eda.mean-inc", "grid_type": ['cs', 'latlon']},
	"state_type": "ensemble"
	}),

Method-2. create an override section in override.yaml

models:
  geos_atmosphere:
    ensmeanvariance_spec:
    diffstates_spec:
      state1:
        fn_input: geos.prior.mean.%yyyy%mm%dd_%hh%MM%ssz.nc4
      state2:
        fn_input: eda.ana.mean.%yyyy%mm%dd_%hh%MM%ssz.nc4
      state_diff:
        fn_output: eda.mean-inc
        grid_type: [ 'cs', 'latlon' ]  
      state_type: ensemble

and run: swell create eda_atmos -o override.yaml to override the config default

Finally to compute the diff states, run
swell task RunJediDiffstates PATHTO/experiment.yaml -d $date -m geos_atmosphere

How to run this block EDA code:

swell create eda_atmos -o override.yaml  -p nccs_discover_sles15 -k
swell launch   your_eda_suite
[override.yaml contains experiment_root and experiment_id]

Example output: T mean increment at 500 hPa :

The work is advised by @rtodling.
@mranst makes the flow.cylc loop work for ensemble members.

Timing results

Resolution	Random obs. thinning	NX/NY	Nnode	Perhost	CPU time Per member (s)
C90	80%	4, 5	1	NA	347.0

In progress

Additional features:

1. In analogy to fv3jedi ctest for 3d-EDA, we have implemented a choice to use hybrid-B matrix in 3D-Var. Since this method is not tested, we have coded this function as an option for future use.

https://github.com/GEOS-ESM/swell/blob/35e846ce2371c9ac360ed13cf852082928bf922e/src/swell/configuration/jedi/interfaces/geos_atmosphere/model/background_error_eda_hybridB.py

A technical interesting part of the code is that I am using a string template in .py file, which is very useful here because it avoids writing 72 lines of pure real numbers into python code, and because we will need to tune and update this set of coefficient quite often later on, presumably. This is viewed as a transition back from python code to yaml style.

Issues

1. This code skips obs perturbation when the obs space contains covariance for R, because we did not find the corresponding source code in JEDI for this function.

swell/src/swell/tasks/run_jedi_eda_executable.py

Lines 159 to 173 in 35e846c

	if imember > 1:
	obs_cov_model = observer.get('obs error', {}).get('covariance model')
	print(f'{observation}: obs_cov_model = {obs_cov_model}')
	if obs_cov_model and 'cross variable covariances' in obs_cov_model:
	print(f"Found cross covariance obs: {obs_cov_model}, skip perturbation")
	else:
	print(f"No cross varaible covariance found: {observation}, Obs Error Diagonal")
	obs_error_dict = {
	'covariance model': 'diagonal',
	'zero-mean perturbations': True,
	'member': imember,
	'number of members': nmember
	}
	observer.update({'obs error': obs_error_dict})
	observer['obs space'].update({'obs perturbations seed': imember})

2. People should note that:
   Adding air_pressure_levels to the list of analysis variable will break the 3D-Var code with hybrid-B. (I assume there is an issue with bump code) So I did not include it in the interface model config setup for EDA. Here my assumption is that always having air_pressure_at_surface included in analysis variable can be sufficient.

To do list:

1. We can eliminate hofx output, increment output in the 3D-Var run for each member, but we want to compare with output from more advanced EDA method.

[mer-a-o]

I tried running this suite and EvaIncrement-geos_atmosphere task fails with this error:

FileNotFoundError: [Errno 2] No such file or directory: '/discover/nobackup/mabdiosk/SwellExperiments/swell-eda_atmos/run/20231010T000000Z/geos_atmosphere/swell-eda_atmos.increment-iter1.20231010_000000z.nc4'

I did not use any override file and simply ran:
swell create eda_atmos

and then launched the experiment.

[metdyn]

Hi All, as a rule the override file is not needed for PRs (excuse me for the confusion), otherwise the usage is diverging and thanks @jeromebarre for pointing this out. (My ensuring PRs will not contain it)

[metdyn]

@metdyn, Do you expect the observations to be fetched from R2D2 for this experiment?

I am unfamiliar with R2D2 so far, but it is definitely welcome to add more options. Currently as you may have found the EBKG path is hard coded in py and the obs are obtained from old script pointing to local IODA files on disk.

The flow here only implemented the single point Analysis without any cycle.

[metdyn]

from @jeromebarre: an incorrect r2d2 credential setup to access the data stored for this PR.
Thanks for pointing out this issue. I did not use r2d2 at all in this PR. We are focusing on running the plain EDA algorithm and to pave way for the next EDA Lanczos PR.

[metdyn]

I tried running this suite and EvaIncrement-geos_atmosphere task fails with this error:
FileNotFoundError: [Errno 2] No such file or directory: '/discover/nobackup/mabdiosk/SwellExperiments/swell-eda_atmos/run/20231010T000000Z/geos_atmosphere/swell-eda_atmos.increment-iter1.20231010_000000z.nc4'

Thanks @mer-a-o for testing. There was a typo "==>" issue pointed out from @jeromebarre (corrected) and I think there is a keyword match pattern in eva_increment.py that I found yesterday when running the next PR. I will fix this and run a test.

[jeromebarre]

@metdyn, Do you expect the observations to be fetched from R2D2 for this experiment?

I am unfamiliar with R2D2 so far, but it is definitely welcome to add more options. Currently as you may have found the EBKG path is hard coded in py and the obs are obtained from old script pointing to local IODA files on disk.

The flow here only implemented the single point Analysis without any cycle.

Currently the the suite without overrides fetches obs from r2d2...

[metdyn]

@mer-a-o, the workflow now proceeds to generate all the graphics files, hence the eva filename not found problem is solved.

[metdyn]

the suite without overrides fetches obs from r2d2...

That is nice and thanks for letting me know!

[Dooruk]

@metdyn, while it is not a strict requirement it would be nice to get this in before next week's SWELL documentation sprint. Do you think this is ready to be merged in soon?

[metdyn]

Doruk, I believe it will be ready shortly! Thanks.

…

On Mon, Aug 3, 2026 at 3:26 PM Doruk Ardağ ***@***.***> wrote: *Dooruk* left a comment (GEOS-ESM/swell#804) <#804 (comment)> @metdyn <https://github.com/metdyn>, while it is not a strict requirement it would be nice to get this in before next week's SWELL documentation sprint. Do you think this is ready to be merged in soon? — Reply to this email directly, view it on GitHub <#804?email_source=notifications&email_token=AIVBWODKTMXTMPTLNWJ6XNL5IDRNJA5CNFSNUABFM5UWIORPF5TWS5BNNB2WEL2JONZXKZKDN5WW2ZLOOQXTKMJXGA3TQNBYGE4KM4TFMFZW63VHNVSW45DJN5XKKZLWMVXHJLDGN5XXIZLSL5RWY2LDNM#issuecomment-5170784818>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AIVBWOGT5ZDOTYOCTKN5YQL5IDRNJAVCNFSNUABFKJSXA33TNF2G64TZHM2DKMRQGI4DGNJTHNEXG43VMU5TINRTGY2TGNBVHA32C5QC> . Triage notifications, keep track of coding agent tasks and review pull requests on the go with GitHub Mobile for iOS <https://github.com/notifications/mobile/ios/AIVBWOBBFSIUPGGG6GGAO5L5IDRNJA5CNFSNUABFM5UWIORPF5TWS5BNNB2WEL2JONZXKZKDN5WW2ZLOOQXTKMJXGA3TQNBYGE4KM4TFMFZW63VHNVSW45DJN5XKKZLWMVXHJKTGN5XXIZLSL5UW64Y> and Android <https://github.com/notifications/mobile/android/AIVBWOEUCQJW2DSFNO5PNF35IDRNJA5CNFSNUABFM5UWIORPF5TWS5BNNB2WEL2JONZXKZKDN5WW2ZLOOQXTKMJXGA3TQNBYGE4KM4TFMFZW63VHNVSW45DJN5XKKZLWMVXHJLTGN5XXIZLSL5QW4ZDSN5UWI>. Download it today! You are receiving this because you were mentioned.Message ID: ***@***.***>

[metdyn]

@Dooruk,
This PR is ready.
(The workflow is tested on discover and has run through EVA plot steps.)

[Dooruk]

@Dooruk, This PR is ready. (The workflow is tested on discover and has run through EVA plot steps.)

Sounds good. @jeromebarre you requested some changes and I don't want to dismiss them without your input. When you get a chance please take a final look so we can merge this.

[Dooruk]

Looks good from my end.

[mranst]

Looks good!

[Dooruk]

Will merge this tomorrow if we don't hear back from Jerome.

[jeromebarre]

I am away right now in a conference. Please go ahead and have this merged. Thanks!

Old